library(rgdal)
library(foreign)
library(spdep)
library(ggplot2)
library(RColorBrewer)
library(classInt)
library(sp)
library(maptools)
require('RANN')
library(coda)
##set wd and load New York leukemia data
setwd("C:\\Users\\BennettA2\\Documents\\Teaching\\Spatial Epi\\Session7")
NY8<-readOGR("NY_data","NY8_utm18")

TCE<-readOGR("NY_data","TCE")

#cities<-readOGR("NY_data","NY8cities")

##load existing neighborhood matrix
#NY_nb<-read.gal("NY_data\\NY_nb.gal",region.id=rownames(as(NY8,"data.frame")))

coords<-coordinates(NY8)
#plot(NY_nb,coords)

# Contiguity neighbors - all that share a boundary point
NY8_1 <- poly2nb(NY8)  #queen

NY8_2 <- poly2nb(NY8,queen=F)  #rook

par(mfrow=c(1,2))
plot(NY8)
plot(NY8_1,coords,col="blue",add=T)
plot(NY8)
plot(NY8_2,coords,col="green",add=T)

par(mfrow=c(1,1))
plot(NY8)
plot(NY8_1,coords,col="red",add=T)
plot(NY8_2,coords,add=T)

# Distance-based neighbors - returns the k nearest points as neighbors
IDs<-row.names(as(NY8, "data.frame"))
NY8_nb_Dist <- knn2nb(knearneigh(coords, k=1),row.names=IDs)
par(mfrow=c(1,1))
plot(NY8)
plot(NY8_nb_Dist, coords, col="green",add=T)

NY8_nb_Dist2<- knn2nb(knearneigh(coords, k=2),row.names=IDs)
NY8_nb_Dist3<- knn2nb(knearneigh(coords, k=3),row.names=IDs)

par(mfrow=c(1,3))
plot(NY8)
plot(NY8_nb_Dist, coords,add=T)
plot(NY8)
plot(NY8_nb_Dist2, coords,add=T)
plot(NY8)
plot(NY8_nb_Dist3, coords,add=T)

# Find the minimum distance at which all areas have a distance-based neighbor
dsts<-unlist(nbdists(NY8_nb_Dist, coords))
summary(dsts)

max_1nn<-max(dsts)
max_1nn

# use dnearneigh to find neighbors with an interpoint distance, with distances settings
# upper and lower distance bounds
NY8_nb0.5 <- dnearneigh(coords, d1=0, d2=0.5*max_1nn, row.names=IDs)
NY8_nb1 <- dnearneigh(coords, d1=0, d2=1*max_1nn, row.names=IDs)
NY8_nb1.5 <- dnearneigh(coords, d1=0, d2=1.5*max_1nn, row.names=IDs)

nb<-list(d1=NY8_nb0.5, d2=NY8_nb1, d3=NY8_nb1.5)

sapply(nb, function(x) is.symmetric.nb(x, verbose=F, force=T))
sapply(nb, function(x) n.comp.nb(x)$nc)

par(mfrow=c(1,3))
plot(NY8)
plot(NY8_nb0.5, coords, col="green",add=T)
plot(NY8)
plot(NY8_nb1, coords, col="green",add=T)
plot(NY8)
plot(NY8_nb1.5, coords, col="green",add=T)

##plot the number of neighbors by distance type
par(mfrow=c(1,1))

dS <- c(0.5, 1, 1.5)*max_1nn
res <- sapply(nb, function(x) table(card(x)))
mx <- max(card(NY8_nb1.5))
res1 <- matrix(0, ncol=(mx+1), nrow=3)
rownames(res1) <- names(res)
colnames(res1) <- as.character(0:mx)
res1[1, names(res$d1)] <- res$d1
res1[2, names(res$d2)] <- res$d2
res1[3, names(res$d3)] <- res$d3
library(RColorBrewer)
pal <- grey.colors(3, 0.95, 0.55, 2.2)
barplot(res1, col=pal, beside=TRUE, legend.text=FALSE, xlab="numbers of neighbours", ylab="tracts")
legend("topright", legend=format(dS, digits=1), fill=pal, bty="n", cex=0.8, title="max. distance")

