Loaloa prevalence with FEM Matérn

library(adlaplace)
library(adlaplaceGrf)

This vignette fits a binomial geostatistical model for loaloa prevalence using a village-level IID effect and a FEM Matérn spatial field (shape = 2, ν = 2), matching the Matérn shape used in geostatsp examples. When INLA is available, the same model is also fit with geostatsp::glgm for comparison.

Data

data("loaloa", package = "geostatsp")

loaloa <- terra::unwrap(loaloa)
elev_r <- terra::unwrap(elevationLoa)
evi_r <- terra::unwrap(eviLoa)

elev_s_r <- (elev_r - 500) / 100
evi_s_r <- (
  evi_r - mean(terra::values(evi_r), na.rm = TRUE)
) / sd(terra::values(evi_r), na.rm = TRUE)

loaloa_for_fit <- geostatsp:::gm.dataSpatial(
  formula = y ~ elev + evi,
  data = loaloa,
  covariates = list(elev = elev_s_r, evi = evi_s_r),
  grid = geostatsp::squareRaster(loaloa, cells = 200, buffer = 1e4)
)

Knot lines

Use a coarse knot raster over a buffered region around the villages; pass it directly as knots to matern() (extent endpoints plus interior cell centers; see knots_from_spatraster()). (geostatsp / terra may mask matern, so call it with ::.)

Fit

geostatsp Matérn shape = 2 is ν = 2 (SPDE α = 3 in 2D), so we pass shape = 2 to matern() (cubic B-splines / random_fem_ssq_3). Locations come from the geometry column (WKT) in the model data frame.

knots_grid <- geostatsp::squareRaster(loaloa, cells = 6, buffer = 300 * 1000)
loaloa_df <- as.data.frame(loaloa_for_fit$data, geom = "WKT")

loaloa_model <- adlaplace::binomial(y, size = N) ~
  elev + evi +
  adlaplace::iid(villageID, init = 0.2) +
  adlaplaceGrf::matern(
    geometry,
    knots = knots_grid,
    shape = 2L,
    # range = practical rho = sqrt(8*nu)/kappa (m); sd = field SD
    init = c(100 * 1000, 0.2),
    lower = c(1 * 1000, 0.001)
  )

fit <- adlaplace::adlaplace(
  loaloa_model,
  data = loaloa_df,
  num_shards = 20L,
  control = list(maxit = 200, trace = 3, REPORT = 1),
  config = list(num_threads = 2L),
  verbose = FALSE
)
#> N = 6, M = 5 machine precision = 2.22045e-16
#> At X0, 0 variables are exactly at the bounds
#> At iterate     0  f=        11024  |proj g|=       1810.4
#> At iterate     1  f =        10802  |proj g|=        1704.3
#> At iterate     2  f =        10280  |proj g|=        880.56
#> At iterate     3  f =        10047  |proj g|=        323.65
#> At iterate     4  f =       9961.9  |proj g|=        121.73
#> At iterate     5  f =       9914.7  |proj g|=        90.542
#> At iterate     6  f =         9877  |proj g|=         41.35
#> At iterate     7  f =       9869.2  |proj g|=        16.928
#> At iterate     8  f =       9866.4  |proj g|=        11.085
#> At iterate     9  f =         9865  |proj g|=        2.8793
#> At iterate    10  f =       9864.8  |proj g|=        15.332
#> At iterate    11  f =       9864.5  |proj g|=        2.0543
#> At iterate    12  f =       9864.5  |proj g|=        2.3233
#> At iterate    13  f =       9864.3  |proj g|=        1.5686
#> At iterate    14  f =       9864.2  |proj g|=        2.0568
#> At iterate    15  f =       9863.9  |proj g|=        2.3491
#> At iterate    16  f =       9863.9  |proj g|=        3.4063
#> At iterate    17  f =       9863.8  |proj g|=        2.4097
#> At iterate    18  f =       9863.7  |proj g|=        4.3066
#> At iterate    19  f =       9863.7  |proj g|=        2.3779
#> At iterate    20  f =       9863.6  |proj g|=        3.0863
#> At iterate    21  f =       9863.3  |proj g|=        4.1074
#> At iterate    22  f =       9863.1  |proj g|=        7.4599
#> At iterate    23  f =         9863  |proj g|=        2.3925
#> At iterate    24  f =       9862.9  |proj g|=        3.0329
#> At iterate    25  f =       9862.9  |proj g|=       0.86045
#> At iterate    26  f =       9862.9  |proj g|=       0.70823
#> At iterate    27  f =       9862.9  |proj g|=       0.52402
#> At iterate    28  f =       9862.9  |proj g|=      0.027066
#> At iterate    29  f =       9862.9  |proj g|=        0.1095
#> At iterate    30  f =       9862.9  |proj g|=      0.036299
#> At iterate    31  f =       9862.9  |proj g|=      0.022447
#> 
#> iterations 31
#> function evaluations 37
#> segments explored during Cauchy searches 31
#> BFGS updates skipped 0
#> active bounds at final generalized Cauchy point 0
#> norm of the final projected gradient 0.0224474
#> final function value 9862.91
#> 
#> F = 9862.91
#> final  value 9862.910184 
#> converged
fit$optim$convergence
#> [1] 0
coef(fit)
#>               elev_linear                evi_linear                 intercept 
#>             -2.588382e-01              8.046617e-01             -2.351296e+00 
#>             villageID_iid geometry_matern_log_range    geometry_matern_log_sd 
#>              6.587140e-01              2.615201e+05              1.771369e+00

