---
title: "Anthrax in Africa: a systematic review and meta-analysis of ecological niche and species distribution model performance, and human behavioural risk factors"
author: "Marta M. Dolcet-Negre -- Universidad de Salamanca"
date: "`r format(Sys.Date(), '%d %B %Y')`"
output:
  pdf_document:
    toc: true
    number_sections: true
    fig_width: 9
    fig_height: 6
    latex_engine: pdflatex
    keep_tex: false
geometry: margin=2.5cm
bibliography: references.bib
csl: vancouver-brackets.csl
link-citations: true
header-includes:
  - \usepackage[utf8]{inputenc}
  - \usepackage[T1]{fontenc}
  - \usepackage[english]{babel}
  - \usepackage{booktabs}
  - \usepackage{longtable}
  - \usepackage{array}
  - \usepackage{float}
  - \usepackage{caption}
  - \usepackage{pdflscape}
  - \captionsetup[table]{skip=5pt}
  - \renewcommand{\arraystretch}{1.2}
---

```{r setup, include=FALSE}
# Create figures folder
if (!dir.exists("figuras")) dir.create("figuras")

knitr::opts_chunk$set(
  echo = FALSE,
  warning = FALSE,
  message = FALSE,
  fig.align = "center",
  dpi = 600,
  fig.pos = "H",
  out.width = "95%",
  fig.path = "figuras/",
  dev = c("pdf", "png")
)
options(knitr.kable.NA = "--")
```

# Framework and rationale

## Context

This document executes the **partial quantitative synthesis** outlined in
Section 10.2 of the master protocol for the systematic review. The protocol
establishes that, given the heterogeneity of designs, outcomes, and metrics, the
primary synthesis is **narrative and methodological**. However, where sufficiently
comparable performance metrics exist, a quantitative pooling with heterogeneity
assessment is conducted. Two such syntheses are presented here:

1. **Model discrimination (AUC)** for *ecological niche modelling* and
   *species distribution modelling* (Family A) studies.
2. **Adjusted odds ratios (aOR)** for human-anthrax risk factors reported by the
   analytical outbreak-epidemiology (Family C) studies.

## Scope

- **Corpus**: 33 studies included after title/abstract and full-text screening
  (Batch 1, 63 deduplicated records; screening 1 April 2026, audited 30 May 2026).
- **Quantitative AUC pooling**: restricted to Family A studies
  (ENM/SDM/suitability) reporting an AUC with an estimable standard error.
- **Quantitative OR pooling**: restricted to Family C exposures reported by
  $\geq 3$ studies on a comparable individual, disease-outcome basis.
- **Complementary syntheses**: Families B (Bayesian/INLA),
  D (hotspot/operational mapping), and E (climate projections) are synthesised
  narratively and with descriptive tables.

## Dataset status

> **Extraction completed**: AUC and SE values were extracted from the full texts
> of the 19 Family A studies (see `family_A_metrics.xlsx`). Unreported SEs were
> approximated using the Hanley & McNeil [@hanley1982] formula or, for one study, from the reported
> range. Test AUC was prioritised over training AUC when both were available, and
> the value used for pooling was reconciled against the curated extraction sheet.
> Family C adjusted odds ratios were extracted from `family_C_OR_extraction.xlsx`.

# Packages

All analyses were performed in R. The pooled estimates and forest/funnel plots
were produced with the **meta** package [@balduzzi2019]; the mixed-effects
meta-regressions were fitted with **metafor** [@viechtbauer2010]. The full set of
packages used is listed below.

```{r packages, echo=TRUE, results='hide', message=FALSE, warning=FALSE}
pkgs <- c("meta", "metafor", "ggplot2", "dplyr", "tidyr",
          "knitr", "kableExtra", "forcats", "irr", "psych")

for (pkg in pkgs) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    install.packages(pkg, quiet = TRUE, repos = "https://cloud.r-project.org")
  }
}

library(meta)
library(metafor)
library(ggplot2)
library(dplyr)
library(tidyr)
library(knitr)
library(kableExtra)
library(forcats)
```

# Dataset: 33 included studies (Batch 1)

## Dataset construction

The 33 studies come from the title/abstract and full-text screening documented in
`screening_results_MMDN.md` and `screening_title_abstract.xlsx`, with the
methodological-family classification taken from `clasificacion-familias.docx`
and the AUC extraction from `family_A_metrics.xlsx`. Each record is identified by
its screening number (`rec_id`) for traceability.

