diff options
| -rwxr-xr-x | R Scripts/sim-analysis.R | 61 |
1 files changed, 33 insertions, 28 deletions
diff --git a/R Scripts/sim-analysis.R b/R Scripts/sim-analysis.R index 5010485..e368419 100755 --- a/R Scripts/sim-analysis.R +++ b/R Scripts/sim-analysis.R @@ -12,12 +12,28 @@ birthyears = as.numeric(format(as.Date(hyp_lcc_verts$dob[vic_ids]),'%Y')) # nbrs = neighborhood(graph, nodes=vic_ids, order=1) load('Results/vic-nbrs.RData') -n = 150 -mean.time = matrix(0,1,n) -med.time = matrix(0,1,n) -mean.50 = matrix(0,1,n) -mean.100 = matrix(0,1,n) -n.vicpairs = matrix(0,1,n) +# find all pairs of vic neighbors to check +vic.nbrs = matrix(NA,14885,2) +idx = 1 +for (i in 1:n.infections){ + u = vic_ids[i] + nbhd = nbrs[[i]] + nbhd = intersect(vic_ids,nbhd) + nbhd = nbhd[u<nbhd] # only looks at pairs with u<v to avoid double-counting + nbhd = setdiff(nbhd,u) + nbhd = nbhd[hyp_lcc_verts$ir_no[u] != hyp_lcc_verts$ir_no[nbhd]] # don't count infections of the same person + if (length(nbhd)==0) next + vic.nbrs[idx:(idx+length(nbhd)-1),] = cbind(rep(match(u,vic_ids),length(nbhd)),match(nbhd,vic_ids)) + idx = idx + length(nbhd) +} + +# run simulation +n = 10000 +mean.time = rep(NA,n) +med.time = rep(NA,n) +mean.50 = rep(NA,n) +mean.100 = rep(NA,n) +n.vicpairs = rep(NA,n) hyp_lcc = upgrade_graph(hyp_lcc) @@ -27,28 +43,15 @@ for(q in 1:n){ graph = hyp_lcc vics = vic_ids + # shuffle infection dates sim.dates = rep(0,length(vic_ids)) for(y in unique(birthyears)){ year_ids = which(birthyears==y) if (length(year_ids)>1) sim.dates[year_ids] = sample(inf.dates[year_ids]) } -# {vics = vic_ids; sim.dates = inf.dates} # data +# vics = vic_ids; sim.dates = inf.dates # data - vic.time = rep(NA,14885) - idx = 1 - for (i in 1:n.infections){ - u = vics[i] - nbhd = nbrs[[i]] - nbhd = intersect(vic_ids,nbhd) - nbhd = nbhd[u<nbhd] # only looks at pairs with u<v to avoid double-counting - nbhd = setdiff(nbhd,u) - nbhd = nbhd[hyp_lcc_verts$ir_no[u] != hyp_lcc_verts$ir_no[nbhd]] # don't count infections of the same person - if (length(nbhd)==0) next - tu = sim.dates[i] - tvs = sim.dates[match(nbhd,vic_ids)] - vic.time[idx:(idx+length(nbhd)-1)] = abs(tu-tvs) - idx = idx + length(nbhd) - } + vic.time = abs(sim.dates[vic.nbrs[,1]] - sim.dates[vic.nbrs[,2]]) n.vicpairs[q] = length(vic.time) mean.time[q] = mean(vic.time) @@ -59,15 +62,17 @@ for(q in 1:n){ } print(proc.time()-ptm) +# save(mean.time,med.time,mean.50,mean.100,file='Results/sims-Aug28.RData') ####### plot simulation data ####### -data = 634 # 774.5483, 634, 0.07252472, 0.1176774 -simdata = med.time -xl = c(0.9*min(min(simdata),data),1.1*max(max(simdata),data)) -yl = c(0,max(hist(simdata,15)$density)) -hist(simdata,15,xlim=xl,ylim=yl,col=rgb(0,1,0,1/4),freq=F, - xlab='Neighbors Infected Within 100 Days', main=NULL) +data = 0.1179039 # 778.2056, 638, 0.07248908, 0.1179039 +simdata = mean.100 +xl = c(0.95*min(min(simdata),data),1.05*max(max(simdata),data)) +yl = c(0,max(hist(simdata,25)$density)) +hist(simdata,25,xlim=xl,ylim=yl,col=rgb(0,1,0,1/4),freq=F, + xlab='Neighbor Infections Within 100 Days', main=NULL) abline(v=data,lwd=4,col='red') +box() ##### inf dates histogram |