##set weights - contiguity
NY8_w<-nb2listw(NY8_1)
NY8_w

NY8_wB<-nb2listw(NY8_1,style="B")
NY8_wB

##set weights - distance: add zero polic=T if there is a 'no neighbor' error
NY8_dist_w<-nb2listw(NY8_nb1)
NY8_dist_w

##moran's tests of global spatial autocorrelation
moran.test(NY8$Cases,listw=NY8_w)
moran.test(NY8$Cases,listw=NY8_dist_w)

NY8$r<-NY8$Cases/NY8$POP8
NY8$Zi<-log((1000*(NY8$Cases+1))/NY8$POP8)

moran.test(NY8$Zi, listw=NY8_w)
lm.morantest(lm(Zi~1,data=NY8),listw=NY8_w)

lm_Zi<-lm(Zi~1,data=NY8,weights=POP8)
lm.morantest(lm_Zi,listw=NY8_w)

##adding covariates
nylm<-lm(Zi~PEXPOSURE+PCTAGE65P+PCTOWNHOME,data=NY8)
nylm<-lm(Zi~PEXPOSURE+PCTAGE65P+PCTOWNHOME,data=NY8,weights=POP8)
summary(nylm)

NY8$lmresid<-residuals(nylm)
lm.morantest(nylm,NY8_w)

##SAR model
nysar<-spautolm(Zi~PEXPOSURE+PCTAGE65P+PCTOWNHOME,data=NY8,listw=NY8_w)
summary(nysar)
nysard<-spautolm(Zi~PEXPOSURE+PCTAGE65P+PCTOWNHOME,data=NY8,listw=NY8_dist_w)
summary(nysard)
nysarw<-spautolm(Zi~PEXPOSURE+PCTAGE65P+PCTOWNHOME,data=NY8,listw=NY8_w,weights=POP8)
summary(nysarw)


##CAR model - not used
##nycar<-spautolm(Zi~PEXPOSURE+PCTAGE65P+PCTOWNHOME,data=NY8,listw=NY8_wB,family="CAR")
##summary(nycar)
##nycarw<-spautolm(Zi~PEXPOSURE+PCTAGE65P+PCTOWNHOME,data=NY8,listw=NY8_wB,weights=POP8,family="CAR")
##summary(nycarw)


##mixed effects model
library(nlme)
NY8$x<-coordinates(NY8)[,1]/1000
NY8$y<-coordinates(NY8)[,2]/1000
sp1<-corSpatial(1,form=~x+y,type="gaussian")
scor<-Initialize(sp1,as(NY8,"data.frame")[,c("x","y")],nugget=FALSE)

spmodel.ns<-lme(Zi~PEXPOSURE+PCTAGE65P+PCTOWNHOME,random=~1|AREAKEY,data=as(NY8,"data.frame"),
             method="ML")
spmodel.s<-lme(Zi~PEXPOSURE+PCTAGE65P+PCTOWNHOME,random=~1|AREAKEY,data=as(NY8,"data.frame"),
             correlation=scor,method="ML")
summary(spmodel)
summary(spmodel.s)

##Bayesian models R2WinBUGS
library(R2WinBUGS)

#create expected
##we are going to model the log relative risk
NY8$EXP<-NY8$POP8*(sum(NY8$Cases)/sum(NY8$POP8))
NY_nb<-poly2nb(NY8)
NY_nb_wb<-nb2WB(NY_nb)

d<-c(list(O=NY8$Cases,E=NY8$EXP),N=281,list(PCTAGE65P=NY8$PCTAGE65P,PCTOWNHOME=NY8$PCTOWNHOME,
                                            PEXPOSURE=NY8$PEXPOSURE))

