###session 9 - space-time data
library(maptools)
library(INLA)
library(spdep)
library(foreign)
library(sp)
library(ggplot2)
library(lattice)
library(zoo)

#read in shape file for districts in Zambia
geoZAM<-readShapePoly("/Users/abennett1/Documents/Teaching/Spatial Epi 2017/Nov 15 17/data/admin2.shp")
plot(geoZAM)

setwd("/Users/abennett1/Documents/Teaching/Spatial Epi 2017/Nov 15 17/data/")

#read in data in district-time long format - data are malaria cases per district and month from 
#Jan 2009 - Dec 2011... covariates are standardized and include temperature, rainfall, vegetation, bed nets
zaminla<-read.table(file="newHMIS3_27.txt",sep=",",header=TRUE)
zam.df<-as.data.frame(zaminla)

adjpoly<-poly2nb(geoZAM)
neiZM<-nb2WB(adjpoly)

nb2INLA(file="ZM.adj",adjpoly)

#add log E for offset 
zam.df$logE<-log(zam.df$E)
zam.df$monthid2<-zam.df$monthid

#run an inla model with space and time effects #we use a random walk prior - this could also be 'ar1' 
#see inla latent models here for more info: http://www.r-inla.org/models/latent-models 
formula.ST<-cases ~ 1 + offset(logE) + itnperhhanom + as.factor(evim1cat) + f(id,model="bym",graph="ZM.adj") + f(monthid,model="rw1")
model.inla.ST<-inla(formula.ST,family="nbinomial",data=zam.df,control.compute=list(dic=T,graph=T))

#check the fixed results
summary(model.inla.ST)
exp(model.inla.ST$summary.fixed)
model.inla.ST$dic$dic

#look at the random effects by area
sum<-model.inla.ST$summary.random
RE1<-sum$id[1:length(unique(zam.df$id)),2]        #uncorrelated
RE2<-sum$id[(length(unique(zam.df$id))+1):(2*(length(unique(zam.df$id)))),2]    #correlated

#add the random effects to the spdf
geoZAM$RE1<-RE1
geoZAM$RE2<-RE2

#plot spatial effects using spplot - you can also try ggplot
library(gridExtra)
p1<-spplot(geoZAM,"RE1")
p2<-spplot(geoZAM,"RE2")

grid.arrange(p1,p2,nrow=1)   #doesn't appear to be any differences - suggests limited spatial effect at this scale

#plot temporal effect
plot(sum$monthid[1:36,1],sum$monthid[1:36,2],type="l")
points(sum$monthid[1:36,4],type="l",lty=2)
points(sum$monthid[1:36,6],type="l",lty=2)

#add an interaction effect
zam.df$ID.area.month <- seq(1,length(unique(zam.df$id)))

#run a model with the interaction + an additional time term
formula.STint<- cases ~ 1 + itnperhhanom + as.factor(evim1cat) + f(id,model="bym",graph="ZM.adj") +
  f(monthid,model="rw1") + f(monthid2,model="iid") +f(ID.area.month, model="iid")
model.inla.STint <- inla(formula.STint,family="nbinomial",data=zam.df, E=E, verbose=TRUE,control.compute=list(dic=T,graph=T))

#check the fixed results - does the model improve? does the interpretation change at all?
summary(model.inla.STint)
exp(model.inla.STint$summary.fixed)
model.inla.STint$dic$dic

#check the random results
sum<-model.inla.STint$summary.random  
RE1int<-sum$id[1:length(unique(zam.df$id)),2]           #uncorrelated
RE2int<-sum$id[(length(unique(zam.df$id))+1):(2*(length(unique(zam.df$id)))),2]     #correlated

#add random effects to teh spdf
geoZAM$RE1int<-RE1int
geoZAM$RE2int<-RE2int

#plot random effects 
library(gridExtra)
p1<-spplot(geoZAM,"RE1int")
p2<-spplot(geoZAM,"RE2int")

grid.arrange(p1,p2,nrow=1)  #no change

#plot temporal effect
plot(sum$monthid[1:36,1],sum$monthid[1:36,2],type="l")
points(sum$monthid[1:36,4],type="l",lty=2)
points(sum$monthid[1:36,6],type="l",lty=2)



###############################################################################
##using facility tesselations###################################################

setwd("/Users/abennett1/Documents/Teaching/Spatial Epi 2017/Nov 15 17/data/")

#add facility tesselation - polygon created for every health facility point
geoZAM2<-readShapePoly("zamvoronoiclip3sjoin2.shp")

#plot facility tesselation
plot(geoZAM2)

#prepare inla objects
adjpoly2<-poly2nb(geoZAM2)

nb2INLA(file="ZMhf.adj",adjpoly2)

#add facility level data - case data by facility and month for 2009-2012
zam.hf<-read.table(file="zamHF09_12confimp.txt",sep=",",header=TRUE)
zam.hf.df<-as.data.frame(zam.hf)

