library(igraph) setwd("~/Documents/Violence Cascades/") load('Raw Data/hyper-lcc.RData') # load data hyp_lcc = upgrade_graph(hyp_lcc) vic_ids = which(V(hyp_lcc)$vic==TRUE) non_vic_ids = which(V(hyp_lcc)$vic==FALSE) hist(as.numeric(hyp_lcc_verts$vic.day[vic_ids]),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(hyp_lcc_verts$vic.day[vic_ids],levels=days)) counts = as.vector(t) infs = data.frame(days,counts) # 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.9,col='#377EB8', xlab='Day',ylab='Number of Infections') curve(fit(x, lambda=co["lambda"], A=co["A"], phi=co["phi"]), add=TRUE ,lwd=5, col="steelblue")