#session 8 - logistic regression
rm(list=ls())

library(raster)
library(mgcv)
library(sp)
library(spdep) 
library(ape) 
library(CARBayes)
library(INLA)

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

#load data - malaria data in Burkina Faso
bf_data<-read.table("data_bf2_binomial.txt",sep=",",header=T)
bf_data<-data.frame(bf_data)

#add admin layers 
BF_admin1<-raster::getData(name="GADM",country='BFA',level=1,download=T)

#create a prevalence variable
bf_data$prev<-bf_data$positives/bf_data$examined

#subset out lat/long
bf_data.xy<-bf_data[c("longitude","latitude")]
coordinates(bf_data.xy)<-~longitude+latitude

# first we'll check the relationship with a covariate
# run a glm model with altitude and check residuals
outcome1<-cbind(bf_data$positives,bf_data$examined-bf_data$positives)
m.0<-glm(outcome1~alt,data=bf_data,family="binomial")

print(summary(m.0))  #does it look like there's a relationship?
confint(m.0)
exp(coef(m.0))    #for reporting odds ratios
exp(confint(m.0))   #and confidence intervals

#get the residuals (these are deviance residuals)
resid<-residuals(m.0)

#we can plot them first - do they look autocorrelated?
bf_data$resid<-resid
plot(bf_data$longitude,bf_data$latitude,cex=bf_data$resid/10)

#test for autocorrelation on residuals - just some basic diagnostics.... 
#get a matrix of distances between points
dists<-as.matrix(dist(cbind(bf_data$longitude, bf_data$latitude)))

# Take the inverse of the matrix values so that closer values have a larger weight and vs vs
dists.inv <- 1/dists
diag(dists.inv) <- 0 # replace the diagonal values with zero

# Computes Moran's I autocorrelation coefficient of x given a matrix of weights (here based on distance) 
# this is not recommended for determining whether or not a 'spatial' model is required
# but can help one view changes in residuals as covariates are added
Moran.I(resid, dists.inv)   # what do you conclude?

#now add more covariates and test 
m.1<-glm(outcome1~alt+temp+rain+ndvi,data=bf_data,family="binomial") # what are the model results?

print(summary(m.1))
confint(m.1)
exp(coef(m.1))
exp(confint(m.1))

resid2<-residuals(m.1)

#how do they look now?
bf_data$resid2<-resid2
plot(bf_data$longitude,bf_data$latitude,cex=bf_data$resid2/10)

# Computes Moran's I again
Moran.I(resid2, dists.inv)  # what do you think?


# we can also create neighborhood matrices - k nearest or based on distance
# and work with these data the way we did areal data
coords<-coordinates(bf_data.xy)

#setting k nearest neighbor matrix 
IDs<-seq(1:length(bf_data[,1]))
dist_nbdist <- knn2nb(knearneigh(coords, k=10),row.names=IDs)   # 5 nearest neighbor points 

#plotting points and neighbors - try this for a few different values of k
par(mfrow=c(1,1))
plot(bf_data.xy)
plot(dist_nbdist, coords, col="green",add=T)

#using this neighbor matrix we can use the function moran.test or moran.mc
BF_w<-nb2listw(dist_nbdist)
moran.test(resid2,listw=BF_w)   # again - use with caution for non-linear models!  