```{r dataset}
# ===================================================================
# MAIN DATASET — 33 included studies (Batch 1, T/A + full-text)
# ===================================================================
# Classification by families A–E according to master protocol (Section 8)
# and clasificacion-familias.docx.
#
# CONVENTIONS:
#   auc      = NA if not reported or not applicable to the design
#   se_auc   = NA if not estimable
#   auc_verified = TRUE if the value appears explicitly in the full text
#   n_obs    = NA if not reported
# ===================================================================

dat <- data.frame(
  rec_id = c(5, 8, 9, 11, 19, 21, 25, 26, 28, 29, 30, 31, 32, 33,
             37, 38, 39, 40, 41, 42, 43, 44, 45, 47, 48, 50, 51, 52,
             54, 57, 58, 61, 63),

  author = c(
    "Witt et al.",                   # R5
    "Mutesi et al.",                 # R8
    "Tumusiime et al.",              # R9
    "Kwizera et al.",                # R11
    "Kimaro et al.",                 # R19
    "John et al.",                   # R21
    "Dougherty et al.",              # R25
    "Gachohi et al.",                # R26
    "Ndolo et al. (Kenya)",          # R28
    "Hamutyinei Dhliwayo et al.",    # R29  (compound surname; corrected)
    "Deka et al.",                   # R30
    "Aminu et al.",                  # R31
    "Pittiglio et al.",              # R32
    "Ndolo et al. (Uganda)",         # R33
    "Makurumidze et al.",            # R37
    "Otieno et al. (b)",             # R38
    "Otieno et al. (a)",             # R39
    "Migisha et al.",                # R40
    "Nderitu et al.",                # R41
    "Assefa et al.",                 # R42
    "Driciru et al.",                # R43
    "Nakanwagi et al.",              # R44
    "Romero-Alvarez et al.",         # R45
    "Yousuf et al.",                 # R47
    "Mwakapeje et al.",              # R48
    "Kamboyi et al.",                # R50
    "Carlson et al.",                # R51
    "Steenkamp et al.",              # R52
    "Driciru et al.",                # R54
    "Kracalik et al.",               # R57
    "Lehman et al.",                 # R58
    "Nsoh et al.",                   # R61 (published as Nsoh AE et al., PAMJ 2016)
    "Blackburn et al."               # R63
  ),

  year = c(2026, 2025, 2025, 2025, 2024, 2024, 2022, 2022, 2022, 2022,
           2022, 2022, 2022, 2022, 2021, 2021, 2021, 2021, 2021, 2020,
           2020, 2020, 2020, 2020, 2019, 2019, 2019, 2018, 2018, 2017,
           2017, 2016, 2015),

  country = c(
    "Ethiopia",            # R5
    "Uganda",              # R8
    "Uganda",              # R9
    "Uganda",              # R11
    "Tanzania",            # R19
    "Zimbabwe",            # R21
    "Namibia",             # R25
    "Kenya",               # R26
    "Kenya",               # R28
    "Zimbabwe",            # R29
    "Global/Africa",       # R30
    "Tanzania",            # R31
    "West Africa",         # R32
    "Uganda",              # R33
    "Zimbabwe",            # R37
    "Kenya",               # R38
    "Kenya",               # R39
    "Uganda",              # R40
    "Kenya",               # R41
    "Ethiopia",            # R42
    "Uganda",              # R43
    "Uganda",              # R44
    "West/Central Africa", # R45
    "Ethiopia",            # R47
    "Tanzania",            # R48
    "Zambia",              # R50
    "Global/Africa",       # R51
    "South Africa",        # R52
    "Uganda",              # R54
    "Ghana",               # R57
    "Zambia",              # R58
    "Ghana",               # R61
    "Nigeria/Cameroon/Chad" # R63
  ),

  region = c(
    "East",        "East",        "East",        "East",
    "East",        "Southern",    "Southern",    "East",
    "East",        "Southern",    "Pan-African", "East",
    "West",        "East",        "Southern",    "East",
    "East",        "East",        "East",        "East",
    "East",        "East",        "West",        "East",
    "East",        "Southern",    "Pan-African", "Southern",
    "East",        "West",        "Southern",    "West",
    "West"
  ),

  # ---------------------------------------------------------------
  # Methodological families (protocol Section 8)
  #   A = ENM/SDM/suitability
  #   B = Bayesian/INLA/spatiotemporal
  #   C = Case-control / analytical outbreak epidemiology
  #   D = Hotspot / operational risk mapping
  #   E = Climate projections (always shared with A)
  # ---------------------------------------------------------------
  family = c(
    "D",     # R5:  Getis-Ord Gi*
    "C",     # R8:  case-control outbreak
    "C",     # R9:  case-control outbreak
    "C",     # R11: case-control outbreak
    "A",     # R19: MaxEnt
    "A+E",   # R21: ensemble SDM + climate
    "A",     # R25: MaxEnt + movement
    "A+D",   # R26: BRT + KDE + LISA
    "A+B",   # R28: ENM component + Bayesian INLA
    "C",     # R29: case-control
    "A",     # R30: ensemble ENM global
    "B+D",   # R31: Bayesian spatial GLMM + mapping
    "A",     # R32: MaxEnt + seasonality
    "A+B",   # R33: ENM component + Bayesian INLA
    "C",     # R37: case-control
    "A+E",   # R38: BRT + RCP
    "A",     # R39: BRT
    "C",     # R40: case-control
    "D",     # R41: SaTScan clustering
    "A",     # R42: ensemble biomod2
    "A",     # R43: MaxEnt
    "C",     # R44: retrospective cohort + case-control
    "A",     # R45: MaxEnt (kuenm)
    "A+E",   # R47: MaxEnt + climate (preprint)
    "A",     # R48: MaxEnt
    "C+D",   # R50: logistic regression + mapping
    "A",     # R51: BRT global
    "A",     # R52: MaxEnt
    "D",     # R54: SaTScan spatiotemporal
    "A",     # R57: Random Forest ENM
    "C",     # R58: logistic regression outbreak
    "D",     # R61: kriging
    "A"      # R63: GARP ENM
  ),

  # Specific method
  method = c(
    "Getis-Ord Gi*",                          # R5
    "Case-control + logistic regression",      # R8
    "Case-control + logistic regression",      # R9
    "Case-control + logistic regression",      # R11
    "ENM (MaxEnt)",                            # R19
    "Ensemble SDM (biomod2) + BCC-CSM",        # R21
    "MaxEnt + movement (SSF/HMM)",             # R25
    "BRT + KDE + LISA",                        # R26
    "Bayesian INLA (ZIP/GAM) + ENM",           # R28
    "Case-control 1:2 + logistic reg.",        # R29
    "Ensemble ENM (biomod2)",                  # R30
    "Bayesian spatial GLMM (glmmfields)",      # R31
    "MaxEnt + climate seasonality",            # R32
    "Bayesian INLA GAM/SPDE + ENM",            # R33
    "Case-control 1:1 + logistic reg.",        # R37
    "BRT + RCP 4.5/8.5",                       # R38
    "BRT",                                     # R39
    "Case-control unmatched + logistic",       # R40
    "SaTScan + agro-ecological",               # R41
    "Ensemble (biomod2: 9 algorithms)",        # R42
    "MaxEnt",                                  # R43
    "Retrospective cohort + case-control",     # R44
    "MaxEnt (kuenm)",                          # R45
    "MaxEnt + climate projections",            # R47
    "MaxEnt",                                  # R48
    "Multivariable logistic + mapping",        # R50
    "BRT (global)",                            # R51
    "MaxEnt",                                  # R52
    "SaTScan spatial scan statistics",         # R54
    "Random Forest ensemble",                  # R57
    "Cross-sectional + logistic reg.",         # R58
    "Kriging interpolation",                   # R61
    "GARP ecological niche model"              # R63
  ),

  # ---------------------------------------------------------------
  # AUC — from family_A_metrics.xlsx
  # Value used for pooling. NA = not applicable / not reported in a
  # comparable format (pROC, sensitivity, no explicit AUC).
  # ---------------------------------------------------------------
  auc = c(
    NA,      # R5:  hotspot, no AUC
    NA,      # R8:  case-control (OR)
    NA,      # R9:  case-control (OR)
    NA,      # R11: case-control (OR)
    0.912,   # R19: mean test AUC = 0.912 +/- 0.026
    NA,      # R21: only ROC>0.7, no exact AUC
    NA,      # R25: no evaluation metric reported
    NA,      # R26: AUC reused from R38/R39 (not independent)
    NA,      # R28: Bayesian INLA, sensitivity 75% (no AUC)
    NA,      # R29: case-control (OR)
    NA,      # R30: pROC ratio = 1.57, not comparable
    NA,      # R31: Bayesian GLMM
    0.93,    # R32: training AUC = 0.93 (SD 0.01); no separate test AUC
    NA,      # R33: Bayesian INLA, sensitivity 90% (no AUC)
    NA,      # R37: case-control (OR)
    0.929,   # R38: mean test AUC = 0.929 (+/- SE reported)
    0.80,    # R39: mean ensemble AUC = 0.80
    NA,      # R40: case-control (OR)
    NA,      # R41: clustering, no AUC
    0.991,   # R42: ensemble ROC = 0.991 (TSS=0.926, Kappa=0.857)
    0.936,   # R43: mean test AUC = 0.936 +/- 0.015
    NA,      # R44: cohort/case-control (AR/RR)
    NA,      # R45: pROC + omission + AICc, no standard AUC
    0.75,    # R47: average test AUC = 0.75 (preprint)
    0.93,    # R48: AUC = 0.93
    NA,      # R50: logistic (OR), no AUC
    0.9244,  # R51: mean submodel AUC = 0.9244
    0.857,   # R52: AUC = 0.857 (Table 2 overall metric; see Limitations)
    NA,      # R54: SaTScan, no AUC
    0.87,    # R57: external test AUC = 0.87 (95% CI 0.85-0.90)
    NA,      # R58: case-control (OR)
    NA,      # R61: kriging, no AUC
    0.71     # R63: average AUC = 0.71 (range 0.57-0.86)
  ),

  # AUC SE — reported, from 95% CI, Hanley-McNeil (1982), or range
  se_auc = c(
    NA, NA, NA, NA,
    0.026,   # R19: from 95% CI width
    NA, NA, NA, NA, NA, NA, NA,
    0.010,   # R32: SD across replicates (training)
    NA, NA,
    0.004,   # R38: reported +/- SE
    0.048,   # R39: Hanley-McNeil
    NA,
    NA,      # R41
    0.008,   # R42: Hanley-McNeil
    0.015,   # R43: from 95% CI
    NA, NA,
    0.034,   # R47: Hanley-McNeil
    0.038,   # R48: Hanley-McNeil
    NA,
    0.004,   # R51: Hanley-McNeil
    0.024,   # R52: Hanley-McNeil (on AUC = 0.857)
    NA,
    0.013,   # R57: from 95% CI
    NA, NA,
    0.074    # R63: approximated from reported range 0.57-0.86
  ),

  # AUC verification (TRUE = confirmed in full text)
  auc_verified = c(
    NA, NA, NA, NA,
    TRUE,    # R19
    NA, NA, NA, NA, NA, NA, NA,
    TRUE,    # R32
    NA, NA,
    TRUE,    # R38
    TRUE,    # R39
    NA, NA,
    TRUE,    # R42
    TRUE,    # R43
    NA, NA,
    TRUE,    # R47
    TRUE,    # R48
    NA,
    TRUE,    # R51
    TRUE,    # R52
    NA,
    TRUE,    # R57
    NA, NA,
    TRUE     # R63
  ),

  # Sample size (occurrence points / presences used for the model)
  n_obs = c(
    NA, NA, NA, NA,
    636,     # R19
    112,     # R21
    40,      # R25
    69,      # R26
    95,      # R28
    NA,      # R29
    874,     # R30
    NA,      # R31 (5000 grid points; not an occurrence n)
    59,      # R32
    294,     # R33
    NA,      # R37
    178,     # R38
    69,      # R39
    NA,      # R40
    NA,      # R41
    128,     # R42
    471,     # R43
    NA,      # R44
    57,      # R45
    158,     # R47
    44,      # R48
    NA,      # R50
    5018,    # R51
    219,     # R52
    NA,      # R54
    67,      # R57
    NA,      # R58
    NA,      # R61
    NA       # R63 (village-level, n unclear)
  ),

  # DOI
  doi = c(
    "10.1016/j.sste.2026.100801",       # R5
    "10.21203/rs.3.rs-6854973/v1",        # R8  preprint
    "10.1186/s42522-025-00151-x",         # R9
    "10.1371/journal.pone.0336769",       # R11
    "10.1111/aje.13350",                  # R19
    "10.1186/s12889-024-17856-9",         # R21
    "10.1186/s40462-022-00331-8",         # R25
    "10.1038/s41598-022-24000-3",         # R26
    "10.1038/s41598-022-24589-5",         # R28
    "10.1371/journal.pone.0278537",       # R29
    "10.1111/tbed.14602",                 # R30
    "10.1038/s41598-022-14081-5",         # R31
    "10.3390/ani12091146",                # R32
    "10.1038/s41598-022-24281-8",         # R33
    "10.1186/s12889-021-10275-0",         # R37
    "10.3390/ijerph18084176",             # R38
    "10.1371/journal.pntd.0009301",       # R39
    "10.11604/pamj-oh.2021.5.13.29385",   # R40
    "10.1186/s12879-021-05871-9",         # R41
    "10.1016/j.prevetmed.2020.105155",    # R42
    "10.1371/journal.pone.0237223",       # R43
    "10.1371/journal.pntd.0008026",       # R44
    "10.1371/journal.pntd.0008131",       # R45
    "10.1101/2020.11.27.400879",          # R47 preprint
    "10.1016/j.ijid.2018.11.367",         # R48
    "10.1002/vms3.168",                   # R50
    "10.1038/s41564-019-0435-4",          # R51
    "10.1371/journal.pone.0191704",       # R52
    "10.1371/journal.pone.0206922",       # R54
    "10.1371/journal.pntd.0005885",       # R57
    "10.3201/eid2309.161597",             # R58
    "10.11604/pamj.supp.2016.25.1.6205",  # R61
    "10.1371/journal.pntd.0003931"        # R63
  ),

  stringsAsFactors = FALSE
)

# ---------------------------------------------------------------
# Derived columns
# ---------------------------------------------------------------

# Primary family (first letter before +)
dat$family_primary <- sub("\\+.*", "", dat$family)

# Includes an ENM/SDM (Family A) component?
dat$has_enm <- grepl("A", dat$family)

# Is it a case-control / outbreak study (Family C)?
dat$is_casecontrol <- grepl("C", dat$family)

# Label for plots
dat$label <- paste0(dat$author, " (", dat$year, ")")

cat("Dataset constructed: ", nrow(dat), " studies\n")
cat("  Family A (ENM/SDM):           ", sum(dat$has_enm), "\n")
cat("  Family B (Bayesian/INLA):     ", sum(grepl("B", dat$family)), "\n")
cat("  Family C (case-control):      ", sum(dat$is_casecontrol), "\n")
cat("  Family D (hotspot/mapping):   ", sum(grepl("D", dat$family)), "\n")
cat("  Family E (climate projection):", sum(grepl("E", dat$family)), "\n")
cat("  With usable AUC + SE:         ", sum(dat$has_enm & !is.na(dat$auc) & !is.na(dat$se_auc)), "\n")
```

