wrote script to find background rate using nls

1 files changed, 32 insertions, 0 deletions

diff --git a/R Scripts/fit-background.R b/R Scripts/fit-background.R
new file mode 100644
index 0000000..d385e7f
--- /dev/null
+++ b/R Scripts/fit-background.R
@@ -0,0 +1,32 @@
+library(igraph)
+setwd("~/Documents/Violence Cascades/")
+load('Results/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, mu, A, phi) {mu + A*(sin((pi/365)*x+phi))^2}
+fit_form = counts ~ mu + A*(sin((pi/365)*days+phi))^2
+
+# explore data
+plot(days,counts)
+curve(fit(x, mu=1, A=12, phi=2.8), add=TRUE ,lwd=2, col="steelblue")
+
+res = nls(formula=fit_form, data=infs, start=list(mu=1, A=12, phi=2.8))
+co = coef(res); co
+
+plot(t)
+plot(days,counts,pch=20,cex=0.4,xlab='Day',ylab='Number of Infections')
+curve(fit(x, mu=co["mu"], A=co["A"], phi=co["phi"]), add=TRUE ,lwd=5, col="steelblue")