library(spatstat)
library(raster)
library(sp)
library(lgcp)
library(geoR)
library(gtools)
library(plotGoogleMaps)
library(lme4)

setwd("C:\\Users\\sturrockh\\Documents\\Work\\teaching\\Spatial Epi UCSF")
# Open Namibia malaria case data
CaseControl<-read.csv("CaseControl.csv")
Cases<-CaseControl[CaseControl$case==1,]
Controls<-CaseControl[CaseControl$case==0,]

# get NAmibia boundary data
NAM_Adm0<-getData('GADM',country='NAM',level=0)
Nam_Owin <- as(as(NAM_Adm0, "SpatialPolygons"), "owin")
Nam_Owin<-owin(xrange=range(CaseControl$long),yrange=range(CaseControl$lat))

# Generate point pattern object of points
points<-ppp(runif(50,0,1),runif(50,0,1))
plot(points)

# Plot kernel density estimate
plot(density(points))

# Try with different bandwidths
plot(density(points,0.02))
plot(density(points,0.08))
plot(density(points,bw.diggle)) # automatic bandwidth selection based on cross-validation


# Get Uganda outlines
UGA_Adm0<-getData('GADM',country='UGA',level=0)
owin_obj <- as(as(UGA_Adm0, "SpatialPolygons"), "owin")

# Read in hookworm data
HK<-read.table("C:\\Users\\sturrockh\\Documents\\Work\\teaching\\3135\\Kriging prac\\tanzania_uganda_hkprev.txt",
               header=T)

# Visualize
coordinates(HK)<-~x+y
proj4string(HK)<- CRS('+init=epsg:4326') # Lat/long WGS84 

m<-bubbleGoogleMaps(HK,
                    filename='HK.htm',
                    mapTypeId='TERRAIN',
                    zcol='Hookworm_prev',
                    layerName = 'Hookworm prevalence',
                    max.radius=8000) # in meters

# Inverse distance weighting
HK_Window<-owin(xrange=range(HK$x),yrange=range(HK$y))
HK_ppp<-ppp(HK$x,HK$y,marks=HK$Hookworm_prev,window=HK_Window)

par(mfrow=c(2,2))
plot(idw(HK_ppp, power=0.2, at="pixels"),col=heat.colors(20))
plot(idw(HK_ppp, power=0.5, at="pixels"),col=heat.colors(20))
plot(idw(HK_ppp, power=1, at="pixels"),col=heat.colors(20))
plot(idw(HK_ppp, power=2, at="pixels"),col=heat.colors(20)) # Larger power puts more weight on nearer values


# Kriging
# Before Kriging, good to transform prevalence data to something vaguely normal
# Here we use the logistic transformation (log odds)
HK$LogOdds<-logit((HK$Hookworm_prev/100)+0.001) # +0.001 for values of 0

# First have to create a geodata object with the package GeoR
# this wants dataframe of x,y and data 
HK.geo<-as.geodata(cbind(HK$x,HK$y,HK$LogOdds))

# We can plot a summary plot
plot(HK.geo, lowes=T) # the lowes option gives us lowes curves for relationship with x and y
# Seems like there is a non-linear trend with x

plot(HK.geo, lowes=T,trend="2nd") # Trend option regresses on x and y

# Now generate and plot a variogram
MaxDist=max(dist(cbind(HK$x,HK$y)))/2 # the max distance you should estimate is half max interpoint distance 
VarioCloud<-variog(HK.geo,max.dist=MaxDist,option="cloud",
                   trend="2nd")
plot(VarioCloud) # all pairwise comparisons

# To produce binned variogram
Vario<-variog(HK.geo,max.dist=MaxDist)
plot(Vario)
Vario<-variog(HK.geo,max.dist=MaxDist,trend="2nd",uvec=seq(0.01,MaxDist,0.1))
Vario$n # Shows you the number in each bin
min(Vario$n)# should be at least 30 pairs in each bin
plot(Vario,pch=16)