## Dataset summary

```{r dataset-table}
dat %>%
  select(rec_id, author, year, country, family, method, auc) %>%
  mutate(
    auc = ifelse(is.na(auc), "--", sprintf("%.3f", auc)),
    Status = case_when(
      grepl("A", dat$family) & !is.na(dat$auc) & !is.na(dat$se_auc) ~ "Poolable",
      grepl("A", dat$family) & !is.na(dat$auc) ~ "AUC only",
      grepl("A", dat$family) ~ "No comparable AUC",
      TRUE ~ ""
    )
  ) %>%
  kable(
    col.names = c("Rec", "Author", "Year", "Country/Region",
                  "Family", "Method", "AUC", "Status"),
    caption = "Included studies (Batch 1, 33 studies). Family A: AUC status for pooling.",
    booktabs = TRUE, longtable = TRUE
  ) %>%
  kable_styling(latex_options = c("striped", "hold_position", "repeat_header"),
                font_size = 8) %>%
  kableExtra::landscape()
```

# Descriptive overview

## Distribution by methodological family

```{r descriptive-family, fig.height=5, fig.cap="Distribution of the 33 included studies by methodological family (A: ENM/SDM; B: Bayesian/INLA; C: outbreak analytical; D: hotspot/mapping; E: climate projection). Counts use non-exclusive family membership, so studies tagged with more than one family are counted in each."}
# Expand shared families
family_expanded <- dat %>%
  mutate(families = strsplit(family, "\\+")) %>%
  unnest(families) %>%
  mutate(
    family_name = case_when(
      families == "A" ~ "A: ENM/SDM",
      families == "B" ~ "B: Bayesian/INLA",
      families == "C" ~ "C: Outbreak analytical",
      families == "D" ~ "D: Hotspot/mapping",
      families == "E" ~ "E: Climate projection",
      TRUE ~ families
    )
  )

ggplot(family_expanded, aes(x = fct_infreq(family_name))) +
  geom_bar(fill = "steelblue", alpha = 0.8, width = 0.6) +
  geom_text(stat = "count", aes(label = after_stat(count)),
            vjust = -0.5, size = 4) +
  labs(x = NULL, y = "Number of studies") +
  theme_minimal(base_size = 13) +
  theme(axis.text.x = element_text(angle = 15, hjust = 1))
```

Ecological niche modelling and species distribution modelling (Family A) constituted the dominant analytical approach, comprising 19 of the 33 included studies once the corpus is counted with non-exclusive family membership (R28 and R33 [@ndolo2022kenya; @ndolo2022uganda], the two INLA models, retain an ENM component and are therefore counted in both A and B; R31 is now classified outside A as a Bayesian spatial GLMM). This finding reflects the strong influence of correlative presence-only modelling frameworks in the African anthrax literature, driven largely by the availability of MaxEnt and boosted regression tree (BRT) algorithms and by the practical appeal of mapping environmental suitability from opportunistic occurrence data. Analytical outbreak epidemiology (Family C, *K* = 9) represented the second most frequent approach, indicating that case-control designs remain a central tool for investigating individual-level risk factors during anthrax events. Hotspot and operational risk mapping (Family D, *K* = 7), Bayesian spatiotemporal models (Family B, *K* = 3), and climate-change projection studies (Family E, *K* = 3, all shared with Family A) were less represented. The scarcity of Bayesian hierarchical approaches is particularly notable given their capacity for formally integrating spatial dependence, temporal dynamics, and covariate uncertainty -- capabilities that remain central to the doctoral thesis framework.

## Geographic distribution

```{r descriptive-geo, fig.height=5, fig.cap="Geographic distribution of the 33 included studies by African region, coloured by primary methodological family."}
ggplot(dat, aes(x = fct_infreq(region), fill = family_primary)) +
  geom_bar(alpha = 0.85, width = 0.6) +
  scale_fill_brewer(palette = "Set2", name = "Family") +
  labs(x = NULL, y = "Number of studies") +
  theme_minimal(base_size = 13)
```

The geographic distribution of included studies reveals a pronounced concentration in Eastern Africa, which accounts for over half of the corpus. Kenya, Uganda, Tanzania, and Ethiopia together contribute the majority of quantitative anthrax studies, reflecting both the ecological significance of the East African Rift Valley system -- where alkaline soils and seasonal flooding create conditions favourable for *Bacillus anthracis* spore persistence -- and the comparatively stronger research infrastructure and surveillance capacity in these countries. Southern Africa (South Africa, Zimbabwe, Namibia, Zambia) is represented primarily through ecological modelling and case-control studies, while West Africa (Ghana, Nigeria) and pan-African or global-scale studies remain underrepresented. This geographic imbalance has important implications for the generalisability of pooled estimates, as environmental drivers of anthrax suitability may differ substantially between the humid equatorial landscapes of East Africa and the semi-arid Sahelian zone of West Africa.

## Temporal distribution

```{r descriptive-time, fig.height=5, fig.cap="Temporal distribution of the 33 included studies by publication year, coloured by primary methodological family."}
ggplot(dat, aes(x = year, fill = family_primary)) +
  geom_bar(alpha = 0.85, width = 0.7) +
  scale_fill_brewer(palette = "Set2", name = "Family") +
  scale_x_continuous(breaks = 2015:2026) +
  labs(x = "Publication year", y = "Number of studies") +
  theme_minimal(base_size = 13)
```

Publication activity in quantitative anthrax modelling in Africa has increased markedly since 2017, with the highest concentration of studies appearing between 2020 and 2022. This temporal pattern likely reflects the combined effect of improved access to remote sensing data (e.g., MODIS, Sentinel), the growing availability of open-source modelling platforms (MaxEnt, biomod2, R-INLA), and increased investment in One Health surveillance frameworks across sub-Saharan Africa. The relative scarcity of publications before 2017 -- with only two pioneering studies [@blackburn2015; @nsoh2016] -- underscores how recent the systematic application of geospatial and statistical modelling to African anthrax ecology truly is. The post-2020 surge includes the first Bayesian spatiotemporal analyses (Family B) and climate projection studies (Family E), suggesting a methodological diversification that may deepen further in the coming years.

# Partial quantitative synthesis: AUC pooling (Family A)

## Data preparation

```{r auc-prep}
# Studies eligible for AUC pooling: Family A, AUC reported, SE estimable.
# Because the family classification already excludes the Bayesian INLA models
# (R28/R33 not poolable as AUC), the AUC+SE filter selects
# the 12 poolable studies directly.
dat_auc <- dat %>%
  filter(has_enm, !is.na(auc), !is.na(se_auc)) %>%
  mutate(
    # Logit transformation to keep CIs within (0,1) and stabilise variance
    logit_auc = log(auc / (1 - auc)),
    # SE on the logit scale (delta method)
    se_logit  = se_auc / (auc * (1 - auc)),
    method_type = case_when(
      grepl("MaxEnt", method)        ~ "MaxEnt",
      grepl("BRT|Boosted", method)   ~ "BRT",
      grepl("Random Forest", method) ~ "RF",
      grepl("[Ee]nsemble|biomod", method) ~ "Ensemble",
      grepl("GARP", method)          ~ "GARP",
      TRUE                           ~ "Other"
    )
  )

cat("Studies eligible for AUC pooling:", nrow(dat_auc), "\n")
cat("Family A studies with AUC but without estimable SE:",
    sum(dat$has_enm & !is.na(dat$auc) & is.na(dat$se_auc)), "\n")
cat("Family A studies without a comparable AUC:",
    sum(dat$has_enm & is.na(dat$auc)), "\n\n")

if (nrow(dat_auc) >= 3) {
  cat("[OK] Sufficient data for pooling (K =", nrow(dat_auc), ")\n\n")
  cat("Studies included in pooling:\n")
  for (i in 1:nrow(dat_auc)) {
    cat(sprintf("  R%-3d %-28s AUC=%.3f  SE=%.4f  n=%s  %s\n",
                dat_auc$rec_id[i], dat_auc$label[i],
                dat_auc$auc[i], dat_auc$se_auc[i],
                ifelse(is.na(dat_auc$n_obs[i]), "NR",
                       as.character(dat_auc$n_obs[i])),
                dat_auc$method_type[i]))
  }
}
```

## Studies included in the pooling

```{r auc-table, results='asis'}
if (nrow(dat_auc) >= 3) {
  dat_auc %>%
    arrange(year) %>%
    select(rec_id, label, method_type, auc, se_auc, n_obs, region) %>%
    mutate(n_obs = ifelse(is.na(n_obs), "NR", as.character(n_obs))) %>%
    kable(
      col.names = c("Rec", "Study", "Method", "AUC", "SE", "n", "Region"),
      caption = "Family A studies included in the AUC pooling (K = 12)",
      digits = 4, booktabs = TRUE
    ) %>%
    kable_styling(latex_options = c("striped", "hold_position"), font_size = 9)
}
```