#bugs directory
BUGSDir<-"C:\\Users\\BennettA2\\Documents\\Teaching\\Spatial Epi\\winbugs14\\WinBUGS14"

#model file and output directory
mfile<-paste(getwd(),"/NY_BYM.txt",sep="",collapse="")
tdir<-paste(getwd(),"/NYoutput2",sep="",collapse="")
dir.create(tdir)

#initial values
inits1<-list(alpha=1,beta=c(0,0,0),u=rep(0,281),
             v=rep(0,281),precu=1,precv=1)
inits2<-list(alpha=10,beta=c(1,1,1),u=rep(1,281),
             v=rep(1,281),precu=0.1,precv=0.1)

#model call
res<-bugs(data=c(d,NY_nb_wb),inits=list(inits1,inits2),
          parameters.to.save=c("sigmau","sigmav","beta","prob","theta","u","v"),
          model.file=mfile,working.directory=tdir,n.thin=3,n.chains=2,
          bugs.directory=BUGSDir,n.iter=6000,n.burnin=3000,debug=T)

#useWINE=TRUE,
#WINE="/Applications/Wine.app/Contents/Resources/bin/wine",
#WINEPATH="/Applications/Wine.app/Contents/Resources/bin/winepath"

b1<-res$sims.array[,,"beta[1]"][,1]
b2<-res$sims.array[,,"beta[1]"][,2]

plot(b1,type="l",col="red")
lines(b2,type="l",col="blue")

plot(density(b1),lwd=2,col="red")
lines(density(b2),lwd=2,col="blue")             
             
#output summary
res$summary[c("beta[1]","beta[2]","beta[3]"),]
res$summary[c("sigmau","sigmav"),]
res$summary

#DIC
res$DIC


#add outputs to datafile
NY8$prob<-res$mean$prob
NY8$theta<-res$mean$theta
NY8$u<-res$mean$u
NY8$v<-res$mean$v
NY8$SMR<-NY8$Cases/NY8$EXP

##ggplot map
require('plyr')
library(ggmap)
NY8@data$id<-rownames(NY8@data)
NY8.points<-fortify(NY8,region="id")
NY8.df<-join(NY8.points,NY8@data,by="id")
ggplot(NY8.df)+aes(long,lat,group=group,fill=SMR)+geom_polygon()+scale_fill_continuous("Relative risk")


##alternative map
nclr <- 6
plotclr <- brewer.pal(nclr, "BuPu")
nsteps <- 5
step1<-(max(NY8$prob)-min(NY8$prob))/nsteps
brks1<-min(NY8$prob)+(0:nsteps)*step1
brks1[1]<-0

step2<-(max(NY8$theta)-min(NY8$theta))/nsteps
brks2<-min(NY8$theta)+(0:nsteps)*step2
brks2[1]<-0

step3<-(max(NY8$v)-min(NY8$v))/nsteps
brks3<-min(NY8$v)+(0:nsteps)*step3

step4<-(max(NY8$u)-min(NY8$u))/nsteps
brks4<-min(NY8$u)+(0:nsteps)*step4

prob_plot <- classIntervals(NY8$prob, style = "fixed", fixedBreaks = brks1)
theta_plot <- classIntervals(NY8$theta, style = "fixed", fixedBreaks = brks2)
v_plot<-classIntervals(NY8$v, style = "fixed", fixedBreaks = brks3)
u_plot<-classIntervals(NY8$u, style = "fixed", fixedBreaks = brks4)
prob_colcode <- findColours(prob_plot, plotclr)
theta_colcode <- findColours(theta_plot, plotclr)
v_colcode <- findColours(v_plot, plotclr)
u_colcode <- findColours(u_plot, plotclr)

grps1<-(as.ordered(cut(NY8$prob, brks1, include.lowest=TRUE)))

par(mfrow=c(1,2))
plot(NY8, col=prob_colcode, main="Probability RR>1")
plot(NY8, col=theta_colcode, main="Relative Risk")

