# install and load relevant packages. 
pacman::p_load(Metrics, spatstat, raster, sp, 
  geoR, gtools, lme4, leaflet, oro.nifti, tidyverse,
  rangeBuilder, maptools)
# load data from local source.
nairobi_cc <- read_csv("Assignment/cases_nairobi.csv")

# format cases and controls as binary integers.
nairobi_cc <- 
  nairobi_cc %>% 
  mutate(case = case_when(
    case == "Yes" ~ 1,
    case == "No" ~ 0,))

# split datasets on case vs. control outcome.
nairobi_cases <- nairobi_cc  %>% 
  filter(case == 1)

nairobi_controls <- nairobi_cc %>% 
  filter(case == 0)

# import GADM file for crs reference. 
ken_adm0 <- raster::getData('GADM', 
                            country = "KEN", 
                             level=0)
# format as spatial data.
case_control_SPDF <- SpatialPointsDataFrame(
  coords = nairobi_cc[,c("lng", "lat")], 
  data = nairobi_cc[,"case"])

# palette for case vs. color outcome
case_color_scheme <- colorFactor(c("blue", "red"),
                                  case_control_SPDF$case)

# Map outcome. 
leaflet() %>% addTiles() %>%  
  addCircleMarkers(data = case_control_SPDF, 
                   color = case_color_scheme(case_control_SPDF$case),
                   radius = .1) %>% 
  addLegend(pal = case_color_scheme, values = case_control_SPDF$case, title = "Cases & Controls <p> Simulated : Nairobi, Kenya", 
            position = "bottomleft")
# format owin object for density plot
nai_owin <- owin(xrange = range(nairobi_cc$lng), 
                                yrange=range(nairobi_cc$lat))

# format point pattern object with owin spatial range. 
cases_ppp <- ppp(nairobi_cases$lng, nairobi_cases$lat, 
                 window = nai_owin)

# plot point pattern object.
plot(cases_ppp)

# create raster layer from point pattern object.
case_density <- density(cases_ppp)
plot(case_density)

# set graphic output.
par(mfrow = c(2,2), mar= c(1,1,1,1))

# explore bandwidth parameters. 
plot(density(cases_ppp, .01), main = "BNDWT = .01")
plot(density(cases_ppp, .0075), main = "BNDWT = .0075")
plot(density(cases_ppp, .005), main = "BNDWT = .005")
plot(density(cases_ppp, bw.ppl), main = "BNDWT = Automatic")

# select bandwidth and format to raster layer
density_raster <- raster(density(cases_ppp, .005), crs = crs(ken_adm0))

# create palette.
dens_colPal <- colorNumeric(tim.colors(), 
                       density_raster[], na.color = NA)

# map density raster. 
leaflet() %>% 
  addProviderTiles("CartoDB.Positron") %>%  
  addRasterImage(density_raster, opacity = .5) 
# point pattern object with cases and controls.
case_control_ppp <- ppp(nairobi_cc$lng, nairobi_cc$lat,
                       window = nai_owin, 
                       marks = as.factor(nairobi_cc$case))


plot(case_control_ppp)

# optimize bandwidth parameter. 
optimal_sig <- bw.relrisk(case_control_ppp)

# estimate spatial variance of relative risk. 
rel_risk_est <- relrisk.ppp(case_control_ppp, sigma = optimal_sig, relative = T)
plot(rel_risk_est)

# create raster layer of estimation
risk_raster <- raster(rel_risk_est, crs = crs(ken_adm0))

# palette.
pal <- colorNumeric(palette = tim.colors(64), 
                    domain = values(risk_raster), na.color = NA)

# map raster layer. 
leaflet() %>%  addProviderTiles("CartoDB.Positron") %>%  
  addRasterImage(risk_raster, opacity = .6, col = pal) %>% 
    addLegend(pal = pal, 
              values = risk_raster[], 
              title = "Spatial Risk Ratio, <p> Simulated: 
              <p> Nairobi, Kenya", 
              position = "bottomleft")