#we can also calculate distance based weights matrices
dsts<-unlist(nbdists(dist_nbdist, 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
BF_nb0.5 <- dnearneigh(coords, d1=0, d2=0.5*max_1nn, row.names=IDs)
BF_nb1 <- dnearneigh(coords, d1=0, d2=1*max_1nn, row.names=IDs)

#plot these - how does this compare to the k method?
par(mfrow=c(1,1))
plot(bf_data.xy)
plot(BF_nb0.5,coords,col="green",add=T)

BF_wd<-nb2listw(BF_nb0.5)
moran.test(resid2,list=BF_wd)

#using one of these neighborhood matrices, we can run a CAR model on the point data, using a binomial distribution
#weights matrix for CAR model
W.mat<-nb2mat(BF_nb0.5,style="B")

form<-positives~alt

model.spatial<-S.CARbym(formula=form,family="binomial",data=bf_data,trials=bf_data$examined,W=W.mat,burnin=2000,n.sample=10000,thin=10)

model.nonspatial<-S.glm(formula=form,family="binomial",data=bf_data,trials=bf_data$examined,burnin=2000,n.sample=10000,thin=10)



#########
#now let's try a GAM model
m.3<-gam(outcome1~s(alt)+s(temp)+s(rain)+s(ndvi),data=bf_data,family="binomial")
summary(m.3)
plot(m.3,pages=1,residuals=T)

m.3b<-gam(outcome1~te(alt)+te(temp)+te(rain)+te(ndvi),data=bf_data,family="binomial")
summary(m.3b)
plot(m.3b,pages=1,residuals=T)

#check the residuals and k specification
gam.check(m.3b)

m.3c<-gam(outcome1~te(alt,bs="tp",k=8)+te(temp,bs="tp",k=8)+te(rain,bs="tp",k=8)+te(ndvi,bs="tp",k=8),data=bf_data,family="binomial")
summary(m.3c)
plot(m.3c,pages=1,residuals=T)

#how does the gam model compare to the glm model?
AIC(m.1,m.3)

#show does it compare to the glm model?


#let's add a spatial effect by including a 2d smoothing parameter on longitude and latitude 
m.4<-gam(outcome1~s(alt)+s(temp)+s(rain)+s(ndvi)+s(longitude,latitude),data=bf_data,family="binomial")

summary(m.4)
plot(m.4,pages=1)
vis.gam(m.4,view=c("longitude","latitude"),plot.type='contour')

#which model is better?
AIC(m.3,m.4)

dev.off()


##predicting from glm and gam models 

#fitted vs response
bf_data$pred<-predict.glm(m.1,type="response")
plot(bf_data$prev,bf_data$pred)
abline(0,1,col="red")

#new prediction sites
#load prediction data - 25245 locations over 3x3km grid
bf_predict<-read.table("prediction_bf2_binomial.txt",sep=",",header=T)
predict<-cbind(alt=bf_predict$alt.pred,temp=bf_predict$temp.pred,rain=bf_predict$rain.pred,ndvi=bf_predict$ndvi.pred)
predict<-data.frame(predict)

#predict to glm model
pred<-predict.glm(m.1,newdata=predict,type="response")
pred<-data.frame(pred)

bf_xyz<-cbind(bf_predict$x.pred,bf_predict$y.pred,pred)
bf_xyz<-data.frame(bf_xyz)
names(bf_xyz)<-c("longitude","latitude","pred")

bf_raster1<-rasterFromXYZ(bf_xyz, digits=2)

plot(bf_raster1)
points(bf_data.xy,cex=bf_data$prev)

#prediction based on GAM model

#fitted vs response
bf_data$pred.gam<-predict(m.4,type="response")
plot(bf_data$prev,bf_data$pred.gam)
abline(0,1,col="red")

#new prediction sites
bf_predict<-read.table("prediction_bf2_binomial.txt",sep=",",header=T)
predict<-cbind(alt=bf_predict$alt.pred,temp=bf_predict$temp.pred,rain=bf_predict$rain.pred,ndvi=bf_predict$ndvi.pred,
               longitude=bf_predict$x.pred,latitude=bf_predict$y.pred)
predict<-data.frame(predict)

pred<-predict(m.4,newdata=predict,type="response")
pred<-data.frame(pred)

bf_xyz<-cbind(bf_predict$x.pred,bf_predict$y.pred,pred)
bf_xyz<-data.frame(bf_xyz)
names(bf_xyz)<-c("longitude","latitude","pred")

bf_raster2<-rasterFromXYZ(bf_xyz, digits=2)

plot(bf_raster2)
points(bf_data.xy,cex=bf_data$prev)


#sample for validation
nsample<-round(nrow(bf_data)*.15)
bf_data.val<-bf_data

#data we will compare to
set.seed(15)
bf_data.pred<-bf_data.val[sample(nrow(bf_data.val),nsample),]

#data we will run model on
set.seed(15)
bf_data.val$positives[sample(nrow(bf_data.val),nsample)]<-NA
bf_data.val<-na.omit(bf_data.val)
bf_data.val<-data.frame(bf_data.val)

#run glm model again on test sample
outcome1.v<-cbind(bf_data.val$positives,bf_data.val$examined-bf_data.val$positives)
m.1v<-glm(outcome1.v~alt+temp+rain+ndvi,data=bf_data.val,family="binomial") # what are the model results?

glm.testpred<-predict.glm(m.1v,newdata=bf_data.pred,type="response")
plot(glm.testpred,bf_data.pred$prev)

library(Metrics)
rmse.glm<-rmse(bf_data.pred$prev,glm.testpred)

#run gam model again on test sample 
m.4v<-gam(outcome1.v~s(alt)+s(temp)+s(rain)+s(ndvi)+s(longitude,latitude),data=bf_data.val,family="binomial")

gam.testpred<-predict.gam(m.4v,newdata=bf_data.pred,type="response")
plot(gam.testpred,bf_data.pred$prev)

rmse.gam<-rmse(bf_data.pred$prev,gam.testpred)
#which model has gives better predictions?


###########
# run the model in INLA model
library(INLA)
library(lattice)

#SPDE models

bound2 <- inla.nonconvex.hull(as.matrix(coordinates(bf_data.xy)), convex=0.5, concave=-0.15)
mesh=inla.mesh.create.helper(as.matrix(coordinates(bf_data.xy)),boundary=bound2,offset=c(.1,.2),min.angle=c(15,15),max.edge=c(1,1))
plot(mesh)
points(bf_data.xy,pch=19,col=2)


##----SPDE model

spde<-inla.spde2.matern(mesh,alpha=2)

##--- Observation structure for estimation data

coords<-as.matrix(coordinates(bf_data.xy)) 

A.est =inla.spde.make.A(mesh,loc=coords)

st.b<-inla.stack(data=list(y=bf_data$positives,n=bf_data$examined),A=list(A.est,1),
                 effects=list(i=1:spde$n.spde,data.frame(b0=1,bf_data)),tag='est')

#prediction data set-up
coords.pred<-as.matrix(cbind(bf_predict$x.pred,bf_predict$y.pred))

A.pred =inla.spde.make.A(mesh,loc=coords.pred)

st.pred<-inla.stack(data=list(y=NA,n=NA),A=list(A.pred,1),
                    effects=list(i=1:spde$n.spde,data.frame(b0=1,predict)),tag='pred')


#test model
res.s<-inla(y ~ 0 + (b0) + alt + temp + rain + ndvi + f(i,model=spde),family="binomial",
            control.predictor=list(A=inla.stack.A(st.b)),Ntrials=examined,control.compute=list(dic=T),control.fixed=list(expand.factor.strategy="inla"),
            data=inla.stack.data(st.b))


#summary of fixed effects + dic
exp(res.s$summary.fix)
res.s$dic$dic

#compare non-spatial
res.ns<-inla(y ~ 0 + (b0) + alt + temp + rain + ndvi,family="binomial",
            control.predictor=list(A=inla.stack.A(st.b)),Ntrials=examined,control.compute=list(dic=T),control.fixed=list(expand.factor.strategy="inla"),
            data=inla.stack.data(st.b))

exp(res.ns$summary.fix)
res.ns$dic$dic


#full prediction model
stk.full<-inla.stack(st.b,st.pred)

res.p<-inla(y ~ 0 + b0 + alt + temp + rain + ndvi + f(i,model=spde),family="binomial",
            control.predictor=list(A=inla.stack.A(stk.full),compute=TRUE),control.fixed=list(expand.factor.strategy="inla"),
            data=inla.stack.data(stk.full),Ntrials=n,control.results=list(return.marginals.random=T,
                                                                          return.marginals.predictor=T),control.family=list(link="logit"))

#results
rf<-inla.spde.result(inla=res.p,name='i',spde=spde)
#rf$marginals.var[[1]]

#plot the spatial random effect
gproj<-inla.mesh.projector(mesh)
g.mean<-inla.mesh.project(gproj,res.p$summary.random$i$mean)

levelplot(g.mean,col.regions=topo.colors(99),mean='latent field mean',xlab='',ylab='',
          scales=list(draw=FALSE))

#plot the linear predictor 
igr<-inla.stack.index(stk.full,'pred')$data

mean<-res.p$summary.fitted.values[igr,"mean"]

invmean<-exp(mean)/(1+exp(mean))

rmean.build<-cbind(coords.pred,invmean)

raster.mean<-as.matrix(rmean.build)

rastermean.out<-rasterFromXYZ(raster.mean,digits=1)

plot(rastermean.out)
points(bf_data.xy,cex=bf_data$prev)

writeRaster(rastermean.out,file="meanpf.tif",format="GTiff",overwrite=T)



#1. Assignments  - test the rmse for the INLA model. HOw does it compare to the gam model?

#2. Add a covariate to the model - does it help predictions?
# check previous lectures for potential covariates...

#3. Test a glm and gam model for the magude data set - adding at least 2 potential covariates
# ones you know how to access:  altitude, distance to rivers,  any other? 
# which model is best and why (which predictors and glm or gam)
# extra points for mapping the predictions!