#clear the working environment - this will remove everything!
rm(list=ls())

#load libraries
library(rgdal)
library(foreign)
library(spdep)
library(ggplot2)
library(RColorBrewer)
library(sp)
library(maptools)
library(coda)

#set working directory
setwd("/Users/abennett1/Documents/Teaching/Spatial Epi 2017/")

#nc data set
sids= readShapePoly(system.file("shapes/sids.shp", package="maptools")[1], IDvar="FIPSNO", proj4string=CRS("+proj=longlat +ellps=clrk66"))

##set wd and load New York leukemia data
NY8<-readOGR("./Spatial regression 10_25_17 data/NY_data","NY8_utm18")

TCE<-readOGR("./Spatial regression 10_25_17 data/NY_data","TCE")

cities<-readOGR("./Spatial regression 10_25_17 data/NY_data","NY8cities")

#view data - what variables are included?
head(NY8@data)

#information on variables is here
?NY_data

#coordinates
coords<-coordinates(NY8)

#plot the census tracks and hazardous waste sites
plot(NY8)
plot(TCE,pch=16,cex=1,add=T)

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

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

#view the neighbors
par(mfrow=c(1,2))
plot(NY8)
plot(NY8_nb,coords,col="blue",add=T)
plot(NY8)
plot(NY8_nbr,coords,col="green",add=T)

#highlight the different neighbors
par(mfrow=c(1,1))
plot(NY8)
plot(NY8_nb,coords,col="red",add=T)
plot(NY8_nbr,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 distance between points
dsts<-unlist(nbdists(NY8_nb_Dist, coords))
summary(dsts)

#get the max (or min) distance
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)

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)


##set weights - contiguity
#weights style W - row standardized
NY8_w<-nb2listw(NY8_nb)
NY8_w

#weights style B - binary
NY8_wB<-nb2listw(NY8_nb,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)

#can also use moran.mc or lm.morantest
moran.mc(NY8$Cases,listw=NY8_w,nsim=1000)

#transforming case outcome to linear incidence - why do we want to do this?
#allows us to run linear (normal) models and standardizes to population 
NY8$Zi<-log((1000*(NY8$Cases+1))/NY8$POP8)

#we check spatial autocorrelation of the transformed values
moran.test(NY8$Zi, listw=NY8_w)

#we can use a linear model formulation and test the residuals
lm.morantest(lm(Zi~1,data=NY8),listw=NY8_w)

#adjust for population size - different population sizes can create appearance of spatial autocorrelation 
lm_Zi<-lm(Zi~1,data=NY8,weights=POP8)
lm.morantest(lm_Zi,listw=NY8_w)

##what if we now add some covariates- which factors look to be associated with log incidence?
nylm<-lm(Zi~PEXPOSURE+PCTAGE65P+PCTOWNHOME,data=NY8)
summary(nylm)

#we can test the residuals of the linear model with covariates - what do you find?  Did adding covariates 
#affect residual autocorrelation?
lm.morantest(nylm,NY8_w)

#we can add the residuals to the dataset for plotting
NY8$lmresid<-residuals(nylm)
moran.test(NY8$lmresid,listw=NY8_w)

#what if we adjust for population size - does this change the interpretation?
nylm_w<-lm(Zi~PEXPOSURE+PCTAGE65P+PCTOWNHOME,data=NY8,weights=POP8)
summary(nylm_w)

lm.morantest(nylm_w,NY8_w)




##now let's try a SAR model
nysar<-spautolm(Zi~PEXPOSURE+PCTAGE65P+PCTOWNHOME,data=NY8,listw=NY8_w)
#Check the results - Lambda is the spatial coefficient;  AIC gives you model fit
summary(nysar)

#let's try with distance neighborhood matrix
nysard<-spautolm(Zi~PEXPOSURE+PCTAGE65P+PCTOWNHOME,data=NY8,listw=NY8_dist_w)
summary(nysard)

#with population weights
nysarw<-spautolm(Zi~PEXPOSURE+PCTAGE65P+PCTOWNHOME,data=NY8,listw=NY8_w,weights=POP8)
summary(nysarw)


##let's run similar CAR models - for CAR we have to use symmetric weights (B)
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)

#what do we conclude about covariate associations?

##we can also try a mixed effects model using package nlme
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.ns)
summary(spmodel.s)




#######disease mapping - mapping observed against expected

#create expected
##we are going to model the log relative risk
#first we create the expected number of cases per census tract
NY8$EXP<-NY8$POP8*(sum(NY8$Cases)/sum(NY8$POP8))