par(mfrow=c(1,2))
plot(NY8, col=v_colcode, main="Spatial RE")
plot(NY8, col=u_colcode, main="Uncorrelated RE")



##choose the best model for the Scottish Lip Cancer dataset
##does it change your interpretation of the main effect?

##predicting the log relative risk
Scot<-readShapePoly("scot")

scot_dat <- read.csv("scotland.csv", header=T, , stringsAsFactors=FALSE)

names(scot_dat) <- c("District","Observed", "Expected", "PcAFF",
                     "Latitude", "Longitude")

row.names(scot_dat) <- formatC(scot_dat$District, width = 2,
                               flag = "0")
ID <- formatC(Scot$ID, width = 2, flag = "0")
scot_LL <- spChFIDs(Scot, ID)
scot_LL <- spCbind(scot_LL, scot_dat[match(ID, row.names(scot_dat)), ])
plot(scot_LL)

# Look at maps of observed and expected counts
nclr <- 6
plotclr <- brewer.pal(nclr, "BuPu")
brks <- c(0, 5, 10, 25, 50, 100)

O_CI <- classIntervals(scot_LL$Observed, style = "fixed", fixedBreaks = brks)
E_CI <- classIntervals(scot_LL$Expected, style = "fixed", fixedBreaks = brks)
O_colcode <- findColours(O_CI, plotclr)
E_colcode <- findColours(E_CI, plotclr)

par(mfrow=c(1,2))
plot(scot_LL, col=O_colcode, main="Observed")
plot(scot_LL, col=E_colcode, main="Expected")

#create SMR
scot_LL$SMR<-scot_LL$Observed/scot_LL$Expected

Scot_nb<-poly2nb(scot_LL)
Scot_nb_wb<-nb2WB(Scot_nb)

d<-c(list(O=scot_LL$Observed,E=scot_LL$Expected),N=56,PcAFF=list(scot_LL$PcAFF))

#model file and output directory- correlated heterogeneity (spatial model)
mfile<-paste(getwd(),"/Scot_BYM_CH.txt",sep="",collapse="")
tdir<-paste(getwd(),"/Scotoutput2",sep="",collapse="")
dir.create(tdir)

#initial values
inits1<-list(alpha=1,beta1=0,u=rep(0,56),
             v=rep(0,56),precu=1,precv=1)
inits2<-list(alpha=10,beta1=1,u=rep(1,56),
             v=rep(1,56),precu=0.1,precv=0.1)

#model call
res<-bugs(data=c(d,Scot_nb_wb),inits=list(inits1,inits2),
          parameters.to.save=c("sigmau","sigmav","beta1","prob","theta","u","v"),
          model.file=mfile,working.directory=tdir,n.thin=3,n.chains=2,
          bugs.directory=BUGSDir,n.iter=6000,n.burnin=3000,debug=T)

res$summary
scotCH<-res$DIC

mfile<-paste(getwd(),"/Scot_BYM_UH.txt",sep="",collapse="")
tdir<-paste(getwd(),"/Scotoutput3",sep="",collapse="")
dir.create(tdir)

#initial values
inits1<-list(alpha=1,beta1=0,u=rep(0,56),
             precu=1)
inits2<-list(alpha=10,beta1=1,u=rep(1,56),
             precu=0.1)

#model call (uncorrelated heterogeneity)
res2<-bugs(data=d,inits=list(inits1,inits2),
          parameters.to.save=c("sigmau","beta1","prob","theta","u"),
          model.file=mfile,working.directory=tdir,n.thin=3,n.chains=2,
          bugs.directory=BUGSDir,n.iter=6000,n.burnin=3000,debug=T)

res2$summary
scotUH<-res2$DIC

#add outputs to datafile
scot_dat$prob<-res$mean$prob
scot_dat$theta<-res$mean$theta
scot_dat$u<-res$mean$u
scot_dat$v<-res$mean$v

row.names(scot_dat) <- formatC(scot_dat$District, width = 2,
                               flag = "0")