IDW interpolation:

# import data from local source.
tza_uga_hk <- read_csv("Assignment/tanzania_uganda_hkprev.csv")

# clean and format prevalence as percent. 
tza_uga_hk <- 
  tza_uga_hk %>% 
  rename(long = x,
         lat = y) %>% 
  mutate(Hookworm_prev_perc = (Hookworm_prev_perc / 100))
# import GADM for crs reference. 
tza_adm0 <- raster::getData('GADM', 
                            country = "TZA", 
                            level=0)

# format data as spatial object. 
tza_uga_SPDF <- SpatialPointsDataFrame(
  coords = tza_uga_hk[,c("long", "lat")], 
  data = tza_uga_hk[,"Hookworm_prev_perc"])
# create bounding box owin object.
tza_uga_window <- owin(tza_uga_SPDF@bbox[1,], tza_uga_SPDF@bbox[2,])

# create point pattern object. 
tza_uga_hk_ppp <- ppp(tza_uga_hk$long, 
                            tza_uga_hk$lat, 
                            marks = tza_uga_hk$Hookworm_prev_perc,
                            window = tza_uga_window)


# set graphic output. 
par(mfrow = c(2,2), mar= c(1,1,1,1))

# explore power parameters.
plot(idw(tza_uga_hk_ppp, power=0.1, at="pixels"),
     col=heat.colors(20), main="power = 0.1") 
plot(idw(tza_uga_hk_ppp, power=0.5, at="pixels"),
     col=heat.colors(20), main="power = 0.5")
plot(idw(tza_uga_hk_ppp, power=1, at="pixels"),
     col=heat.colors(20), main="power = 1")
plot(idw(tza_uga_hk_ppp, power=2, at="pixels"),
     col=heat.colors(20), main="power = 2")

Optimzing power parameter:

# set range of test values. 
powers <- seq(.1, 2.5, .1) 

# create placeholder for mean squared error. 
mse_result <- NULL 

# function to evaluate mean squared error at sequential intervals defined in "powers" object.
for(power in powers){cv_idw <- 
  idw(tza_uga_hk_ppp, power=power, at="points") 
  mse_result <- c(mse_result, mse(tza_uga_hk_ppp$marks,cv_idw)) 
}

# create optimal minimum mse object. 
optimal_power <- powers[which.min(mse_result)] 

# plot output of function. 
plot(powers, mse_result) 

# apply inverse distance weights to point pattern data with optimal power value. 
cv_idw_opt <- idw(tza_uga_hk_ppp, power=optimal_power, at="points") 

# plot to explore differences between observed and predicted prevalence values. 
plot(tza_uga_hk_ppp$marks, cv_idw_opt, xlab="Observed prevalence",
     ylab="Predicted prevalence")

# create raster layer with optimized idw object.
tza_uga_hk_raster <- raster(idw(tza_uga_hk_ppp, 
                                power=optimal_power, 
                                at="pixels"), crs = crs(tza_adm0))
# palette.
colPal <- colorNumeric(tim.colors(), 
                       tza_uga_hk_raster[], na.color = NA)