```{r meta-analysis-main}
# ===================================================================
# Random-effects meta-analysis (DerSimonian-Laird), logit(AUC) scale.
# Matches family_A_metrics.xlsx sheet 'Meta_analysis_results'.
# ===================================================================
m_auc <- metagen(
  TE      = logit_auc,
  seTE    = se_logit,
  studlab = label,
  data    = dat_auc,
  sm      = "MRAW",
  method.tau = "DL",      # DerSimonian-Laird (protocol)
  prediction = TRUE,
  title   = "AUC Pooling (Family A)"
)

# Back-transform logit -> AUC
pooled_auc   <- plogis(m_auc$TE.random)
pooled_lower <- plogis(m_auc$lower.random)
pooled_upper <- plogis(m_auc$upper.random)
pred_lower   <- plogis(m_auc$lower.predict)
pred_upper   <- plogis(m_auc$upper.predict)

cat("===================================================\n")
cat("AUC META-ANALYSIS -- Family A (ENM/SDM)\n")
cat("===================================================\n")
cat(sprintf("Studies included (K):         %d\n", m_auc$k))
cat(sprintf("Pooled AUC:                   %.3f\n", pooled_auc))
cat(sprintf("95%% CI:                       %.3f - %.3f\n",
            pooled_lower, pooled_upper))
cat(sprintf("95%% prediction interval:      %.3f - %.3f\n",
            pred_lower, pred_upper))
cat(sprintf("I^2 (heterogeneity):          %.1f%%\n", m_auc$I2 * 100))
cat(sprintf("tau^2 (DL):                   %.4f\n", m_auc$tau2))
cat(sprintf("Q test: Q = %.2f, df = %d, p = %.4f\n",
            m_auc$Q, m_auc$df.Q, m_auc$pval.Q))

# --- Raw-scale sensitivity (no transformation) ---
m_auc_raw <- metagen(TE = auc, seTE = se_auc, studlab = label,
                     data = dat_auc, sm = "MRAW", method.tau = "DL")
cat(sprintf("\nRaw-scale sensitivity pooled AUC: %.3f (95%% CI %.3f - %.3f), I^2 = %.0f%%\n",
            m_auc_raw$TE.random, m_auc_raw$lower.random,
            m_auc_raw$upper.random, m_auc_raw$I2 * 100))
cat("(Raw scale produces upper limits > 1 for the highest-AUC studies;\n")
cat(" the logit scale is therefore the primary analysis.)\n")
```

## Forest plot

```{r forest, fig.width=10, fig.height=9, fig.cap="Forest plot of the random-effects (DerSimonian-Laird) meta-analysis of model discrimination for the 12 poolable Family A (ENM/SDM) studies. Effects are pooled on the logit(AUC) scale; the diamond is the pooled estimate and the red bar the 95\\% prediction interval. The column headed logit/(AUC) gives each study's effect on the logit scale and SE its standard error."}
forest(m_auc,
       print.tau2 = TRUE,
       print.I2   = TRUE,
       prediction = TRUE,
       leftlabs   = c("Study", "logit\n(AUC)", "SE"),
       xlab       = "logit(AUC)",
       smlab      = "Pooled AUC\n(logit, DL)",
       col.diamond = "steelblue",
       col.predict = "darkred",
       colgap.left = "5mm",
       fontsize   = 10)

# AUC on original scale (back-transformed) for readability
cat("\n--- AUC on original scale (back-transformed) ---\n")
for (i in 1:nrow(dat_auc)) {
  ci_lo <- plogis(dat_auc$logit_auc[i] - 1.96 * dat_auc$se_logit[i])
  ci_hi <- plogis(dat_auc$logit_auc[i] + 1.96 * dat_auc$se_logit[i])
  cat(sprintf("  %-30s AUC = %.3f [%.3f - %.3f]\n",
              dat_auc$label[i], dat_auc$auc[i], ci_lo, ci_hi))
}
cat(sprintf("\n  POOLED                         AUC = %.3f [%.3f - %.3f]\n",
            pooled_auc, pooled_lower, pooled_upper))
```

The random-effects meta-analysis of 12 Family A studies yielded a pooled AUC of 0.892 (95% CI: 0.859--0.917), indicating that ecological niche and species distribution models for *B. anthracis* in Africa achieve, on average, good-to-excellent discriminatory capacity (AUC values in the 0.8--0.9 range are conventionally classed as good and those above 0.9 as excellent). Pooling was performed on the logit scale under a DerSimonian & Laird random-effects model [@dersimonian1986] and back-transformed for presentation, which keeps every confidence limit within the admissible (0,1) range. The use of AUC as the primary discrimination metric, and the cautions attached to interpreting it across heterogeneous ecological niche and species distribution models, follow standard guidance for this study type [@araujo2006]. A raw-scale sensitivity analysis gave a slightly higher pooled AUC of 0.903 (95% CI 0.881--0.926) but produced upper limits exceeding 1 for the highest-AUC studies, confirming the logit scale as the appropriate primary analysis.