ID <- formatC(Scot$ID, width = 2, flag = "0")
scot_LL <- spChFIDs(Scot, ID)
scot_LL <- spCbind(scot_LL, scot_dat[match(ID, row.names(scot_dat)), ])

O_CI <- classIntervals(scot_LL$theta, style = "sd")
E_CI <- classIntervals(scot_LL$prob, style= "sd")
SMR_CI <- classIntervals(scot_LL$SMR, style= "sd")
O_colcode <- findColours(O_CI, plotclr)
E_colcode <- findColours(E_CI, plotclr)
SMR_colcode <- findColours(SMR_CI, plotclr)

par(mfrow=c(1,2))
plot(scot_LL, col=O_colcode, main="Relative risk")
plot(scot_LL, col=E_colcode, main="Prob  RR>1")

par(mfrow=c(1,2))
plot(scot_LL, col=O_colcode, main="Relative risk")
plot(scot_LL, col=SMR_colcode, main="SMR")

require('plyr')
scot_LL@data$id<-rownames(scot_LL@data)
scot.points<-fortify(scot_LL,region="id")
scot.df<-join(scot.points,scot_LL@data,by="id")
ggplot(scot.df)+aes(long,lat,group=group,fill=theta)+geom_polygon()+scale_fill_continuous("RR")


##INLA MODELS

library(spdep)
library(INLA)
setwd("C:\\Users\\BennettA2\\Documents\\Teaching\\Spatial Epi\\Session7")

Scot<-readShapePoly("scot")

scot_dat <- read.csv("scotland.csv", header=T, , stringsAsFactors=FALSE)

names(scot_dat) <- c("District","Observed", "Expected", "PcAFF",
                     "Latitude", "Longitude")

row.names(scot_dat) <- formatC(scot_dat$District, width = 2,
                               flag = "0")
ID <- formatC(Scot$ID, width = 2, flag = "0")
scot_LL <- spChFIDs(Scot, ID)
scot_LL <- spCbind(scot_LL, scot_dat[match(ID, row.names(scot_dat)), ])
plot(scot_LL)

source("fillmap.R")
maintitle1<-paste("Observed")
fillmap(scot_LL,maintitle1,scot_LL$Observed,n.col=6)

Scot_nb<-poly2nb(scot_LL)
nb_scot<-nb2INLA(file="Scot_nb.adj",Scot_nb)

Observed<-scot_LL$Observed
PcAFF<-scot_LL$PcAFF
Expected<-scot_LL$Expected
ID<-scot_LL$ID
data<-list(Observed=Observed,PcAFF=PcAFF,Expected=Expected,ID=ID)

#which model is best?  
#bym model  = UH + CH
#besag model = correlated heterogeneity (spatial model)
#iid model = UH

formula.bym<-Observed ~ 1 + PcAFF + f(ID,model="iid",graph="Scot_nb.adj") 
model.inla<-inla(formula.bym,family="poisson",data=data,control.compute=list(dic=T,graph=T),
                 E=Expected)

summary(model.inla)
model.inla$dic$dic



sum<-model.inla$summary.random
RE1<-sum$ID[1:56,2]   # uncorrelated RE
RE2<-sum$ID[57:112,2]   #correlated RE


#add outputs to datafile
scot_dat$UH<-RE1
scot_dat$CH<-RE2

row.names(scot_dat) <- formatC(scot_dat$District, width = 2,
                               flag = "0")
ID <- formatC(Scot$ID, width = 2, flag = "0")
scot_LL <- spChFIDs(Scot, ID)
scot_LL <- spCbind(scot_LL, scot_dat[match(ID, row.names(scot_dat)), ])

require('plyr')
scot_LL@data$id<-rownames(scot_LL@data)
scot.points<-fortify(scot_LL,region="id")
scot.df<-join(scot.points,scot_LL@data,by="id")
ggplot(scot.df)+aes(long,lat,group=group,fill=CH)+geom_polygon()+scale_fill_continuous("CH")