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
|
library(igraph)
setwd("~/Documents/Violence Cascades/")
# load data
load('Raw Data/lcc.RData')
vic_days = as.numeric(unlist(lcc_verts[,16:21]))
vic_days = vic_days[!is.na(vic_days)]
hist(vic_days,2000,col='lightblue',
xlab='Day of Study Period',main='Infections During the Study Period')
# get infection counts per day
days = 1:3012
t = table(factor(vic_days,levels=days))
counts = as.vector(t)
start_date = as.Date("2005-12-31")
dates = start_date + days
infs = data.frame(days,counts,dates)
# define background function
fit = function(x, lambda, A, phi) {lambda + A*(sin((2*pi/365.24)*x+phi))}
fit_form = counts ~ lambda + A*(sin((2*pi/365.24)*days+phi))
# explore data
plot(days,counts)
curve(fit(x, lambda=3, A=2, phi=4), add=TRUE ,lwd=4, col="steelblue")
res = nls(formula=fit_form, data=infs, start=list(lambda=3, A=2, phi=4))
co = coef(res); co
plot(t)
plot(days,counts,pch=20,cex=0.5,col='#1f78b4',xaxt='n',
xlab='',ylab='Number of shootings',cex.axis=0.6,cex.lab=0.6)
curve(fit(x, lambda=co["lambda"], A=co["A"], phi=co["phi"]),
add=TRUE ,lwd=3, col="#1a9850")
axis(1,at=seq(1,max(days),365),lab=2006:2014,cex.axis=.6)
|