Substantial heterogeneity was detected across studies ($I^2$ = 91%, Cochran's *Q* test *p* < 0.0001), indicating that the observed variation in AUC values cannot be attributed to sampling error alone and reflects genuine differences in study design, geographic context, modelling algorithm, sample size, and covariate selection. The wide prediction interval spans from acceptable to near-perfect discrimination, meaning that a new study conducted under broadly similar conditions could plausibly report an AUC anywhere within this range. This level of heterogeneity is expected given the diversity of algorithms (MaxEnt, BRT, ensemble, Random Forest, GARP), spatial scales (park-level to global), and occurrence sample sizes (n = 44 to n = 5,018) represented in the corpus.

The forest plot shows that the highest-performing study was a biomod2 ensemble [@assefa2020] (AUC = 0.991), though this extreme value warrants caution as it may reflect overfitting in a relatively small sample (n = 128). At the other end, a MaxEnt model [@yousuf2020] (AUC = 0.750) and a BRT [@otieno2021a] (AUC = 0.800) exhibited the lowest discrimination, still within the range classed as acceptable to good. Studies employing large, multi-country datasets with rigorous spatial thinning [@carlson2019; @otieno2021b] achieved AUC values close to 0.93 with very narrow confidence intervals, reflecting the precision gains afforded by larger sample sizes.

## Subgroup analyses

### By ENM method type

```{r subgroup-method, fig.width=10, fig.height=10, fig.cap="Subgroup forest plot of pooled discrimination for the Family A studies stratified by ENM/SDM algorithm (x-axis = logit(AUC); DerSimonian-Laird random effects). Only algorithm families with at least two studies (MaxEnt, K = 6; boosted regression trees, K = 3) are pooled separately; single-study categories are shown but not pooled. The chi-square tests for subgroup differences are printed at the foot of the plot."}
if (length(unique(dat_auc$method_type)) > 1) {
  m_sub_method <- metagen(
    TE = logit_auc, seTE = se_logit, studlab = label,
    data = dat_auc, sm = "MRAW", method.tau = "DL",
    subgroup = method_type, tau.common = FALSE
  )
  cat("=== SUBGROUPS BY METHOD ===\n")
  print(summary(m_sub_method))

  # Pre-specified formal comparison: only MaxEnt and BRT have >= 2 studies.
  # The omnibus 5-category test is inflated by single-study categories
  # (Ensemble AUC 0.991, GARP AUC 0.710) and is NOT a valid between-algorithm test.
  d2 <- subset(dat_auc, method_type %in% c("MaxEnt", "BRT"))
  m_2grp <- metagen(TE = logit_auc, seTE = se_logit, studlab = label,
                    data = d2, sm = "MRAW", method.tau = "DL",
                    subgroup = method_type, tau.common = FALSE)
  cat(sprintf("\nFormal MaxEnt vs BRT test (2 groups): Q = %.2f, df = %d, p = %.4f\n",
              m_2grp$Q.b.random, m_2grp$df.Q.b, m_2grp$pval.Q.b.random))
  cat(sprintf("Omnibus 5-category test (incl. single-study Ensemble/RF/GARP): Q = %.2f, df = %d, p = %.4f\n",
              m_sub_method$Q.b.random, m_sub_method$df.Q.b, m_sub_method$pval.Q.b.random))

  forest(m_sub_method,
         print.subgroup.name = TRUE,
         col.diamond = "steelblue",
         xlab = "",          # omitted: the axis label collided with the subgroup-difference test line; scale stated in the caption
         fontsize = 9)
} else {
  cat("Subgroup by method: insufficient categories.\n")
}
```

A pre-specified subgroup analysis was conducted for the two algorithm families with at least two studies. MaxEnt models (*K* = 6) yielded a pooled AUC of 0.895 (95% CI 0.826--0.939) and boosted regression trees (*K* = 3) a pooled AUC of 0.914 (95% CI 0.890--0.934); the Ensemble, Random Forest, and GARP categories each contained a single study and were not pooled separately. The pre-specified formal comparison between the two multi-study algorithm families (MaxEnt vs boosted regression trees) was non-significant (test for subgroup differences $Q$ = 0.45, df = 1, *p* = 0.50); the omnibus test across all five method categories is significant (*p* < 0.0001) only because of the single-study Ensemble (AUC 0.991) and GARP (AUC 0.710) categories, which are statistically unstable and provide no valid basis for between-algorithm inference, so only the two-group comparison is reported in the manuscript. The two pooled subgroups overlap substantially, suggesting that the choice of algorithm per se does not systematically determine discriminatory performance in African anthrax ENM, and that other factors -- occurrence data quality, spatial extent, environmental predictor selection, and validation strategy -- exert a stronger influence on AUC. Heterogeneity remained high within the MaxEnt subgroup, reinforcing the interpretation that within-algorithm variation exceeds between-algorithm variation. The practical implication is that methodological quality and the ecological relevance of predictor variables are likely more important determinants of model accuracy than the specific machine-learning algorithm employed.

### By geographic region

```{r subgroup-region}
if (length(unique(dat_auc$region)) > 1) {
  m_sub_region <- metagen(
    TE = logit_auc, seTE = se_logit, studlab = label,
    data = dat_auc, sm = "MRAW", method.tau = "DL",
    subgroup = region, tau.common = FALSE
  )
  cat("=== SUBGROUPS BY REGION ===\n")
  print(summary(m_sub_region))
} else {
  cat("Subgroup by region: insufficient categories.\n")
}
```

The region subgroups must be interpreted with caution: East Africa contributes the majority of poolable studies, while the Southern, Pan-African, and (in the AUC pool) West African strata each contain only one or two studies. As reported in the meta-regression below, region did not explain the between-study heterogeneity once examined formally (all coefficients non-significant), and the single-study strata are statistically unstable. Any apparent regional differences may reflect genuine ecological variation in the predictability of anthrax suitability across biomes, but they may equally reflect differences in data quality, surveillance intensity, and modelling methodology between regions.

## Meta-regression

Sources of heterogeneity were explored with mixed-effects (moment-based) meta-regression on the logit-transformed AUC, fitting each moderator separately, with residual $I^2$ recomputed for each model. These analyses are exploratory given the small number of studies (*K* = 12; *K* = 11 for sample size, as one study did not report *n*).

### AUC vs. publication year

```{r metareg-year, fig.cap="Bubble plot of the meta-regression of logit(AUC) on publication year for the 12 Family A studies. Each bubble is one study, sized by its inverse-variance weight; the line is the weighted mixed-effects fit with its 95\\% confidence band."}
dat_rma <- escalc(measure = "GEN", yi = logit_auc, sei = se_logit, data = dat_auc)
m_reg_year <- rma(yi, vi, mods = ~ year, data = dat_rma, method = "DL")
cat("=== META-REGRESSION: Publication year ===\n")
print(m_reg_year)

bubble_dat <- data.frame(year = dat_auc$year, logit = dat_auc$logit_auc,
                         w = 1 / dat_auc$se_logit^2)
ggplot(bubble_dat, aes(year, logit, size = w)) +
  geom_point(alpha = 0.5, colour = "darkblue") +
  geom_smooth(method = "lm", se = TRUE, colour = "darkblue",
              aes(weight = w), linewidth = 1) +
  labs(x = "Publication year", y = "logit(AUC)", size = "Weight") +
  theme_minimal(base_size = 12)
```

Publication year showed a borderline positive association with discrimination (logit-AUC slope $\approx$ +0.131 per year, 95% CI -0.004 to +0.266, *p* = 0.057): more recent models tended to report higher AUCs, but the effect did not reach significance and residual heterogeneity remained high (residual $I^2 \approx$ 90%). This near-significant trend could reflect methodological maturation (better regularisation, ensemble techniques, higher-resolution covariates), but equally a shift toward easier-to-model settings or simple reporting drift; with only 12 studies it should be treated as hypothesis-generating rather than conclusive.

### AUC vs. sample size

```{r metareg-n, fig.cap="Bubble plot of the meta-regression of logit(AUC) on log10 of the number of presence records (K = 11; one study did not report n). Each bubble is one study, sized by its inverse-variance weight; the line is the weighted mixed-effects fit with its 95\\% confidence band."}
dat_auc_n <- dat_auc %>% filter(!is.na(n_obs)) %>%
  mutate(log10_n = log10(n_obs))
dat_rma_n <- escalc(measure = "GEN", yi = logit_auc, sei = se_logit, data = dat_auc_n)
m_reg_n <- rma(yi, vi, mods = ~ log10_n, data = dat_rma_n, method = "DL")
cat("=== META-REGRESSION: log10(N presence records) ===\n")
print(m_reg_n)

bub_n <- data.frame(log10_n = dat_auc_n$log10_n, logit = dat_auc_n$logit_auc,
                    w = 1 / dat_auc_n$se_logit^2)
ggplot(bub_n, aes(log10_n, logit, size = w)) +
  geom_point(alpha = 0.5, colour = "darkred") +
  geom_smooth(method = "lm", se = TRUE, colour = "darkred",
              aes(weight = w), linewidth = 1) +
  labs(x = "log10(presence records)", y = "logit(AUC)", size = "Weight") +
  theme_minimal(base_size = 12)
```

Sample size was not associated with AUC (slope per log10 *n* $\approx$ +0.220, 95% CI -0.366 to +0.805, *p* = 0.462; residual $I^2 \approx$ 91%). This counters the intuitive expectation that more data should improve fit, but the relationship between sample size and AUC in ecological niche modelling is complex and non-monotonic: very small datasets may achieve high AUC through overfitting, while large multi-context datasets may average over diverse ecological settings and so reduce apparent discrimination. The inclusion of both park-level studies (n = 44--219) and continental/global studies (n = 874--5,018) likely flattens any linear trend. The pattern reinforces the importance of evaluating performance not only through AUC but also through omission rates, the true skill statistic (TSS), and spatial transferability.

### Region

Region was also examined as a moderator (East as reference). None of the regional contrasts was significant (West $\approx$ -0.42, Southern $\approx$ -0.48, Pan-African $\approx$ +0.23 on the logit scale; all *p* > 0.4) and residual heterogeneity remained $\approx$ 92%. The Southern and Pan-African strata each contained a single study and are therefore unstable. Overall, none of the examined moderators accounted for the substantial heterogeneity, which remains largely unexplained.

# Publication bias assessment

Small-study effects and publication bias were assessed visually with funnel plots (standard and contour-enhanced), formally with Egger's regression test, and with the trim-and-fill procedure.

```{r pub-bias, fig.height=7, fig.cap="Funnel plots for the AUC meta-analysis. Each point is one of the 12 poolable Family A (ENM/SDM) studies, plotted by its effect (logit(AUC), x-axis) against its standard error (y-axis, inverted); the dashed vertical line is the pooled estimate and the diagonal lines the pseudo-95\\% confidence limits. Top: standard funnel. Bottom: contour-enhanced funnel, with shaded bands marking conventional significance regions (the studies are individually labelled in the forest plot above)."}
par(mar = c(5, 4.5, 1.5, 2))   # extra bottom/left margin so axis labels are not clipped
# studlab omitted: point labels overlap/clip on a funnel; each point = one Family A study (see forest plot)
meta::funnel(m_auc, xlab = "logit(AUC)",
             col = "navy", bg = "steelblue")
legend("topright", legend = "Each point = one Family A study",
       pch = 21, pt.bg = "steelblue", col = "navy", bg = "white", bty = "o")

par(mar = c(5, 4.5, 1.5, 2))
meta::funnel(m_auc, xlab = "logit(AUC)",
             contour = c(0.9, 0.95, 0.99),
             col.contour = c("gray75", "gray85", "gray95"))
legend("topright",
       legend = c("Each point = one Family A study",
                  "p < 0.1", "p < 0.05", "p < 0.01"),
       pch = c(21, NA, NA, NA), pt.bg = c("black", NA, NA, NA),
       fill = c(NA, "gray75", "gray85", "gray95"),
       border = NA, bg = "white")
```

```{r pub-bias-tests, fig.cap="Trim-and-fill funnel plot for the AUC meta-analysis. Filled points are the 12 observed Family A studies and any open points are studies imputed by the trim-and-fill algorithm to restore funnel symmetry; the x-axis is logit(AUC) and the y-axis the standard error."}
cat("==================================================\n")
cat("PUBLICATION BIAS TESTS\n")
cat("==================================================\n\n")

egger <- metabias(m_auc, method.bias = "linreg", k.min = 5)
cat("Egger test (linear regression):\n")
cat(sprintf("  intercept = %.3f, t = %.3f, p = %.4f\n",
            egger$estimate[1], egger$statistic, egger$p.value))
cat(sprintf("  => %s\n\n",
    ifelse(egger$p.value < 0.05,
           "SIGNIFICANT: possible small-study / publication bias",
           "Not significant: no strong evidence of bias")))

tf <- trimfill(m_auc)
cat("=== TRIM AND FILL ===\n")
cat(sprintf("Imputed studies: %d\n", tf$k0))
print(summary(tf))

par(mar = c(5, 4.5, 1.5, 2))
meta::funnel(tf, xlab = "logit(AUC)")
legend("topright",
       legend = c("Observed study", "Imputed (trim-and-fill)"),
       pch = c(16, 1), col = "black", bg = "white", bty = "o")
```

Egger's regression test for funnel-plot asymmetry was not statistically significant (intercept $\approx$ -2.15, $t$ = -1.45, df = 10, *p* = 0.179), providing no formal evidence of small-study effects; the funnel was visually roughly symmetric, and the contour-enhanced plot showed most studies in regions of statistical significance, as is typical when AUC values sit well above 0.5. The trim-and-fill analysis imputed few or no missing studies, indicating that the observed distribution does not suggest suppression of low-performing models.

Taken together, these results indicate that the pooled AUC estimate of 0.892 is unlikely to be materially inflated by publication bias. With only 12 studies, however, both the meta-regressions and Egger's test have low statistical power, so non-significant findings do not exclude moderate moderator effects or small-study bias. The concept of publication bias in ecological modelling also differs from clinical epidemiology: authors of ENM/SDM studies may simply not report a discrimination metric (rather than selectively reporting high ones), and seven Family A studies were excluded from pooling precisely because they reported pROC ratios, sensitivity, or no comparable AUC. This reporting heterogeneity may introduce a selection effect that standard publication-bias tests do not capture.

# Partial quantitative synthesis: risk-factor odds ratios (Family C)

Nine Family C studies (analytical observational designs: outbreak case-control, matched case-control, and cross-sectional surveys) reported adjusted odds ratios (aOR) for human anthrax. Two were not poolable: a herd/animal-level ecological analysis [@kamboyi2019] and a study reporting only a source-specific meat-origin odds ratio [@nakanwagi2020]. Exposures reported by at least three studies on a comparable individual, disease-outcome basis were pooled on the natural-log scale (SE from the reported 95% CI) under a DerSimonian-Laird random-effects model and back-transformed; heterogeneity was quantified with $I^2$.

```{r family-c-data}
# --- Exposure 1: consumption of meat from dead/sick animals ---
exp1 <- data.frame(
  study = c("Mutesi et al. (2024)", "Tumusiime et al. (2025)",
            "Hamutyinei Dhliwayo et al. (2022)", "Migisha et al. (2021)"),
  aor   = c(5.80, 6.19, 12.20, 6.40),
  lo    = c(2.70, 2.76, 1.41, 1.53),
  hi    = c(12.00, 13.90, 105.74, 26.70)
)
# --- Exposure 2: slaughtering / skinning / handling carcasses ---
exp2 <- data.frame(
  study = c("Mutesi et al. (2024)", "Kwizera et al. (2025)",
            "Makurumidze et al. (2021)", "Migisha et al. (2021)",
            "Lehman et al. (2017)"),
  aor   = c(9.30, 6.00, 4.32, 5.74, 12.00),
  lo    = c(2.50, 1.80, 1.25, 1.39, 4.30),
  hi    = c(15.00, 19.00, 14.94, 23.80, 36.50)
)

pool_or <- function(d, label) {
  logOR <- log(d$aor)
  seLog <- (log(d$hi) - log(d$lo)) / (2 * 1.96)
  m <- metagen(TE = logOR, seTE = seLog, studlab = d$study,
               sm = "OR", method.tau = "DL", prediction = FALSE)
  cat("===", label, "===\n")
  cat(sprintf("  Pooled aOR = %.2f (95%% CI %.2f - %.2f); I^2 = %.0f%%; k = %d\n\n",
              exp(m$TE.random), exp(m$lower.random), exp(m$upper.random),
              m$I2 * 100, m$k))
  m
}

m_exp1 <- pool_or(exp1, "Exposure 1: meat from dead/sick animals")
m_exp2 <- pool_or(exp2, "Exposure 2: slaughtering/skinning/handling carcasses")
```

## Forest plots — Family C

```{r family-c-forest1, fig.height=4, fig.cap="Forest plot of adjusted odds ratios (random-effects, DerSimonian-Laird) for consumption of meat from dead or sick animals as a risk factor for human anthrax (Family C, K = 4). Effects are pooled on the log-OR scale and shown back-transformed; the diamond is the pooled aOR."}
forest(m_exp1, xlab = "Adjusted OR (log scale)",
       smlab = "Meat from dead/sick animals",
       col.diamond = "firebrick", leftlabs = c("Study", "log(OR)", "SE"))
```

```{r family-c-forest2, fig.height=4.5, fig.cap="Forest plot of adjusted odds ratios (random-effects, DerSimonian-Laird) for slaughtering, skinning or handling carcasses as a risk factor for human anthrax (Family C, K = 5). Effects are pooled on the log-OR scale and shown back-transformed; the diamond is the pooled aOR."}
forest(m_exp2, xlab = "Adjusted OR (log scale)",
       smlab = "Slaughtering / handling carcasses",
       col.diamond = "firebrick", leftlabs = c("Study", "log(OR)", "SE"))
```

Consumption of meat from dead or sick animals was strongly associated with human anthrax (pooled aOR = 6.26, 95% CI 3.80--10.30; *K* = 4), with no detectable heterogeneity ($I^2$ = 0%). Slaughtering, skinning, or handling carcasses showed a comparable and consistent association (pooled aOR = 7.55, 95% CI 4.58--12.46; *K* = 5; $I^2$ = 0%). Both exposures therefore increase the odds of infection roughly six- to eight-fold, and the effect was highly consistent across countries (Uganda, Zimbabwe, Zambia) and years, pointing to two dominant, modifiable transmission pathways for targeted prevention (safe carcass disposal, meat inspection, and community education). The very wide interval for one study [@hamutyinei2022] reflects a small sample, and the absence of heterogeneity, while reassuring, rests on few studies (*K* = 4--5).

# Complementary syntheses: Families B, D and E (narrative)

Families B and D do not support a discrimination meta-analysis because their reported metrics are non-comparable (Bayesian INLA models validated by sensitivity; hotspot statistics yielding location-specific cluster outputs rather than a common effect size). The three Family B studies (a Kenyan INLA model [@ndolo2022kenya], sensitivity 75%; a Ugandan INLA model [@ndolo2022uganda], sensitivity 90%; and a Bayesian spatial GLMM [@aminu2022], AUC 0.94) consistently identified proximity to water and soil chemistry (calcium, organic carbon) as drivers and revealed previously unmapped hotspots. The seven Family D studies (Getis-Ord Gi*, KDE/LISA, SaTScan, kriging, spatial clustering) consistently found anthrax to be strongly spatially clustered -- repeatedly concentrating in a small fraction of the territory -- and linked to water bodies, alkaline soils, and livestock density. The three Family E studies project suitability under climate scenarios; their discrimination AUCs are already represented in the Family A analysis, and a stand-alone pool of the two reporting an exact AUC (AUC 0.929 [@otieno2021b]; AUC 0.750 [@yousuf2020]) is uninformative ($I^2 \approx$ 98%) and shown only for completeness. The substantive Family E output is the projected change in suitable area, which is scenario- and region-dependent (general northward/eastward shift and expansion under intermediate pathways in East Africa, contraction in parts of Ethiopia and southern Africa under high-emission pathways). Full extraction tables for these families are provided in `Families_BDE_summary.docx`.

# Sensitivity analysis: exclusion of preprints

Two of the 33 included studies are preprints that have not undergone formal peer review: R8 [@mutesi2025] (Research Square) and R47 [@yousuf2020] (bioRxiv). Preprints were included because they had accessible full texts with verifiable methods, and excluding them would risk a selection effect that systematically removes studies reporting less favourable results or originating from settings with less access to high-impact journals. R47 (AUC = 0.750) is the lowest-performing study in the AUC pool, so its exclusion would inflate the pooled estimate.

To assess robustness, the random-effects AUC meta-analysis was repeated excluding R47 (R8 is Family C and does not contribute to the AUC pool).

```{r sensitivity-preprint}
cat("===================================================\n")
cat("SENSITIVITY: Excluding preprints from AUC pooling\n")
cat("===================================================\n\n")

dat_auc_np <- dat_auc %>% filter(rec_id != 47)
cat(sprintf("Studies after excluding R47: K = %d\n", nrow(dat_auc_np)))
cat("Excluded: R47 Yousuf et al. (2020) -- bioRxiv preprint, AUC = 0.750\n\n")

m_np <- metagen(TE = logit_auc, seTE = se_logit, studlab = label,
                data = dat_auc_np, sm = "MRAW", method.tau = "DL",
                prediction = TRUE)

cat(sprintf("Pooled AUC (no preprints): %.3f (95%% CI %.3f - %.3f), I^2 = %.0f%%\n\n",
            plogis(m_np$TE.random), plogis(m_np$lower.random),
            plogis(m_np$upper.random), m_np$I2 * 100))

cat("--- Comparison ---\n")
cat(sprintf("Main analysis (K=%d): AUC = %.3f [%.3f - %.3f], I^2 = %.0f%%\n",
            m_auc$k, pooled_auc, pooled_lower, pooled_upper, m_auc$I2 * 100))
cat(sprintf("No preprints  (K=%d): AUC = %.3f [%.3f - %.3f], I^2 = %.0f%%\n",
            m_np$k, plogis(m_np$TE.random), plogis(m_np$lower.random),
            plogis(m_np$upper.random), m_np$I2 * 100))
```

As expected, excluding the lowest-performing study raises the pooled AUC slightly while heterogeneity remains substantial. The direction and magnitude of the change confirm that including preprints yields a more conservative, defensible estimate, and the main conclusions -- good-to-excellent average discrimination with high context-dependent heterogeneity -- are robust to the inclusion or exclusion of preprints. A peer-reviewed-only sensitivity analysis is pre-specified at the pooling stage of the full review.

# Sensitivity analysis: exclusion of high-risk-of-bias studies

Three of the 12 pooled Family A studies were rated at high overall risk of bias on the custom seven-domain tool (R51, Carlson et al. 2019; R52, Steenkamp et al. 2018; R63, Blackburn et al. 2015); the remaining nine were rated moderate and none low. To verify that the pooled discrimination is not driven by lower-quality studies, the random-effects AUC meta-analysis was repeated excluding these three.

```{r sensitivity-rob}
cat("===================================================\n")
cat("SENSITIVITY: Excluding high-risk-of-bias studies\n")
cat("===================================================\n\n")

dat_auc_lowrob <- dat_auc %>% filter(!(rec_id %in% c(51, 52, 63)))
cat(sprintf("Studies after excluding high-RoB (R51, R52): K = %d\n", nrow(dat_auc_lowrob)))
cat("Excluded: R51 Carlson et al. (2019); R52 Steenkamp et al. (2018); R63 Blackburn et al. (2015)\n\n")

m_rob <- metagen(TE = logit_auc, seTE = se_logit, studlab = label,
                 data = dat_auc_lowrob, sm = "MRAW", method.tau = "DL",
                 prediction = TRUE)

cat(sprintf("Pooled AUC (excl. high-RoB): %.3f (95%% CI %.3f - %.3f), I^2 = %.0f%%\n\n",
            plogis(m_rob$TE.random), plogis(m_rob$lower.random),
            plogis(m_rob$upper.random), m_rob$I2 * 100))

cat("--- Comparison ---\n")
cat(sprintf("Main analysis (K=%d):  AUC = %.3f [%.3f - %.3f], I^2 = %.0f%%\n",
            m_auc$k, pooled_auc, pooled_lower, pooled_upper, m_auc$I2 * 100))
cat(sprintf("Excl. high-RoB (K=%d): AUC = %.3f [%.3f - %.3f], I^2 = %.0f%%\n",
            m_rob$k, plogis(m_rob$TE.random), plogis(m_rob$lower.random),
            plogis(m_rob$upper.random), m_rob$I2 * 100))
```

Excluding the three high-risk-of-bias studies left the conclusion unchanged, with the pooled estimate marginally higher (AUC = 0.903, 95% CI 0.859--0.934; *K* = 9) and heterogeneity still high ($I^2 \approx$ 92%). The good-to-excellent average discrimination is therefore robust to the exclusion of lower-quality studies; the slight increase is expected because the excluded set includes the lowest-AUC model (R63, AUC 0.71).

# Sensitivity analysis: leave-one-out and fixed-effect model

To check that no single study drives the pooled estimate, a leave-one-out analysis was performed (each study omitted in turn and the pool recomputed), and a fixed-effect (common-effect) model was fitted for comparison with the primary random-effects analysis.

```{r sensitivity-loo}
# Leave-one-out (random-effects)
loo <- metainf(m_auc, pooled = "random")
loo_auc <- plogis(loo$TE)
loo_auc <- loo_auc[!is.na(loo_auc) & loo$studlab != ""]
cat("===================================================\n")
cat("SENSITIVITY: Leave-one-out and fixed-effect model\n")
cat("===================================================\n\n")
cat(sprintf("Leave-one-out pooled AUC range: %.3f - %.3f (main %.3f)\n",
            min(loo_auc), max(loo_auc), pooled_auc))
cat("  => no single study materially changes the pooled estimate.\n\n")

cat(sprintf("Fixed-effect (common) AUC: %.3f (95%% CI %.3f - %.3f)\n",
            plogis(m_auc$TE.common), plogis(m_auc$lower.common),
            plogis(m_auc$upper.common)))
cat(sprintf("Random-effects     AUC: %.3f (95%% CI %.3f - %.3f)\n",
            pooled_auc, pooled_lower, pooled_upper))
```

The leave-one-out estimates spanned only AUC 0.886--0.902 around the overall value of 0.892, confirming that no individual study -- including the high-weight, large multi-country datasets -- materially influenced the pooled discrimination. A fixed-effect model gave a consistent, slightly higher estimate (AUC = 0.915, 95% CI 0.909--0.920); the higher value is expected because the common-effect model up-weights the few very precise studies, and the substantial heterogeneity makes the random-effects estimate the appropriate primary analysis.

# Integrated synthesis and the persistence-activation framework

The quantitative synthesis demonstrates that ecological niche and species distribution models for anthrax in Africa achieve, on average, good-to-excellent discriminatory performance (pooled AUC = 0.892), while the Family C analysis identifies two dominant, highly consistent behavioural transmission pathways for human infection (eating meat from dead/sick animals, and slaughtering or handling carcasses). Discriminatory performance alone, however, does not guarantee ecological or epidemiological insight. A model that accurately classifies high- versus low-risk areas may do so by leveraging persistent environmental gradients (soil properties, topography, geology) without capturing the dynamic conditions that trigger actual outbreaks (rainfall anomalies, flooding, seasonal vegetation pulses).

This distinction lies at the heart of the persistence-activation framework that structures this doctoral thesis. The persistence layer comprises relatively stable environmental variables -- soil pH, calcium concentrations, clay content, topography, and land cover -- that define where *Bacillus anthracis* spores can remain viable over decades. The activation layer includes dynamic climatic and hydrological variables -- rainfall, flooding, NDVI, soil moisture, and temperature -- that determine when conditions become conducive to spore re-exposure and outbreak emergence. The reviewed ENM/SDM studies predominantly capture the persistence layer, because the bioclimatic and soil variables typically used as MaxEnt or BRT covariates reflect long-term averages rather than event-specific triggers. This is both a strength (persistence-driven suitability is stable and mappable) and a limitation (it cannot predict when an outbreak will occur in a suitable area).

The high heterogeneity observed in this meta-analysis ($I^2$ = 91%) may partly reflect differences in the degree to which studies capture persistence versus activation dynamics. Studies that incorporate seasonal or inter-annual variables (e.g., one study used precipitation seasonality and monthly NDVI [@pittiglio2022]; another used precipitation of the driest month [@kracalik2017]) may implicitly incorporate activation signals, whereas purely bioclimatic MaxEnt models rely almost entirely on persistence-related spatial gradients. Future research should explicitly decompose model performance into its persistence and activation components, ideally through Bayesian spatiotemporal frameworks that can accommodate both layers within a single inferential architecture.

# Validation of the Systematic Review Toolkit: independent reproduction of inter-rater agreement

Screening and risk-of-bias assessment were captured and consolidated with a
purpose-built, browser-based **Systematic Review Toolkit**. Because that
interface is a bespoke data-capture tool rather than formally validated software,
every inter-rater agreement statistic it reports is re-derived here, from the
raw reviewer-level exports (`MMDN_Screening.json`, `Reviewer_Screening.json`,
`MMDN_RoB.json`, `Experts_RoB.json`), using the established **irr** and **psych**
packages. The Toolkit therefore contributes only data capture and consolidation;
all reported agreement statistics are reproducible with standard, peer-reviewed
software. The screening coefficient is the Cohen's $\kappa$ for the
title/abstract validation sample (lead reviewer MMDN vs the independent reviewer
MJRL); the risk-of-bias coefficient is the quadratic-weighted $\kappa$ between the
lead reviewer and the domain experts across all dual-rated study $\times$ domain
ratings, before consensus reconciliation.

```{r toolkit-validation, echo=TRUE}
## ---- Screening: validation sample (MMDN vs independent reviewer MJRL) ----
## Source: SR_Toolkit exports MMDN_Screening.json & Reviewer_Screening.json
screen <- data.frame(
  rec  = c(8, 13, 24, 27, 28, 42, 47, 51, 52, 56, 59, 60, 61),
  MMDN = c("Include","Exclude","Exclude","Uncertain","Include","Include",
           "Include","Include","Include","Exclude","Exclude","Exclude","Include"),
  MJRL = c("Include","Exclude","Exclude","Exclude","Include","Include",
           "Include","Include","Include","Exclude","Exclude","Exclude","Include"),
  stringsAsFactors = FALSE
)
# Title/abstract screening is a binary keep/drop decision: records flagged
# 'Uncertain' are carried forward to full text (i.e. treated as 'Include').
screen$MMDN_bin <- ifelse(screen$MMDN == "Exclude", "Exclude", "Include")
screen$MJRL_bin <- ifelse(screen$MJRL == "Exclude", "Exclude", "Include")
# Numeric coding for psych (Exclude = 1, Include = 2)
screen$MMDN_num <- ifelse(screen$MMDN_bin == "Include", 2L, 1L)
screen$MJRL_num <- ifelse(screen$MJRL_bin == "Include", 2L, 1L)

k_screen   <- irr::kappa2(screen[, c("MMDN_bin", "MJRL_bin")])
ck_screen  <- psych::cohen.kappa(cbind(screen$MMDN_num, screen$MJRL_num))
raw_screen <- mean(screen$MMDN_bin == screen$MJRL_bin)

## ---- Risk of bias: MMDN vs domain-expert ratings (ordinal Low/Moderate/High) ----
## Source: SR_Toolkit exports MMDN_RoB.json & Experts_RoB.json
## 29 dual-rated studies x 7 domains = 203 paired ordinal ratings.
## Confusion matrix (rows = MMDN, cols = experts):
rob_cm <- matrix(c(41, 24,  6,
                   35, 54,  8,
                    3, 29,  3),
                 nrow = 3, byrow = TRUE,
                 dimnames = list(MMDN    = c("Low","Moderate","High"),
                                 Experts = c("Low","Moderate","High")))
# Ordinal numeric coding (Low = 1, Moderate = 2, High = 3). Numeric coding is
# required because irr::kappa2 sorts factor *labels* alphabetically, which would
# scramble the ordinal order (High < Low < Moderate) and invalidate the weighted
# distances; numeric levels preserve the Low < Moderate < High ordering.
mm <- integer(0); ex <- integer(0)
for (i in 1:3) for (j in 1:3) if (rob_cm[i, j] > 0) {
  mm <- c(mm, rep(i, rob_cm[i, j])); ex <- c(ex, rep(j, rob_cm[i, j]))
}
rob_pairs <- data.frame(MMDN = mm, Experts = ex)
k_rob_unw <- irr::kappa2(rob_pairs, weight = "unweighted")
k_rob_lin <- irr::kappa2(rob_pairs, weight = "equal")     # linear weights
k_rob_sq  <- irr::kappa2(rob_pairs, weight = "squared")   # quadratic weights (primary)
raw_rob   <- mean(rob_pairs$MMDN == rob_pairs$Experts)

cat("=====================================================================\n")
cat("TOOLKIT VALIDATION -- inter-rater agreement reproduced in R\n")
cat("=====================================================================\n\n")
cat(sprintf("SCREENING (binary keep/drop; validation sample, n = %d)\n", nrow(screen)))
cat("  Toolkit-reported Cohen's kappa : 0.843\n")
cat(sprintf("  R (irr::kappa2)                : %.3f\n", k_screen$value))
cat(sprintf("  R (psych::cohen.kappa, 95%% CI) : %.3f (%.2f - %.2f)\n",
            ck_screen$kappa, ck_screen$confid[1, 1], ck_screen$confid[1, 3]))
cat(sprintf("  Raw agreement                  : %.1f%% (%d/%d)\n\n",
            100 * raw_screen, sum(screen$MMDN_bin == screen$MJRL_bin), nrow(screen)))
cat(sprintf("RISK OF BIAS (ordinal Low/Moderate/High; %d study x domain ratings)\n",
            nrow(rob_pairs)))
cat("  Toolkit-reported weighted kappa : 0.267 (quadratic), 48.3% raw\n")
cat(sprintf("  R (irr, squared/quadratic)      : %.3f\n", k_rob_sq$value))
cat(sprintf("  R (irr, equal/linear)           : %.3f\n", k_rob_lin$value))
cat(sprintf("  R (irr, unweighted)             : %.3f\n", k_rob_unw$value))
cat(sprintf("  Raw agreement                   : %.1f%%\n", 100 * raw_rob))
```

The independent recomputation reproduced the Toolkit's statistics exactly. The
screening Cohen's $\kappa$ was 0.843 (95% CI 0.55--1.00; 12 of 13 records
concordant, 92.3% raw agreement), and the risk-of-bias quadratic-weighted
$\kappa$ was 0.267 (48.3% raw agreement across 203 study $\times$ domain ratings).
Because these values are identical to those produced by the browser-based
interface, the Toolkit's agreement computations are confirmed to be correct and
fully reproducible in standard, validated R packages; the interface itself
performs only data capture and consolidation, and introduces no bespoke
statistics.

