--- 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() ```