# Fit variogram model by minimized least sqaures
VarioMod_lin<-variofit(Vario, cov.model = "linear")
VarioMod_sph<-variofit(Vario, cov.model = "sph")
VarioMod_exp<-variofit(Vario, cov.model = "exp")

# plot results
lines(VarioMod_lin,col="green",lwd=2)
lines(VarioMod_sph,col="blue",lwd=2)
lines(VarioMod_exp,col="red",lwd=2)


VarioMod_exp
VarioMod_sph # spherical model has lower sum of squares so 'better'

# Use variogram to Krig values at prediction locations 
# First get grid of points from the IDW example for comparison
# could use the expand.grid function 
IDW<-idw(HK_ppp, power=0.2, at="pixels")
pred_grid_x<-rep(IDW$xcol,length(IDW$yrow))
pred_grid_y<-sort(rep(IDW$yrow,length(IDW$xcol)))
pred_grid<-cbind(pred_grid_x,pred_grid_y)


KrigPred <- krige.conv(HK.geo, loc=pred_grid,
                       krige=krige.control(obj.model=VarioMod_sph,
                                           trend.d="2nd",trend.l="2nd"))

# Visualize predictions
image(KrigPred,col=heat.colors(50))

# Back transform to prevalence
KrigPred_prev<-inv.logit(KrigPred$predict)
KrigPred_SPDF<-SpatialPixelsDataFrame(pred_grid,
                data.frame(Predictions=KrigPred_prev),
                proj4string = CRS('+init=epsg:4326'))

# If you want to visualize over Google Maps
plotGoogleMaps(KrigPred_SPDF, mapTypeId='ROADMAP',col=heat.colors(20))


# Random effects
# If we look at the trend with lat and long
plot(HK.geo,lowes=T) # Seems like there might be a trend with longitude

# Load individual level data we can use to run a logistic regression, regressing on 
# longitude (with cluster level random effect)
setwd("C:\\Users\\sturrockh\\Documents\\Work\\teaching\\Spatial Epi UCSF\\")
load("HK_Individual.RData")

# Use lme4 package to run a random effects model
Mod1<-glmer(HK_Pos ~ x + (1 | cluster_ID),family="binomial",data=HK_Individual)
Mod2<-glmer(HK_Pos ~ y + (1 | cluster_ID),family="binomial",data=HK_Individual)
hist(unlist(ranef(Mod1))) # fairly normal..

# Plot and fit variogram
HK_RE.geo<-as.geodata(cbind(HK$x,HK$y,unlist(ranef(Mod1))))
Vario_RE<-variog(HK_RE.geo,max.dist=MaxDist)
plot(Vario_RE)
VarioMod_RE_sph<-variofit(Vario_RE, cov.model = "sph")
VarioMod_RE_exp<-variofit(Vario_RE, cov.model = "exp")
lines(VarioMod_RE_sph,col="blue",lwd=2)
lines(VarioMod_RE_exp,col="red",lwd=2)

# Now make predictions using longitude to predict large scale trend and kriging to predict random effect
LongPred<-predict(Mod1,re.form=NA,newdata=data.frame(x=pred_grid_x))
LongPred_SPDF<-SpatialPixelsDataFrame(pred_grid,
                                      data.frame(LongPred=LongPred))
image(LongPred_SPDF)

# Now Krig the random effect
KrigPred_RE <- krige.conv(HK_RE.geo, loc=pred_grid,
                       krige=krige.control(obj.model=VarioMod_RE_sph))
image(KrigPred_RE)

# Add this to your longitudinal trend
FinalPred_LogOdds<-LongPred+KrigPred_RE$predict

# Back transform to prevalence
FinalPred_prevalence<-inv.logit(FinalPred_LogOdds)

# Visualize
FinalPred_prevalence_SPDF<-SpatialPixelsDataFrame(pred_grid,
                                      data.frame(FinalPred_prevalence=FinalPred_prevalence),
                                      proj4string = CRS('+init=epsg:4326'))
plotGoogleMaps(FinalPred_prevalence_SPDF, mapTypeId='ROADMAP',col=heat.colors(20))