library(rgdal)
library(foreign)
library(spdep)
library(ggplot2)
library(RColorBrewer)
library(classInt)

##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
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<-lme(Zi~PEXPOSURE+PCTAGE65P+PCTOWNHOME,random=~1|AREAKEY,data=as(NY8,"data.frame"),
             correlation=scor,method="ML")
spmodel<-lme(Zi~PEXPOSURE+PCTAGE65P+PCTOWNHOME,random=~1|AREAKEY,data=as(NY8,"data.frame"),
             method="ML")
summary(spmodel)


##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)

#output summary
res$summary
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

##ggplot map
require('plyr')
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=prob)+geom_polygon()+scale_fill_continuous("probability")

##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

prob_plot <- classIntervals(NY8$prob, style = "fixed", fixedBreaks = brks1)
theta_plot <- classIntervals(NY8$theta, style = "fixed", fixedBreaks = brks2)
prob_colcode <- findColours(prob_plot, plotclr)
theta_colcode <- findColours(theta_plot, plotclr)

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

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

#densities
theta_dens<-res$sims.array[,,"theta"][,1]
beta1_dens<-res$sims.array[,,"beta[1]"][,1]
plot(density(theta_dens))


##choose the best model for the Scottish Lip Cancer dataset
##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)

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
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
res2$DIC