1 files changed, 10 insertions, 7 deletions
diff --git a/R Scripts/sim-analysis.R b/R Scripts/sim-analysis.R
index 51659f5..39fd476 100755
--- a/R Scripts/sim-analysis.R
+++ b/R Scripts/sim-analysis.R
@@ -11,6 +11,7 @@ idxs = which(!is.na(vic_days)) %% vcount(lcc)
 vic_days = vic_days[!is.na(vic_days)]
 n.infected = sum(lcc_verts$vic)
 birthyears = as.numeric(format(as.Date(lcc_verts$dob[idxs]),'%Y'))
+districts = as.numeric(lcc_verts$district[idxs])
 
 # find all pairs of vic neighbors to check
 ## nbrs = neighborhood(lcc,order=1,nodes=vic_ids); save(nbrs,file='Raw Data/vic-nbrs.RData')
@@ -33,16 +34,17 @@ n = 10000
 ptm = proc.time()
 simdata = foreach(q=1:n, .combine=rbind) %dopar% {  
   # shuffle infection dates
-  sim.dates = rep(0,length(vic_days))
+  sim.dates = vic_days # data
   for(y in unique(birthyears)){
-    year_ids = which(birthyears==y)
-    if (length(year_ids)>1) sim.dates[year_ids] = sample(vic_days[year_ids])
+    for(d in unique(districts)){
+      ids = which(birthyears==y & districts==d)
+      if (length(ids)>1) sim.dates[ids] = sample(vic_days[ids])
+    }
   }
-#   sim.dates = vic_days # data
   
   vic.dt = rep(0,dim(vic.nbrs)[1])
   for(i in 1:dim(vic.nbrs)[1]){
-    mindt = 5000
+    mindt = 50000
     tus = sim.dates[which(vic.nbrs[i,1]==idxs)]
     tvs = sim.dates[which(vic.nbrs[i,2]==idxs)]
     for(tu in tus){
@@ -73,16 +75,17 @@ xlabs = c('Mean days between infections',
           'Infections within 30 days',
           'Infections within 100 days')
 xlims = matrix(c(660,765,475,625,0.03,0.08,0.1,0.18),ncol=2,byrow=T)
+yl = c(0,0.11*n)
 par(mfrow=c(2,2))
 for(i in 1:4){
   data = d[i] 
   sdata = simdata[,i]
   xl = xlims[i,]
-  h = hist(sdata,50,xlim=xl,col='#1f78b4',freq=T,border=NA,axes=F,
+  h = hist(sdata,50,xlim=xl,ylim=yl,col='#1f78b4',freq=T,border=NA,axes=F,
        xlab=xlabs[i],ylab='Relative frequency',main=NULL)
   if(i %in% c(1,2)) axis(1,at=pretty(xl,5))
   if(i %in% c(3,4)) axis(1,at=pretty(xl,5),lab=paste(pretty(xl,5)*100,'%',sep=''))
-  axis(2,at=pretty(h$counts,3),lab=pretty(h$counts/n,3))
+  axis(2,at=c(0,0.05,0.1)*n,lab=c(0,0.05,0.1))
   abline(v=data,lwd=4,col='#e41a1c')
   box(lwd=1.1)
 }