The contrast between the two coefficients is itself informative and is reported
transparently rather than concealed. Title/abstract screening, a binary
decision, achieved almost perfect agreement ($\kappa$ = 0.843), supporting the
reliability of study selection. Risk-of-bias rating, by contrast, showed only
fair chance-corrected agreement before reconciliation (weighted $\kappa$ = 0.267)
despite moderate raw agreement (48.3%) -- a well-recognised pattern (the
"kappa paradox") that arises when an ordinal instrument is applied to a
methodologically diverse literature and most ratings concentrate in the central
"Moderate" category, which deflates chance-corrected coefficients. This low
pre-consensus agreement is precisely why **every domain of every study was
independently dual-rated and then reconciled to a documented consensus**, rather
than relying on any single rater or on a sampled check; the final risk-of-bias
ratings therefore do not depend on inter-rater agreement at the rating stage.
Reporting both coefficients, together with their exact reproduction in R, makes
the reliability of the whole process auditable end to end.

# Consistency with the manuscript

```{r prisma-counts, results='asis'}
cat("## PRISMA counts\n\n")
cat("| PRISMA node | Count |\n")
cat("|-------------|-------|\n")
cat("| PubMed records | 28 |\n")
cat("| Scopus records | 54 |\n")
cat("| Web of Science records | 40 |\n")
cat("| Scopus Preprints records | 8 |\n")
cat("| Total hits across databases | 130 |\n")
cat("| Duplicates removed | 67 |\n")
cat("| Unique records screened | 63 |\n")
cat("| Excluded at title/abstract | 28 |\n")
cat("| Assessed for eligibility (full text) | 35 |\n")
cat("| Excluded at full text (R27, R53) | 2 |\n")
cat("| Studies included | 33 |\n")
cat("\n")
cat("Counts reconciled with `S1_Screening_supplementary.docx` (PRISMA flow).\n")
cat("Note: 30 records were excluded overall (28 at title/abstract + 2 at full text).\n")
```