#create "SMR" 
NY8$SMR<-NY8$Cases/NY8$EXP

#use spplot to plot SMR or relative risk
spplot(NY8, "SMR")

##plotting using ggplot - we have to use 'fortify' for a spatial polygons data frame
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")+coord_fixed(1)


###let's try some Bayesian models
##CARBayes
NY8$Cases2<-round(NY8$Cases)
W.mat<-nb2mat(NY8_nb,style="B")

#linear model
form1<-Zi~PEXPOSURE+PCTAGE65P+PCTOWNHOME
model.spatial<-S.CARleroux(formula=form1,family="gaussian",data=NY8@data,W=W.mat,burnin=2000,n.sample=10000,thin=10)

#get some outputs - how does this compare to previous?
model.spatial
plot(model.spatial$samples$beta[,2])
plot(model.spatial$samples$beta[,3])

#get the fitted values
NY8$fitted1<-model.spatial$fitted.values

#compare raw and fitted
library(gridExtra)
p1<-spplot(NY8,"Zi")
p2<-spplot(NY8,"fitted1")

grid.arrange(p1,p2,nrow=1)

#let's try a Poisson model - BYM form
form2<-Cases2~PEXPOSURE+PCTAGE65P+PCTOWNHOME+offset(log(EXP))

model.spatial<-S.CARbym(formula=form2,family="poisson",data=NY8@data,W=W.mat,burnin=2000,n.sample=10000,thin=10)

NY8$fitted2<-model.spatial$fitted.values
NY8$RR<-exp(log(NY8$fitted2)-log(NY8$EXP))

#some plotting - how do the smooth rates compare?
p3<-spplot(NY8,"SMR")
p4<-spplot(NY8,"RR")

grid.arrange(p3,p4,nrow=1)



#INLA
row.names(NY8@data)<-seq(1:nrow(NY8@data))
NY8_nb<-poly2nb(NY8)
nb_NY8<-nb2INLA(file="NY8_nb.adj",NY8_nb)

Observed<-NY8$Cases2
PEXPOSURE<-NY8$PEXPOSURE
PCTAGE65P<-NY8$PCTAGE65P
PCTOWNHOME<-NY8$PCTOWNHOME
Expected<-NY8$EXP
ID<-as.numeric(row.names(as(NY8, "data.frame")))
data<-list(Observed=Observed,PEXPOSURE=PEXPOSURE,PCTAGE65P=PCTAGE65P,PCTOWNHOME=PCTOWNHOME,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 + PEXPOSURE + PCTAGE65P + PCTOWNHOME + f(ID,model="bym",graph="NY8_nb.adj") 

model.inla<-inla(formula.bym,family="poisson",data=data,control.compute=list(dic=T,graph=T),
                 E=Expected)

#how does this compare to CARBayes?
summary(model.inla)
model.inla$dic$dic

NY8$fitted3<-model.inla$summary.fitted.values$mean

#how do the fitted compare to CARBayes?
p5<-spplot(NY8,"RR")
p6<-spplot(NY8,"fitted3")

grid.arrange(p5,p6,nrow=1)

#lets compare to a non-spatial model (iid)

formula.iid<-Observed ~ 1 + PEXPOSURE + PCTAGE65P + PCTOWNHOME + f(ID,model="iid",graph="NY8_nb.adj") 

model.inla<-inla(formula.iid,family="poisson",data=data,control.compute=list(dic=T,graph=T),
                 E=Expected)

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


#IN-class and homework
#1. Load the Scottish Lip Cancer dataset

#Observed=number of lip cancer cases
#Expected=expected number of lip cancer cases (population standardized)
#PcAFF=proportion of population working in agriculture, forestry, or fisheries (why might this be important?)

#2. Test for spatial autocorrelation - you should know how to do this

#3. Compare a spatial model to a non-spatial model - which is better and why?

#4. Is the effect of PcAFF the same in both models? What is your interpretation of this covariate
  # in the best model?

#5. Run a spatial model and map the fitted compared to the observed 


##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
O<-spplot(scot_LL,"Observed",main=list(label="Observed"),colorkey=list(at=seq(0,100,10)))
E<-spplot(scot_LL,"Expected",main=list(label="Expected"),colorkey=list(at=seq(0,100,10)))

grid.arrange(O,E,nrow=1)

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

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

#plot in spplot
spplot(scot_LL,"SMR")

#plot in ggplot2
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=SMR)+geom_polygon()+scale_fill_continuous("SMR")