# map. 
leaflet() %>% addProviderTiles("CartoDB.Positron")  %>%
   addCircleMarkers(data = tza_uga_SPDF,
    radius = ~ Hookworm_prev_perc, color = "black",
    stroke = FALSE, fillOpacity = 0.4 ) %>% 
  addRasterImage(tza_uga_hk_raster, col = colPal, opacity=0.5) %>%
  addLegend(pal = colPal, values = tza_uga_hk_raster[], 
            title = "Inverse Distance Weighted <p> Hookworm 
            Prevalence: <p> Tanzania & Uganda Region", position = "bottomleft")

Kriging:

# create "geodata" object. 
tza_uga_hk_geo <- as.geodata(tza_uga_hk[, 
                  c("long", "lat", "Hookworm_prev_perc")])

# set graphics
par(mfrow=c(4,6))

# evaluate data for first and second order trend.  
plot(tza_uga_hk_geo, lowes=T, trend = "1st")

plot(tza_uga_hk_geo, lowes=T, trend = "2nd")

# create variogram
max_dist <- max(dist(tza_uga_hk[,c("long","lat")])) / 2 

# create data cloud
vario_cloud<-variog(tza_uga_hk_geo, option="cloud", 
                    max.dist= max_dist)
## variog: computing omnidirectional variogram
# plot cloud
plot(vario_cloud)

# bin points by distance, set range parameters. 
vario <- variog(tza_uga_hk_geo,max.dist=max_dist, 
# specify 2nd order trend based on Lowes output.                 
                trend = "2nd",
        uvec=seq(0.01,max_dist,0.2)) 
## variog: computing omnidirectional variogram
# plot
plot(vario)

# check if bin vlaues are eacg greater than 30
min(vario$n) > 30
## [1] TRUE
# plot binned variogram. 
plot(vario,pch=16)

# fit spherical model 
variomod_sph <-variofit(vario, cov.model = "sph")
## variofit: covariance model used is spherical 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq kappa
## initial.value "0.04"  "0.83" "0"   "0.5"
## status        "est"   "est"  "est" "fix"
## loss value: 0.820712114225986
# fit exponetial model 
variomod_exp <-variofit(vario, cov.model = "exp")
## variofit: covariance model used is exponential 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq  kappa
## initial.value "0.04"  "0.83" "0.01" "0.5"
## status        "est"   "est"  "est"  "fix"
## loss value: 0.806885323394223
# plot results from each model
plot(vario,pch=16)
lines(variomod_sph,col="blue",lwd=2)
lines(variomod_exp,col="red",lwd=2) 

# summarise spherical model and extract sum.of.squares value. 
sph_summary <- summary(variomod_sph)
sph_summary$sum.of.squares
##     value 
## 0.2489733
# summarise exponential model and extract sum.of.squares value. 
exp_summary <- summary(variomod_exp)
exp_summary$sum.of.squares
##     value 
## 0.3079354
# create prediction grid. 
IDW <- idw(tza_uga_hk_ppp, power=optimal_power, 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)

# "krig" to the points with the spherical model (smaller mse between the two tested). 
krig_pred <- krige.conv(tza_uga_hk_geo, loc=pred_grid,
            krige=krige.control(obj.model=variomod_sph))
## krige.conv: model with mean given by a 2nd order polynomial on the coordinates
## krige.conv: Kriging performed using global neighbourhood
# create raster image of kriged output. 
image(krig_pred,col=heat.colors(50))

# create raster of kriged prediction. 
krig_pred_raster <- rasterFromXYZ(data.frame(x=pred_grid_x,
                                 y=pred_grid_y,
                                 z=krig_pred$predict), 
                                 crs = crs(tza_adm0))

# plot kriged prevalence prediction. 
plot(krig_pred_raster)
points(tza_uga_hk[,c("long","lat")],
       cex =  tza_uga_hk$Hookworm_prev_perc)

# palette for map
pred_colPal <- colorNumeric(tim.colors(), 
                            krig_pred_raster[], na.color = NA)

# map with legend. 
leaflet() %>% addProviderTiles("CartoDB.Positron")  %>% 
  addCircleMarkers(data = tza_uga_SPDF,
    radius = ~ Hookworm_prev_perc, color = "black",
    stroke = FALSE, fillOpacity = 0.3 ) %>% 
  addRasterImage(krig_pred_raster, col = pred_colPal, opacity=0.5) %>%
  addLegend(pal = pred_colPal, values = krig_pred_raster[], 
            title = "Kriged hookworm prevalence <p> Tanzania & Uganda Region", position = "bottomleft")

Evaluating IDW and Krig models with Mean Squared Error:

# cross validate geo object with `xvalid()` function. 
xvalid_result <- xvalid(tza_uga_hk_geo, model = variomod_sph)
## xvalid: number of data locations       = 393
## xvalid: number of validation locations = 393
## xvalid: performing cross-validation at location ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346, 347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, 359, 360, 361, 362, 363, 364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392, 393, 
## xvalid: end of cross-validation
# plot to evalute prediction.
plot(xvalid_result$data,xvalid_result$predicted, asp=1,
     xlab = "Observed prevalence", ylab="Cross-validated predicted prevalence")
abline(0,1)

# "nudge" values for complete positive in data. 
tza_uga_hk$Hookworm_prev_perc_adj <- 
  tza_uga_hk$Hookworm_prev_perc + .001

# compute logit function on adjusted points
tza_uga_hk$Hookworm_prev_perc_logit <- logit(tza_uga_hk$Hookworm_prev_perc_adj)

# reformat as geo object
tza_uga_hk_geo_logit <- as.geodata(tza_uga_hk[,
                c('long','lat',"Hookworm_prev_perc_logit")])

# create variogram from adjusted logit returns. 
vario_logit <- variog(tza_uga_hk_geo_logit, max.dist = max_dist)
## variog: computing omnidirectional variogram
# fit to spherical model
variomod_sph_logit <- variofit(vario_logit, cov.model = 'sph')
## variofit: covariance model used is spherical 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq  kappa
## initial.value "1.67"  "1.65" "0.22" "0.5"
## status        "est"   "est"  "est"  "fix"
## loss value: 2443.57768681399
# cross validate with spherical model. 
xvalid_result_logit <- xvalid(tza_uga_hk_geo_logit, model = variomod_sph_logit)
## xvalid: number of data locations       = 393
## xvalid: number of validation locations = 393
## xvalid: performing cross-validation at location ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346, 347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, 359, 360, 361, 362, 363, 364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392, 393, 
## xvalid: end of cross-validation
# inverse transform, "un-logitify" logit values.
xvalid_result_inv_logit <- inv.logit(xvalid_result_logit$predicted)
# print MSE returns.
paste("MSE with IDW:    ", mse(cv_idw_opt, tza_uga_hk$Hookworm_prev_perc))
## [1] "MSE with IDW:     0.0305233269189452"
paste("MSE with Kriging:", mse(xvalid_result_inv_logit, tza_uga_hk$Hookworm_prev_perc))
## [1] "MSE with Kriging: 0.0327923337495328"
krig_pred_logit <- krige.conv(tza_uga_hk_geo_logit, loc=pred_grid,
    krige=krige.control(obj.model=variomod_sph_logit))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
krig_pred_raster <-rasterFromXYZ(
  data.frame(x=pred_grid_x, y=pred_grid_y,
             z=inv.logit(krig_pred_logit$predict)))

Visualize Differentials in Model Outputs:

# format overlay raster
overlay_raster <- raster::overlay(tza_uga_hk_raster, 
                                  krig_pred_raster, 
                                  
# function to subtract second from first layer object. 
                                  fun = function(x,y){return(x-y)})

# palette
comp_colPal <- colorNumeric(heat.colors(50), 
                            overlay_raster[], na.color = NA)

# map with legend. 
leaflet() %>% addProviderTiles("CartoDB.Positron")  %>% 
  addRasterImage(overlay_raster, col = comp_colPal, opacity=0.8) %>%
   addCircleMarkers(data = tza_uga_SPDF,
    radius = ~ Hookworm_prev_perc, color = "black",
    stroke = FALSE, fillOpacity = 0.4 ) %>% 
  addLegend(pal = comp_colPal, values = overlay_raster[], 
            title = "Comparison Plot, <p> (IDW - Krig) differential: <p> Tanzanina & Uganda Region", position = "bottomleft")

~fin