Tanzania SP-Resistance Marker BYM2 Models

Modeling notes for Main.R

1 Overview

We fit a separate Bayesian spatio-temporal BYM2 model (R-INLA) for each of the three SP-resistance markers, by region and year (2021-2023):

Marker	Gene	Columns used
K540E	pfdhps	found_540E / tested_540
A581G	pfdhps	found_581G / tested_581
I164L	pfdhfr	found_164L / tested_164

Covariates: Malaria Atlas Project (MAP) package pfpr (PfPR2-10, time-varying by region-year), sulfadoxine-pyrimethamine sp (SP_cov, % receiving 3+ IPTp-SP doses, region-level, from the 2022 report), MAP access (travel time to healthcare in minutes, region-level, log-transformed).

2 Why BYM2, and why the full grid

Mutation 581

Roughly half of Tanzania’s 26 regions were only surveyed for these markers in 2023, the 2021 and 2022 years have some missingness, not zero. The script builds the complete 26-region x 3-year grid on purpose and lets INLA produce a posterior prediction for the missing cells, borrowing strength from neighbouring regions (the spatial BYM2 term) and from each region’s own temporal trend and covariates. Regions/years with no survey are marked observed = FALSE in the output and shown with a dashed border on the maps.

This works because pfpr/access/SP have complete 26-region x year values, rather than the mutations data which only have rows for the 52 region-years that were actually surveyed: 26 in 2023, 13 in 2022 and 13 in 2021. INLA requires non-NA covariates even for “prediction” rows with an NA response, so that complete secondary source is what makes the full-grid prediction possible.

3 Inputs

Expected in ~/Desktop/Nicola_DM4DM/:

• Shapefile: region polygons.
• Covariates: region x year PfPR + access + SP (complete for all 26 regions).
• Mutations data — the joined marker-count.

4 Pipeline

1. Loads and harmonises region names across all three files (join key = lower-cased, letters-only, so “Dar-es-salaam” / “Dar-Es-Salaam” match).

2. Builds a queen-contiguity adjacency graph (spdep::poly2nb) and the INLA graph file needed for the BYM2 spatial term.

3. Builds the full 26 x 3 region-year panel and standardises covariates.

4. Runs exploratory spatial analysis: global Moran’s I on 2023 raw prevalence (the one year with full mutation values), a local Moran’s I (LISA) cluster map, and choropleths of raw prevalence.

5. Fits, per marker, three binomial models for comparison: BYM2 (spatial + rw1 temporal), IID-only (no spatial structure), and covariates-only — with a DIC/WAIC/mean-log-CPO table so you can see how much the spatial term actually buys you.

6. Reports each model’s phi (share of the spatial random effect that is spatially structured vs. unstructured noise) — the key BYM2 diagnostic.

7. Extracts fitted/predicted prevalence with 95% credible intervals and exceedance probabilities for every region-year-marker cell, including the unobserved ones. The prevalence cutoff is marker-specific (exceed_thresholds), not one blanket value:

   Marker	Threshold	Why
   K540E (pfdhps)	95%	Half of WHO’s cited IPTp-SP-reconsideration criterion (K540E > 95% AND A581G > 10%)
   A581G (pfdhps)	10%	The other half of that same WHO criterion
   I164L (pfdhfr)	5%	No standalone WHO cutoff exists; still rare in East Africa, so a low bar flags emergence early

   A region-year is flagged likely_exceeds = TRUE (shown with a red outline on the exceedance maps) when the model is more than 80% confident (LIKELY_EXCEEDS_PROB) the true prevalence is over its marker’s threshold — a common convention in Bayesian disease mapping, not a hard rule. These thresholds are defensible starting points, not settled policy — revisit with your collaborator if better-supported cutoffs exist, especially for I164L.

8. Tests whether the covariates-only model still has significant residual spatial autocorrelation (Moran’s I on residuals) — direct evidence for whether the spatial term is earning its keep.

9. Writes maps (raw prevalence, fitted/predicted prevalence, spatial random effect, exceedance probability, LISA) and tables (CSV), plus a geopackage with the wide predictions attached to the region polygons.

10. Held-out validation using 2023. Since 2023 is the only year every region was surveyed, it doubles as a validation set for the same borrowing-strength mechanism that fills in 2021/2022 for the 13 partially-surveyed regions (which have no ground truth to check against at all). Each marker’s model is refit with the 2023 counts hidden, and the resulting purely spatio-temporally predicted 2023 prevalence is compared against what was actually observed: a scatter plot (maps/validation_2023_observed_vs_predicted.png) plus correlation, RMSE, MAE, and 95% credible-interval coverage (tables/validation_2023_heldout_metrics.csv). Coverage near 0.95 means the intervals are well-calibrated; much lower means the model is overconfident.

5 Assumptions worth double-checking with your collaborator

• Separate model per marker rather than a joint multivariate model. Simpler and standard for this kind of surveillance mapping; a joint model (sharing information across markers) is a reasonable extension if the markers turn out to be strongly correlated.
• Temporal structure: a rw1 random walk over the 3 years (mild smoothing) rather than a plain iid year effect or a full space-time interaction. With only 3 time points the difference is small; a Type I interaction term (space_time_id, already built into the panel) is a drop-in extension if you want it.
• sp and access are treated as region-level constants repeated across 2021-2023, because the source data only has one value per region for each (SP coverage: 2022 report; access: 2019 raster).
• Priors: PC-priors, P(phi < 0.5) = 0.5 and P(sigma > 1) = 0.01 — the standard defaults recommended for BYM2; tune if you have stronger prior beliefs.

6 Running it

This script was written and reviewed without a local R/INLA runtime available in the authoring session, so run it section by section the first time (especially the package-install section and the adjacency-graph section) and check for warnings before letting the model-fitting loop run end to end.

# Run the full pipeline from this notes file if you'd rather work from Quarto:
source("Main.R")

7 Outputs

After a successful run, BYM2_outputs/ will contain:

• maps/ — raw & fitted/predicted choropleths per marker/year, spatial random effect maps, LISA cluster map, exceedance probability maps
• tables/ — model comparison (DIC/WAIC/CPO), fixed effects, phi, Moran’s I (raw and residual)
• Tanzania_BYM2_predictions.csv and Tanzania_BYM2_predictions.gpkg — region x year x marker predictions, including the spatio-temporally predicted cells for regions with no direct survey in 2021/2022