The systematic review was conducted and reported in accordance with the PRISMA 2020 statement [@page2021]. The systematic search was conducted in **three bibliographic databases** -- PubMed, Scopus, and Web of Science -- complemented by the Scopus preprint server; Embase and Google Scholar were not used as primary sources. The PRISMA counts above are consistent with these three databases plus the preprint server (28 + 54 + 40 + 8 = 130 hits). The master protocol and manuscript should therefore list these three databases (plus the preprint server) when the PRISMA flow diagram is finalised, superseding any earlier draft that mentioned five.

# Limitations of the analysis

Several limitations should be considered when interpreting these results.

First, for several of the 12 studies in the AUC pool, the standard error was approximated using the Hanley & McNeil formula [@hanley1982] (R42, R48, R51, R52) or from a reported range (R63) rather than extracted directly. These approximations assume a binormal ROC model and may understate true uncertainty, particularly for small samples (n = 44 [@mwakapeje2019]) or non-standard validation designs.

Second, studies vary in validation strategy: some report a test AUC from independent hold-out data [@kimaro2024; @otieno2021b; @driciru2020; @kracalik2017], others a training AUC [@pittiglio2022], and others AUC from internal cross-validation. For one study [@steenkamp2018] the text reports a training AUC of 0.9 and a test AUC of 0.8, while its Table 2 reports an overall AUC of 0.857; the curated extraction (`family_A_metrics.xlsx`) uses 0.857 as the value entering the pool, and this choice should be reconfirmed against the primary publication before submission. This cross-study inconsistency introduces a degree of non-comparability.

