# Cheat sheet #1.
nairobi_cases$case_or_control <- ifelse(nairobi_cases$case=="Yes","case", "control")

Nairobi_Owin <- owin(xrange=range(nairobi_cases$lng),yrange=range(nairobi_cases$lat))

nairobi_cases_ppp <- ppp(nairobi_cases$lng, nairobi_cases$lat, 
                         window = Nairobi_Owin, 
                         marks=as.factor(nairobi_cases$case_or_control))

risk_est <- risk(nairobi_cases_ppp,tolerate = TRUE) # adapt=TRUE uses adaptive smoothing 
# (i.e. different smount of smoothing dependent on where you are)
plot(risk_est)

risk_raster <- raster(risk_est$rr, crs = crs(NAM_Adm0))


# Then plot using the leaflet package
pal = colorNumeric(palette=tim.colors(64), domain=risk_raster[], na.color = NA)
leaflet() %>% addTiles("http://{s}.basemaps.cartocdn.com/light_all/{z}/{x}/{y}.png") %>% 
  addRasterImage(risk_raster, opacity=0.6, col = pal)


## 2.
CV_idw_1<-idw(BF_malaria_data_ppp, power=1, at="points")
plot(BF_malaria_data_ppp$marks, CV_idw_1)

# Calc MSE
library(Metrics)
mse(BF_malaria_data_ppp$marks,CV_idw_1)


best_IDW_power <- function(ppp_data, powers){
  
  results <- NULL
  
  for(i in 1:length(powers)){
    CV_idw <- idw(ppp_data, power=powers[i], at="points")
    mse_CV_idw <- mse(ppp_data$marks, CV_idw)
    
    # add the result to your result string
    results <- c(results, mse_CV_idw)
  }
  
  best_power <- powers[which.min(results)]
  
  plot(powers, results, ylab="MSE")
  
  # Package output into a list
  return(list (power = best_power,
               MSE = results[which.min(results)],
               CV_predictions = idw(ppp_data, best_power, at="points")))
}

best_power <- best_IDW_power(BF_malaria_data_ppp, powers = seq(0.1, 3, 0.1))
best_power$power
best_power$MSE
best_power$CV_predictions

plot(idw(BF_malaria_data_ppp, power=best_power$power, at="pixels"))


## 3.
# 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/100,window=HK_Window)

best_power_HK <- best_IDW_power(HK_ppp, powers = seq(0.1, 3, 0.1))

# Generate and plot
IDW_layer <- idw(HK_ppp, power=best_power_HK$power, at="pixels")
image(IDW_layer,col=heat.colors(20))

# 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)
Vario$n # Shows you the number in each bin
min(Vario$n)# should be at least 30 pairs in each bin
plot(Vario)

# 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_lin
VarioMod_exp
VarioMod_sph # spherical model has lower sum of squares so 'better'

# Use IDW grid to krig to
pred_grid_x<-rep(IDW_layer$xcol,length(IDW_layer$yrow))
pred_grid_y<-sort(rep(IDW_layer$yrow,length(IDW_layer$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))

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

# Back transform to prevalence
KrigPred_prev<-inv.logit(KrigPred$predict)
KrigPred_prev_raster <- rasterFromXYZ(data.frame(x = pred_grid_x,
                                                 y = pred_grid_y,
                                                 z = KrigPred_prev),
                                      crs = crs(NAM_Adm0))

plot(KrigPred_prev_raster)

# Or using leaflet
colPal <- colorNumeric(tim.colors(), KrigPred_prev_raster[], na.color = NA)
leaflet() %>% addTiles() %>% addRasterImage(KrigPred_prev_raster,
                                            col = colPal,
                                            opacity = 0.7) %>%
  addLegend(pal = colPal, values = KrigPred_prev_raster[])


# Cross validate
xvalid_result_HK <- xvalid(HK.geo, model = VarioMod_sph)
plot(xvalid_result_HK$data, xvalid_result_HK$predicted)
abline(0,1, col="red")

# Compare to IDW
par(mfrow=c(2,2))
plot(IDW_layer, col = tim.colors(64))
plot(KrigPred_prev_raster, col = tim.colors(64))
plot(HK_ppp$marks, best_power_HK$CV_predictions, asp=1, 
     xlab="Observed", ylab="Predicted")
abline(0,1,col="red")
plot(inv.logit(xvalid_result_HK$data), inv.logit(xvalid_result_HK$predicted), asp=1,
     xlab="Observed", ylab="Predicted")
abline(0,1,col="red")

# Compare MSE
best_power_HK$MSE
mse(inv.logit(xvalid_result_HK$data), inv.logit(xvalid_result_HK$predicted))