loaFit <- geostatsp::glgm(
  y ~ elev_s + evi_s +
    f(
      villageID,
      prior = "pc.prec",
      param = c(1, 0.5),
      model = "iid"
    ),
  loaloa,
  grid = 100,
  covariates = list(elev_s = elev_s_r, evi_s = evi_s_r),
  family = "binomial",
  Ntrials = loaloa$N,
  shape = 2,
  buffer = 25000,
  prior = list(sd = 1, range = 500 * 1000),
  control.inla = list(strategy = "gaussian")
)

Summary

logLik(fit)
#> 'log Lik.' -9862.91 (df=6)

cf <- coef(fit)
ad_pars <- c(
  intercept = unname(cf[["intercept"]]),
  elev = unname(cf[["elev_linear"]]),
  evi = unname(cf[["evi_linear"]]),
  village_sd = unname(cf[["villageID_iid"]]),
  # matern: practical range rho = sqrt(8*nu)/kappa (m); sd = field SD
  range = unname(cf[["geometry_matern_log_range"]]),
  sd = unname(cf[["geometry_matern_log_sd"]])
)
ad_pars
#>     intercept          elev           evi    village_sd         range 
#> -2.351296e+00 -2.588382e-01  8.046617e-01  6.587140e-01  2.615201e+05 
#>            sd 
#>  1.771369e+00

if (length(loaFit$parameters)) {
  glgm_sum <- loaFit$parameters$summary
  knitr::kable(glgm_sum[, c("mean", "0.025quant", "0.975quant")], digits = 3)

  # Map ad_pars names to rows of loaFit$parameters$summary
  row_for <- function(...) {
    pats <- c(...)
    hits <- unique(unlist(lapply(pats, function(p) {
      grep(p, rownames(glgm_sum), ignore.case = TRUE)
    })))
    if (!length(hits)) {
      return(NA_character_)
    }
    rownames(glgm_sum)[hits[1L]]
  }
  map <- c(
    intercept = row_for("^\\(Intercept\\)$", "^Intercept$"),
    elev = row_for("^elev"),
    evi = row_for("^evi"),
    village_sd = row_for("^sd villageID$", "sd.?villageID", "villageID.*[Ss]d"),
    range = row_for("^range/1000$", "^range$"),
    sd = row_for("^sd$", "SD for space")
  )

  glgm_pars <- setNames(
    vapply(map, function(nm) {
      if (is.na(nm)) {
        return(NA_real_)
      }
      as.numeric(glgm_sum[nm, "mean"])
    }, numeric(1)),
    names(map)
  )
  # summary may store range in km as "range/1000"
  if (!is.na(map[["range"]]) && identical(map[["range"]], "range/1000")) {
    glgm_pars[["range"]] <- glgm_pars[["range"]] * 1000
  }

  compare <- rbind(
    adlaplace = ad_pars[names(map)],
    glgm = glgm_pars
  )
  knitr::kable(compare, digits = 3)
}

Estimated field

my_sim <- matern_est(fit, loaloa_for_fit$grid, n = 10)

Posterior mean and SD of the spatial field from both fits (glgm when available):

data("worldMap", package = "mapmisc")
worldMap <- terra::unwrap(worldMap)
worldMap <- terra::crop(worldMap, terra::ext(c(-2, 10, -10, 5) * 1e6))
worldMapT <- terra::project(worldMap, loaloa)
mapmisc::map.new(loaloa)
if (plot_compare) {
  terra::plot(
    loaFit$raster[["random.mean"]],
    add = TRUE, col = my_col$col, breaks = my_col$breaks, legend = FALSE
  )
}
terra::plot(worldMapT, add = TRUE)
mapmisc::legendBreaks("right", my_col)

Posterior mean and SD of the spatial field (glgm vs adlaplace).


mapmisc::map.new(loaloa)
terra::plot(
  my_sim[["mean"]],
  add = TRUE, col = my_col$col, breaks = my_col$breaks, legend = FALSE
)
terra::plot(worldMapT, add = TRUE)
mapmisc::legendBreaks("right", my_col)

Posterior mean and SD of the spatial field (glgm vs adlaplace).


mapmisc::map.new(loaloa)
if (plot_compare) {
  terra::plot(
    loaFit$raster[["random.sd"]],
    add = TRUE, col = sd_col$col, breaks = sd_col$breaks, legend = FALSE
  )
}
terra::plot(worldMapT, add = TRUE)
mapmisc::legendBreaks("right", sd_col)

Posterior mean and SD of the spatial field (glgm vs adlaplace).


mapmisc::map.new(loaloa)
terra::plot(
  my_sim[["sd"]],
  add = TRUE, col = sd_col$col, breaks = sd_col$breaks, legend = FALSE
)
terra::plot(worldMapT, add = TRUE)
mapmisc::legendBreaks("right", sd_col)

Posterior mean and SD of the spatial field (glgm vs adlaplace).

Four posterior draws of the adlaplace Matérn field:

sim_layers <- paste0("sim", 1:6)
terra::plot(my_sim[[sim_layers]], nc = 2, col = my_col$col, breaks = my_col$breaks, legend = FALSE)
mapmisc::legendBreaks("right", my_col, bty = "n")

- Data
- Knot lines
- Fit
- Summary
- Estimated field