#create an interaction term
zam.hf.df$monthid2<-zam.hf.df$monthid
zam.hf.df$ID.area.month <- seq(1,length(unique(zam.hf.df$id)))

#create expected counts for each facility based on mean over the entire period
A<-c()
for (i in 1:1369) {
  A[i]<-as.numeric(with(subset(zam.hf.df, id==i), mean(totalconfcaseST, na.rm=TRUE)))
}

E<-c()
for (i in 1:1369) {
  E<-c(E, rep(A[i], times=48, na.rm=FALSE))
}

#add expected counts
zam.hf.df$E<-E

#subset data - for speed
zam.hf.sub.df<-subset(zam.hf.df, monthid>24)

#run interaction model - takes some time
formula.ST2int<- totalconfcaseST ~ 1 + f(id,model="bym",graph="ZMhf.adj") +
  f(monthid,model="rw1") + f(monthid2,model="iid") +f(ID.area.month, model="iid")
model.inla.ST2int <- inla(formula.ST2int,family="poisson",data=zam.hf.sub.df, E=E, verbose=TRUE)

#check the random effects
sum<-model.inla.ST2int$summary.random
RE1<-sum$id[1:length(unique(zam.hf.sub.df$id)),2]   # uncorrelated RE
RE2<-sum$id[1370:2738,2]   # correlated RE

#add to the hf spdf
geoZAM2$RE1<-RE1
geoZAM2$RE2<-RE2

#plot using ggplot
require('plyr')
geoZAM2@data$id<-rownames(geoZAM2@data)
geoZAM2.points<-fortify(geoZAM2,region="id")
geoZAM2.df<-join(geoZAM2.points,geoZAM2@data,by="id")
g1<-ggplot(geoZAM2.df)+aes(long,lat,group=group,fill=RE1)+geom_polygon()+scale_fill_continuous("UH")+coord_fixed()
g2<-ggplot(geoZAM2.df)+aes(long,lat,group=group,fill=RE2)+geom_polygon()+scale_fill_continuous("CH")+coord_fixed()

grid.arrange(g1,g2,nrow=1)  #compare the UH and CH




#Plotting space time interaction - exceedence probabilities for rr>2
#duplicate spatial data frame
geoZAM3<-geoZAM2

delta <- data.frame(delta=sum$ID.area.month[,2],month=zam.hf.sub.df$monthid,ID.area=zam.hf.sub.df$id)
delta.matrix <- matrix(delta[,1], 1369 ,24,byrow=FALSE)
rownames(delta.matrix)<- delta[1:1369,3]

#Space time probability>2
a=0.69314
inlaprob.delta<-lapply(model.inla.ST2int$marginals.random[[4]], function(X){
  1-inla.pmarginal(a, X)
})
pp.delta<-unlist(inlaprob.delta)

pp.cutoff.interaction <- c(0,0.5,0.9999,1)
pp.delta.matrix <- matrix(pp.delta, 1369,24,byrow=TRUE)
#pp.delta.matrix <- pp.delta.matrix[1:1369,]
pp.delta.factor <- data.frame(NAME=seq(1,1369))
for(i in 1:24){
  pp.delta.factor.temp <- cut(pp.delta.matrix[,i],breaks=pp.cutoff.interaction,include.lowest=TRUE) 
  pp.delta.factor <- cbind(pp.delta.factor,pp.delta.factor.temp)
}
colnames(pp.delta.factor)<- c("NAME",seq(25,48))

#Maps - plotting exceedence probabilities 
attr(geoZAM3, "data")=data.frame(pp.delta.factor)

trellis.par.set(axis.line=list(col=NA))
spplot(obj=geoZAM3, zcol="X39", col.regions=c(gray(1.5:0.5/2),2),main="Month 39",par.settings=list(fontsize=list(text=17)))


#plotting with space-time 
library(spacetime)
library("maps")
library("plm")
##data("Produc")

#time variable
monthid=1:length(unique(zam.hf.df$monthid))
time=yearmon(2009 + (monthid-1)/12)

#space-time data frame
spacetime=STFDF(geoZAM2,time,zam.hf.df[order(zam.hf.df[1],zam.hf.df[2]),])

library(RColorBrewer)
stplot(spacetime[,1:12,"totalconfcaseST"],names.attr=monthid[1:12],col.regions=brewer.pal(9,"YlOrRd"),cuts=9)



##ST-Plots for district

#time variable
monthid=1:length(unique(zam.df$monthid))
time=yearmon(2009 + (monthid-1)/12)

#space-time data frame - needs to be ordered with spatial index moving fastest
spacetime=STFDF(geoZAM,time,zam.df[order(zam.df[1],zam.df[2]),])

#plotting cases 
stplot(spacetime[,1:12,"cases"],names.attr=monthid[1:12],col.regions=brewer.pal(9,"YlOrRd"),cuts=9)
title("Monthly cases")

#plotting minimum temperature
stplot(spacetime[,1:12,"mintm2cstan"],names.attr=c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec"),col.regions=brewer.pal(9,"YlOrRd"),cuts=9)
title("Standardized min temp")