Pointless spatial modeling
Adam Howes
Abstract
Task We explore the code from the paper “Pointless spatial modeling” by Katie Wilson and Jon Wakefield.
Code
set.seed(10)
load("kenya.rdata")
#' Read in population data
#' Follow instructions in pointless-spatial-modeling README as to how to get this
pop.dat <- raster("gpw_v4_population_count_rev11_2015_2pt5_min.tif")
pop.dat <- setMinMax(pop.dat)
kenya.extent <- extent(33.9, 42, -4.7, 5.1)
pop.dat.kenya <- crop(pop.dat, kenya.extent)
kenya.df1 <- fortify(adm0) %>%
mutate(id = as.numeric(id))## Regions defined for each Polygonskenya.df47 <- fortify(adm1) %>%
mutate(id = as.numeric(id))## Regions defined for each Polygonskenya.df8 <- fortify(adm0.5) %>%
mutate(id = as.numeric(id) + 1)## Regions defined for each Polygons#' Set-up variables for simulation
#' * Create mesh for first 5 scenarios: "mesh.true"
#' * Create D matrix used for scenario 2 and 3: "D" (47 x m)
#' * Create D matrix used for scenario 4 and 5: "D8" (see note above)
#' * Add columns to kenya.data indicating which polygons the EAs reside in
# Create internal mesh points
points.kenya <- spsample(adm0, 2500, "regular")## Warning in proj4string(obj): CRS object has comment, which is lost in output; in tests, see
## https://cran.r-project.org/web/packages/sp/vignettes/CRS_warnings.htmlboundary.kenya <- inla.nonconvex.hull(
points.kenya@coords,
convex = 0.1,
resolution = c(120, 120)
)
#' Determine population near these mesh points
pop.at.meshX <- raster::extract(pop.dat.kenya, points.kenya@coords)
pop.at.meshX[is.na(pop.at.meshX)] <- 0
pop.at.mesh <- data.frame(points.kenya@coords)
pop.at.mesh$pop <- pop.at.meshX
names(pop.at.mesh)[1:2] <- c("mlong", "mlat")
#' Make the mesh
mesh.true <- inla.mesh.2d(
pop.at.mesh[, 1:2],
offset = c(1, 2),
max.edge = 1 * c(1, 5)
)
plot(mesh.true)
#' adm1
#' Determine which polys the mesh points live in
meshgrid <- data.frame(
Longitude = mesh.true$loc[,1],
Latitude = mesh.true$loc[,2]
)
coordinates(meshgrid) <- ~ Longitude + Latitude
proj4string(meshgrid) <- proj4string(adm1)
meshlocs <- over(meshgrid, adm1)[, 1]
tapply(rep(1, length(meshlocs)), meshlocs, sum)## 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
## 48 12 14 7 14 11 186 21 112 97 12 12 11 54 7 6 11 132 35 37 26 27 33 108 324 32 13 2 10 4 32 10 74 3 14 14 91 16 69 167 11 9
## 43 44 45 46 47
## 271 14 2 239 41mesh.df <- data.frame(
mlong = mesh.true$loc[, 1],
mlat = mesh.true$loc[, 2],
loc = meshlocs,
id = 1:length(meshlocs)
)
#' adm0.5
# Determine which polys the mesh points live in
#' There is slight misalignment of the boundaries, which should completely
#' overlap. Thus, we manually enter the provinces based on which of the 47
#' districts the mesh point is in.
meshlocs8 <- over(meshgrid, adm0.5)[, 1]
mesh.df$loc8X <- meshlocs8
mapping47to8 <- as.vector(apply(table(mesh.df[, c("loc", "loc8X")]), 1, which.max))
mesh.df$loc8 <- mapping47to8[mesh.df$loc]
#' Population at mesh including polygon location
pop.at.mesh <- dplyr::right_join(mesh.df, pop.at.mesh)## Joining, by = c("mlong", "mlat")#' Create D matrix
#' adm1
D <- inla.spde.make.A(
mesh = mesh.true,
loc = as.matrix(pop.at.mesh[,1:2]),
block = pop.at.mesh$loc,
weights = pop.at.mesh$pop
)
D.tmp <- list()
D.tmp$D <- D / rowSums(D)
D.tmp$mesh.weights <- colSums(D)
mesh.df$weight.unscaled <- D.tmp$mesh.weights
D <- D.tmp$D
mesh.df$weight.scaled <- colSums(D)
nmesh.area <- 1 / tapply(rep(1, nrow(mesh.df)), mesh.df$loc, sum)
mesh.df$weight.comp <- nmesh.area[mesh.df$loc] #' Comparison weight
#' adm0.5
D8 <- inla.spde.make.A(
mesh = mesh.true,
loc = as.matrix(pop.at.mesh[,1:2]),
block = pop.at.mesh$loc8,
weights = pop.at.mesh$pop
)
D.tmp8 <- list()
D.tmp8$D <- D8 / rowSums(D8)
D.tmp8$mesh.weights <- colSums(D8)
mesh.df$weight.unscaled8 <- D.tmp8$mesh.weights
D8 <- D.tmp8$D
mesh.df$weight.scaled8 <- colSums(D8)
nmesh.area8 <- 1 / tapply(rep(1, nrow(mesh.df)), mesh.df$loc8, sum)
mesh.df$weight.comp8 <- nmesh.area8[mesh.df$loc8] #' Comparison weight
#' Set-up parameters
alpha <- 0
sigma2 <- 0.25
lambda2 <- 1 - sigma2
erange <- sqrt(8) / exp(0.5)
theta2 <- log(sqrt(8) / erange)
kappa <- exp(theta2)
tau2 <- 1 / (4 * pi * kappa^2 * lambda2)
theta1 <- log(sqrt(tau2))
#' Set-up field
spde <- inla.spde2.matern(mesh.true)
Q.true <- inla.spde.precision(spde, c(theta1, theta2))
field.true <- as.numeric(inla.qsample(1, Q = Q.true, seed = 50L))## Warning in inla.qsample(1, Q = Q.true, seed = 50L): Since 'seed!=0', parallel model is disabled and serial model is selected#' Connecting field to data
proj.survey <- inla.mesh.projector(
mesh.true,
loc = as.matrix(kenya.data[, c("LONGNUM", "LATNUM")])
)
field.survey <- inla.mesh.project(proj.survey, field = field.true)
kenya.data$field <- field.survey
kenya.data$mu <- alpha + kenya.data$field
#' Make matrix that will be used to obtain predictions of the surface
pred.info <- MakePredictionAMatrix(mesh.true, adm0, dims=c(110, 110))
A.pred <- pred.info$A.pred
stkgrid <- inla.stack(
data = list(y = NA, N = NA),
A = list(A.pred, 1),
effects = list(i = 1:spde$n.spde, alpha = rep(1, nrow(A.pred))), tag = 'prd.gr'
)
rel.points.clip <- pred.info$rel.points.clip
#' Simulate data
#' Census
#' Data simulated on 1km x 1km grid
pop.data.kenya.points <- data.frame(rasterToPoints(pop.dat.kenya, spatial = F))
names(pop.data.kenya.points) <- c("long", "lat", "pop")
coordinates(pop.data.kenya.points) <- ~ long + lat
pop.data.kenya.points@proj4string <- adm0@proj4string
pop.pt.insideX <- pop.data.kenya.points[adm0, ]
pop.pt.inside <- as.data.frame(pop.pt.insideX)
coordinates(pop.pt.inside) <- ~ long + lat
pop.pt.inside@proj4string <- adm0@proj4stringfunctions.R
# Make a smaller version of the A matrix so that predicting the risk surface
# is less computationally involved
#
# Args:
# mesh.true: The mesh
# adm0: Spatial polygons dataframe that is for the boundary of Kenya
#
# Return: list with prediction matrix
# and points where predictions will be made
MakePredictionAMatrix <- function(mesh.true, adm0, dims = c(120, 120)) {
smallproj <- inla.mesh.projector(mesh.true, dims=dims)
x <- smallproj$x
y <- smallproj$y
all.temp <- (expand.grid(x, y , KEEP.OUT.ATTRS=F))
names(all.temp) <- c("long", "lat")
coordinates(all.temp) <- ~ long + lat
all.temp@proj4string <- adm0@proj4string
rel.points <- all.temp[adm0, ]
rel.points.clip <- as.data.frame(rel.points)
A.pred <- inla.spde.make.A(mesh.true, loc=as.matrix(rel.points.clip))
return(list(A.pred=A.pred, rel.points.clip=rel.points.clip))
}Original computing environment
sessionInfo()## R version 4.2.0 (2022-04-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Monterey 12.6
##
## Matrix products: default
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] dplyr_1.0.9 INLA_22.05.07 foreach_1.5.2 Matrix_1.5-1 knitr_1.39 rgdal_1.5-30 raster_3.5-15
## [8] sp_1.4-7 rstan_2.21.5 ggplot2_3.3.6 StanHeaders_2.21.0-7 sf_1.0-7 bsae_0.2.7
##
## loaded via a namespace (and not attached):
## [1] viridis_0.6.2 sass_0.4.1 jsonlite_1.8.0 viridisLite_0.4.0 ids_1.0.1 splines_4.2.0 bslib_0.3.1 RcppParallel_5.1.5
## [9] assertthat_0.2.1 askpass_1.1 highr_0.9 stats4_4.2.0 yaml_2.3.5 pillar_1.7.0 lattice_0.20-45 glue_1.6.2
## [17] uuid_1.1-0 digest_0.6.29 colorspace_2.0-3 htmltools_0.5.2 plyr_1.8.7 pkgconfig_2.0.3 purrr_0.3.4 scales_1.2.0
## [25] processx_3.5.3 terra_1.5-21 tibble_3.1.7 openssl_2.0.1 proxy_0.4-26 bayesplot_1.9.0 generics_0.1.2 ellipsis_0.3.2
## [33] withr_2.5.0 cli_3.3.0 magrittr_2.0.3 crayon_1.5.1 evaluate_0.15 ps_1.7.0 fansi_1.0.3 class_7.3-20
## [41] pkgbuild_1.3.1 tools_4.2.0 loo_2.5.1 prettyunits_1.1.1 lifecycle_1.0.1 matrixStats_0.62.0 stringr_1.4.0 munsell_0.5.0
## [49] callr_3.7.0 compiler_4.2.0 jquerylib_0.1.4 e1071_1.7-9 rlang_1.0.2 classInt_0.4-3 units_0.8-0 grid_4.2.0
## [57] ggridges_0.5.3 iterators_1.0.14 rstudioapi_0.13 rmarkdown_2.14 orderly_1.4.3 gtable_0.3.0 codetools_0.2-18 inline_0.3.19
## [65] DBI_1.1.2 R6_2.5.1 splancs_2.01-43 gridExtra_2.3 fastmap_1.1.0 utf8_1.2.2 KernSmooth_2.23-20 stringi_1.7.6
## [73] Rcpp_1.0.8.3 vctrs_0.4.1 tidyselect_1.1.2 xfun_0.31