Loaloa prevalence with FEM Matérn

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::geostatData(
  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 hierarchical knot specification: a coarse grid over a region buffered by about 300 km (Neumann padding), then the study extent, then finer rectangles over the west and east village clusters. Pass the hb_knots() result as knots to matern().

loaloa_knots <- adlaplaceFem::hb_knots(
  outer = geostatsp::squareRaster(loaloa, cells = 4, buffer = 100 * 1000),
  inner = list(
    terra::extend(terra::ext(loaloa), 30 * 1000),
    list(
      terra::ext(3.9e5, 5.5e5, 5.4e5, 7.3e5),
      terra::ext(6e5, 8.1e5, 3.7e5, 7.4e5),
      terra::ext(1e6, 1.19e6, 400000, 600000)
    )
  ),
  fact = 2
)

knitr::kable(hb_summary(hb_basis(loaloa_knots, degree = 3L)), digits = 0)

	level	n_active	n_x	n_y	total_active
x	0	77	10	8	312
x1	1	85	20	14	312
x2	2	150	35	23	312

Hierarchical knot regions (coarse extent and west/east refinements) over village locations.

(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.

loaloa_df <- as.data.frame(loaloa_for_fit$data, geom = "WKT")

loaloa_model <- adlaplace::binom(y, size = N) ~
  elev + evi +
  adlaplace::iid(villageID, init = 0.2) +
  adlaplaceFem::matern(
    geometry,
    knots = loaloa_knots,
    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,
  control = list(maxit = 50, trace = 3, REPORT = 1),
  config = list(num_threads = 2L, num_shards = 20L),
  verbose = FALSE
)
#> N = 6, M = 5 machine precision = 2.22045e-16
#> At X0, 0 variables are exactly at the bounds
#> At iterate     0  f=        10811  |proj g|=       1304.5
#> At iterate     1  f =        10468  |proj g|=        1166.2
#> At iterate     2  f =        10090  |proj g|=        497.91
#> At iterate     3  f =       9969.4  |proj g|=        226.01
#> At iterate     4  f =       9902.4  |proj g|=        77.211
#> At iterate     5  f =       9873.2  |proj g|=        38.547
#> At iterate     6  f =       9861.6  |proj g|=         29.53
#> At iterate     7  f =       9858.2  |proj g|=        25.926
#> At iterate     8  f =       9857.3  |proj g|=        8.5979
#> At iterate     9  f =       9857.1  |proj g|=        4.0208
#> At iterate    10  f =       9856.8  |proj g|=        3.6259
#> At iterate    11  f =       9856.7  |proj g|=        4.1873
#> At iterate    12  f =       9856.2  |proj g|=        6.3412
#> At iterate    13  f =       9855.1  |proj g|=        5.8223
#> At iterate    14  f =       9854.4  |proj g|=        9.9145
#> At iterate    15  f =       9853.9  |proj g|=        5.0565
#> At iterate    16  f =       9853.8  |proj g|=        4.6368
#> At iterate    17  f =       9853.6  |proj g|=        1.6763
#> At iterate    18  f =       9853.5  |proj g|=       0.71721
#> At iterate    19  f =       9853.4  |proj g|=       0.31697
#> At iterate    20  f =       9853.4  |proj g|=       0.11512
#> At iterate    21  f =       9853.4  |proj g|=      0.069114
#> At iterate    22  f =       9853.4  |proj g|=      0.022564
#> 
#> iterations 22
#> function evaluations 26
#> segments explored during Cauchy searches 22
#> BFGS updates skipped 0
#> active bounds at final generalized Cauchy point 0
#> norm of the final projected gradient 0.0225638
#> final function value 9853.43
#> 
#> F = 9853.43
#> final  value 9853.426195 
#> converged
fit$optim$convergence
#> NULL
coef(fit)
#>               elev_linear                evi_linear                 intercept 
#>             -2.261879e-01              6.552641e-01             -2.241525e+00 
#>             villageID_iid geometry_matern_log_range    geometry_matern_log_sd 
#>              5.184559e-01              1.235380e+05              1.534549e+00

loaFit <- geostatsp::glgm(
  y ~ elev_s + evi_s +
    f(
      villageID,
      prior = "pc.prec",
      param = c(1, 0.5),
      model = "iid"
    ),
  loaloa,
  grid = 20,
  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.' -9853.426 (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.241525e+00 -2.261879e-01  6.552641e-01  5.184559e-01  1.235380e+05 
#>            sd 
#>  1.534549e+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)
}

	intercept	elev	evi	village_sd	range	sd
adlaplace	-2.242	-0.226	0.655	0.518	123538.0	1.535
glgm	-2.403	-0.140	0.590	0.562	185955.5	1.555

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):

opar <- par(mfrow = c(if (plot_compare) 2L else 1L, 2L))

# map.new() defaults to zero margins, so ask for a strip to hold the title.
# legendBreaks() restores par(no.readonly = TRUE) on exit, and re-applying
# mfrow moves the panel pointer to the last cell, which would send the next
# map.new() to a fresh page; save and restore mfg to keep the mfrow layout.
plot_field <- function(r, cols, main) {
  mapmisc::map.new(loaloa, mar = c(0, 0, 2, 0))
  terra::plot(
    r,
    add = TRUE, col = cols$col, breaks = cols$breaks, legend = FALSE
  )
  terra::plot(worldMapT, add = TRUE, border = "grey20", lwd = 0.7, col = NA)
  mfg <- par("mfg")
  mapmisc::legendBreaks("right", cols)
  par(mfg = mfg)
  title(main = main, line = 0.4, cex.main = 0.95)
}

if (plot_compare) {
  plot_field(loaFit$raster[["random.mean"]], my_col, "glgm mean")
}
plot_field(my_sim[["mean"]], my_col, "adlaplace mean")
if (plot_compare) {
  plot_field(loaFit$raster[["random.sd"]], sd_col, "glgm SD")
}
plot_field(my_sim[["sd"]], sd_col, "adlaplace SD")

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

par(opar)

Six posterior draws of the adlaplace Matérn field:

terra::plot(
  my_sim[[sim_layers]],
  nc = 2, col = sim_col$col, breaks = sim_col$breaks, legend = FALSE
)
mapmisc::legendBreaks("right", sim_col, bty = "n")

Loaloa prevalence with FEM Matérn

Patrick Brown

2026-10-11

Data

Knot lines

Fit

Summary

Estimated field