1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
|
library(igraph)
setwd("~/Documents/Violence Cascades/")
load('Results/hyper-lcc.RData')
vic_ids = which(V(hyp_lcc)$vic==1)
n.infections = length(vic_ids)
n.days = max(hyp_lcc_verts$vic.day,na.rm=T) - min(hyp_lcc_verts$vic.day,na.rm=T)
inf.dates = hyp_lcc_verts$vic.day[vic_ids]
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')
# 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.30 = rep(NA,n)
mean.50 = rep(NA,n)
n.vicpairs = rep(NA,n)
hyp_lcc = upgrade_graph(hyp_lcc)
ptm = proc.time()
for(q in 1:n){
if (q%%250==0) print(paste('run:',q))
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
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)
med.time[q] = median(vic.time)
mean.30[q] = mean(vic.time<=30)
mean.50[q] = mean(vic.time<=50)
# decay.rate[sim,q] = getDecayRate(vic.time) #get lm of histogram. or maybe just cor?
}
print(proc.time()-ptm)
# save(mean.time,med.time,mean.30,mean.50,file='Results/sims-Aug28.RData')
####### plot simulation data #######
data = 0.1179039 # 778.2056, 638, 0.05105811, 0.07248908
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
dates = c(V(lcc)$fatal_date[!is.na(V(lcc)$fatal_date)],
V(lcc)$nonfatal_date_1[!is.na(V(lcc)$nonfatal_date_1)],
V(lcc)$nonfatal_date_2[!is.na(V(lcc)$nonfatal_date_2)],
V(lcc)$nonfatal_date_3[!is.na(V(lcc)$nonfatal_date_3)],
V(lcc)$nonfatal_date_4[!is.na(V(lcc)$nonfatal_date_4)],
V(lcc)$nonfatal_date_5[!is.na(V(lcc)$nonfatal_date_5)])
dates = as.Date(dates)
hist(dates, breaks='months',freq=T,col='lightblue',format='%m/%Y',
xlab='Month',ylab='Number of Shootings',main='')
|