Third, the pool combines AUC values from methodologically distinct algorithms (MaxEnt, BRT, biomod2 ensembles, Random Forest, GARP), each with different assumptions about feature classes, regularisation, and background sampling. The method subgroup analysis partially addresses this, but the overall pooled AUC should be read as a summary across the Family A literature rather than a benchmark for any single algorithm.

Fourth, two studies in the pool [@otieno2021a; @otieno2021b] share overlapping Kenyan occurrence data, introducing statistical non-independence not formally accounted for in the random-effects model; a sensitivity analysis excluding one would be informative.

Fifth, the geographic concentration in Eastern Africa limits the external validity of the pooled estimate for other African regions with distinct ecological settings.

Sixth, the corpus includes two preprints (R8, R47) that have not been peer reviewed. They were retained because they had accessible, verifiable methods, and the sensitivity analysis (Section above) confirmed the AUC conclusions are robust to excluding R47.

Seventh, seven Family A studies could not be pooled: partial-ROC ratios [@deka2022; @romeroalvarez2020], no reported discrimination metric [@dougherty2022], AUC reused from a prior publication [@gachohi2022], only threshold criteria reported [@john2024], and Bayesian INLA models validated by sensitivity rather than AUC [@ndolo2022kenya; @ndolo2022uganda]. These methodologically justified exclusions reduce the representativeness of the quantitative synthesis and are summarised narratively instead.

Eighth, the Family C OR pools rest on few studies (*K* = 4--5) with slightly varying exposure definitions (e.g. "handling" vs "skinning/butchering"), which were grouped pragmatically; each aOR should be reconfirmed against the primary publication and the groupings agreed with co-authors. The zero heterogeneity ($I^2$ = 0%), while reassuring, is estimated imprecisely at this number of studies.

Finally, this analysis is restricted to Batch 1 of the screening (63 records, 33 included). Pooled estimates may change as supplementary searches and manual reference tracking are incorporated.

# Conclusions

This partial quantitative synthesis shows that ecological niche and species distribution models for anthrax in Africa achieve a pooled discriminatory performance of AUC = 0.892 (95% CI 0.859--0.917; *K* = 12), indicating that the spatial distribution of *Bacillus anthracis* suitability can be predicted with good-to-excellent accuracy from environmental covariates. The substantial between-study heterogeneity ($I^2$ = 91%) and the wide prediction interval indicate that model performance is highly context-dependent, varying with algorithm, geographic extent, occurrence-data quality, and covariate selection.

The subgroup analyses suggest that neither algorithm choice (MaxEnt vs. BRT) nor geographic region alone accounts for the observed heterogeneity, implying that study-specific design decisions -- particularly spatial thinning, background sampling, and the balance between persistence and activation covariates -- are the primary drivers of variation. The meta-regressions on sample size and region were non-significant, and publication year showed only a borderline trend, reinforcing the interpretation that data quality and ecological context matter more than algorithmic complexity or dataset size.

In parallel, the Family C synthesis establishes that two modifiable behavioural exposures -- eating meat from dead or sick animals (pooled aOR = 6.26) and slaughtering or handling carcasses (pooled aOR = 7.55) -- drive human infection with remarkable consistency across settings ($I^2$ = 0%), providing a clear evidence base for targeted prevention.

No evidence of publication bias was detected through Egger's test or trim-and-fill, although the limited power at *K* = 12 and the reporting heterogeneity of the ENM literature warrant continued vigilance.

From the perspective of the doctoral thesis, this synthesis establishes the quantitative baseline against which Bayesian spatiotemporal models (Family B) can be evaluated. While ENM/SDM approaches excel at mapping where anthrax can persist, they are structurally limited in modelling when outbreaks emerge, because they lack explicit temporal dynamics, latent spatial processes, and the capacity to decompose environmental effects into persistence and activation components. Integrating both layers within a unified Bayesian spatiotemporal framework -- INLA-SPDE for continuous spatial fields, CHIRPS/ERA5 for dynamic activation signals, and SoilGrids for static persistence layers -- represents the natural methodological frontier identified by this review.

# References

References are numbered in Vancouver style in order of first appearance in the
text. Two conventions coexist, intentionally: in the **narrative text**, every
cited source (methodological references and primary studies) is given as a
Vancouver-numbered citation; in the **tables and forest/funnel plots**, the
included studies are additionally identified by author--year labels, because
those graphics must be readable on their own (a bracketed number would be
illegible inside a forest plot). The list below therefore contains the
methodological/reporting sources and the primary studies cited in the prose.

::: {#refs}
:::

# Session information

```{r session}
sessionInfo()
```
