Modelling results

Author

acp25

Published

April 9, 2026

Purpose

This report reproduces the published results from “A simulation-based policy analysis of anticipatory action for cholera outbreaks, Democratic Republic of the Congo” by Loo, P. S., Rajah, J. K., de Leon, H. J. H., Kopainsky, B., & Milano, L. (2025). https://doi.org/10.2471/BLT.25.293226; but reimplemented as a stochastic compartmental model using the chlaa package (odin2/dust2/monty stack) located at /rds/general/user/acp25/home/MIMIC/Clean_data/CHLAA/chlaa.

Goals:

  • Reproduce the paper’s 4 intervention scenarios using chlaa
  • Fit the stochastic model to Nyiragongo outbreak data
  • Compare stochastic vs deterministic behaviour
  • Evaluate anticipatory action with uncertainty quantification
  • Assess data quality of Combined_data for multi-health-zone modelling

1. Setup

Code 
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)

library(chlaa)
library(dplyr)
library(tidyr)
library(ggplot2)

data_dir <- "/rds/general/user/acp25/home/MIMIC/Clean_data/data"

2. Paper summary and model mapping

2.1 Published model overview

The paper uses a system dynamics SIRS model (Stella Architect) for Nyiragongo health zone (N = 540,000) over Jan 2022 – Jan 2025. The model tracks susceptible, infected (asymptomatic / mild / severe, treated vs untreated), and recovered individuals with waning immunity. Transmission is driven by contaminated water contact. Eight intervention types are modelled: chlorination, hygiene kits, latrines, cholera treatment centres (CTC), oral rehydration corners (ORC), oral cholera vaccination (OCV), surveillance, and case-area targeted interventions (CATI).

2.2 The 4 scenarios

Scenario Trigger Timing Vaccination WASH/CATI Immunity
1. Baseline (reactive) Government declaration (day ~348, 14 Dec 2022) Late 1 dose, 52% coverage Excludes CATI 6 months
2. Anticipatory action 15 cases/day threshold (day ~297, 25 Oct 2022) Early 1 dose, 52% coverage All including CATI 6 months
3. AA + 1-dose OCV 15 cases/day threshold Early, vax 1 month after trigger 1 dose, 52% coverage All including CATI 6 months
4. AA + 2-dose OCV 15 cases/day threshold Early, vax 1 month after trigger 2 doses, 26% coverage All including CATI 3 years

2.3 Mapping paper parameters to chlaa

The chlaa package reimplements the same compartmental structure (S → E → A/M/Sev → treated/untreated → R with waning) as a stochastic discrete-time model in odin2/dust2. Key correspondences:

Code 
# Display chlaa parameter info for reference
param_info <- chlaa_parameter_info()
knitr::kable(
    param_info %>% select(name, default, units, description),
    caption = "chlaa model parameters"
)
chlaa model parameters
name default units description
N 5.4000e+05 persons Total population size.
E0 0.0000e+00 persons Initial exposed (latent).
A0 0.0000e+00 persons Initial asymptomatic infectious.
M0 0.0000e+00 persons Initial mild symptomatic (pre-triage stage).
Sev0 1.0000e+00 persons Initial severe symptomatic (pre-triage stage).
Mu0 0.0000e+00 persons Initial mild untreated.
Mt0 0.0000e+00 persons Initial mild treated.
Sevu0 0.0000e+00 persons Initial severe untreated.
Sevt0 0.0000e+00 persons Initial severe treated.
Ra0 0.0000e+00 persons Initial recovered after asymptomatic infection.
Rs0 0.0000e+00 persons Initial recovered after symptomatic infection.
V10 0.0000e+00 persons Initial vaccinated (1 dose protected).
V20 0.0000e+00 persons Initial vaccinated (2 dose protected).
Du0 0.0000e+00 persons Initial cumulative deaths (untreated).
Dt0 0.0000e+00 persons Initial cumulative deaths (treated).
C0 0.0000e+00 index Initial environmental contamination state.
prop_asym 7.5000e-01 probability Proportion of infections asymptomatic.
incubation_time 4.8450e+00 days Mean incubation time.
duration_asym 5.0000e+00 days Duration of asymptomatic infectiousness.
duration_sym 1.4480e+01 days Duration of symptomatic infection.
time_to_next_stage 1.0000e+00 days Time to next symptomatic stage (triage/progression).
p_progress_severe 3.0000e-01 probability Probability mild progresses to severe.
immunity_asym 2.8000e+02 days Immunity duration after asymptomatic infection.
immunity_sym 1.0950e+03 days Immunity duration after symptomatic infection.
contact_rate 1.0010e+01 contacts/person/day Effective contact rate.
trans_prob 1.2700e-01 probability/contact Per-contact transmission probability.
time_to_contaminate 1.9075e+01 days Time-scale for contamination dynamics.
water_clearance_time 3.0000e+01 days Environmental clearance time.
contam_half_sat 1.0000e+00 index Half-saturation constant for contamination effect.
shed_asym 9.0690e+04 CFU/person/day Shedding rate asymptomatic.
shed_mild 9.5005e+06 CFU/person/day Shedding rate mild.
shed_severe 3.2945e+07 CFU/person/day Shedding rate severe.
contam_scale 1.0000e+10 CFU/index Scaling from CFU to contamination index.
seek_mild 1.0000e-01 probability Care seeking probability (mild).
seek_severe 2.0000e-01 probability Care seeking probability (severe).
orc_capacity 5.0000e+02 persons/day ORC capacity while active.
ctc_capacity 1.0000e+02 persons/day CTC capacity while active.
treated_shed_mult_orc 5.0000e-01 multiplier Shedding multiplier in ORC.
treated_shed_mult_ctc 0.0000e+00 multiplier Shedding multiplier in CTC.
fatality_treated 2.1000e-03 probability Fatality probability treated (severe).
fatality_untreated 5.0000e-01 probability Fatality probability untreated (severe).
chlor_start 0.0000e+00 day Start time for chlorination.
chlor_end 0.0000e+00 day End time for chlorination.
chlor_effect 0.0000e+00 fraction Effect size for chlorination (transmission reduction).
hyg_start 0.0000e+00 day Start time for hygiene promotion.
hyg_end 0.0000e+00 day End time for hygiene promotion.
hyg_effect 0.0000e+00 fraction Effect size for hygiene (transmission reduction).
lat_start 0.0000e+00 day Start time for latrine intervention.
lat_end 0.0000e+00 day End time for latrine intervention.
lat_effect 0.0000e+00 fraction Effect size for latrines (shedding reduction).
cati_start 0.0000e+00 day Start time for CATI.
cati_end 0.0000e+00 day End time for CATI.
cati_effect 0.0000e+00 fraction Effect size for CATI (transmission reduction).
orc_start 0.0000e+00 day Start time for ORC availability.
orc_end 0.0000e+00 day End time for ORC availability.
ctc_start 0.0000e+00 day Start time for CTC availability.
ctc_end 0.0000e+00 day End time for CTC availability.
vax1_start 0.0000e+00 day Start time dose 1 campaign.
vax1_end 0.0000e+00 day End time dose 1 campaign.
vax1_doses_per_day 0.0000e+00 doses/day Dose 1 delivery rate.
vax1_total_doses 0.0000e+00 doses Total dose 1 supply.
vax2_start 0.0000e+00 day Start time dose 2 campaign.
vax2_end 0.0000e+00 day End time dose 2 campaign.
vax2_doses_per_day 0.0000e+00 doses/day Dose 2 delivery rate.
vax2_total_doses 0.0000e+00 doses Total dose 2 supply.
ve_1 4.0000e-01 fraction Efficacy after dose 1 (susceptibility reduction).
ve_2 7.0000e-01 fraction Efficacy after dose 2 (susceptibility reduction).
vax_immunity_1 1.8000e+02 days Protection duration after dose 1.
vax_immunity_2 1.0950e+03 days Protection duration after dose 2.
reporting_rate 2.0000e-01 fraction Reporting fraction for observed cases.
obs_size 2.5000e+01 size Negative binomial size parameter.

2.4 Paper parameter values vs chlaa defaults

Code 
paper_params <- tribble(
    ~parameter, ~paper_value, ~paper_units, ~chlaa_param, ~chlaa_default,
    "Population", 540000, "persons", "N", 540000,
    "Incubation time", 4.845, "days", "incubation_time", 4.845,
    "Symptomatic duration", 14.48, "days", "duration_sym", 14.48,
    "Asymptomatic duration", 5.0, "days", "duration_asym", 5.0,
    "Contact rate", 10.01, "contacts/person/day", "contact_rate", 10.01,
    "Transmission prob", 0.055, "prob/contact", "trans_prob", 0.127,
    "Prop asymptomatic", 0.75, "proportion", "prop_asym", 0.75,
    "Prop severe", 0.30, "proportion", "p_progress_severe", 0.30,
    "Shedding (asym)", 90.69e3, "CFU/person/day", "shed_asym", 90.69e3,
    "Shedding (mild)", 9.5005e6, "CFU/person/day", "shed_mild", 9.5005e6,
    "Shedding (severe)", 32.945e6, "CFU/person/day", "shed_severe", 32.945e6,
    "Water clearance", 30.0, "days", "water_clearance_time", 30.0,
    "CFR (treated)", 0.0021, "proportion", "fatality_treated", 0.0021,
    "CFR (untreated)", 0.50, "proportion", "fatality_untreated", 0.5,
    "VE dose 1", 0.40, "proportion", "ve_1", 0.4,
    "VE dose 2", 0.70, "proportion", "ve_2", 0.7,
    "Vax immunity (1 dose)", 180, "days", "vax_immunity_1", 180,
    "Vax immunity (2 dose)", 1095, "days", "vax_immunity_2", 1095
)

knitr::kable(paper_params, caption = "Paper parameters vs chlaa defaults")
Paper parameters vs chlaa defaults
parameter paper_value paper_units chlaa_param chlaa_default
Population 5.4000e+05 persons N 5.4000e+05
Incubation time 4.8450e+00 days incubation_time 4.8450e+00
Symptomatic duration 1.4480e+01 days duration_sym 1.4480e+01
Asymptomatic duration 5.0000e+00 days duration_asym 5.0000e+00
Contact rate 1.0010e+01 contacts/person/day contact_rate 1.0010e+01
Transmission prob 5.5000e-02 prob/contact trans_prob 1.2700e-01
Prop asymptomatic 7.5000e-01 proportion prop_asym 7.5000e-01
Prop severe 3.0000e-01 proportion p_progress_severe 3.0000e-01
Shedding (asym) 9.0690e+04 CFU/person/day shed_asym 9.0690e+04
Shedding (mild) 9.5005e+06 CFU/person/day shed_mild 9.5005e+06
Shedding (severe) 3.2945e+07 CFU/person/day shed_severe 3.2945e+07
Water clearance 3.0000e+01 days water_clearance_time 3.0000e+01
CFR (treated) 2.1000e-03 proportion fatality_treated 2.1000e-03
CFR (untreated) 5.0000e-01 proportion fatality_untreated 5.0000e-01
VE dose 1 4.0000e-01 proportion ve_1 4.0000e-01
VE dose 2 7.0000e-01 proportion ve_2 7.0000e-01
Vax immunity (1 dose) 1.8000e+02 days vax_immunity_1 1.8000e+02
Vax immunity (2 dose) 1.0950e+03 days vax_immunity_2 1.0950e+03

3. Reproduce the 4 scenarios (single health zone: Nyiragongo)

3.1 Define baseline parameters

Code 
# Simulation period: Jan 1, 2022 to Jan 1, 2025 = 1096 days
sim_days <- 1096
time_vec <- 0:sim_days

# Baseline parameters matching the paper
base_pars <- chlaa_parameters(
    N = 540000,
    Sev0 = 2, # Initial severe cases (paper: 1-3)
    incubation_time = 4.845,
    duration_sym = 14.48,
    duration_asym = 5.0,
    contact_rate = 10.01,
    trans_prob = 0.055, # Paper value (lower than chlaa default)
    prop_asym = 0.75,
    p_progress_severe = 0.30,
    water_clearance_time = 30.0,
    fatality_treated = 0.0021,
    fatality_untreated = 0.5,
    ve_1 = 0.40,
    ve_2 = 0.70,
    vax_immunity_1 = 180, # 6 months
    vax_immunity_2 = 1095, # 3 years
    reporting_rate = 0.2
)

3.2 Define scenario parameters

Code 
# Key time points (days from Jan 1, 2022)
# AA trigger: Oct 25, 2022 = day 297
# Government declaration: Dec 14, 2022 = day 348
aa_trigger <- 297
gov_declare <- 348
vax_delay <- 30 # vaccination starts 1 month after trigger

# Total vaccine doses for 52% and 26% coverage
doses_52pct <- round(540000 * 0.52) # ~280,800
doses_26pct <- round(540000 * 0.26) # ~140,400

# --- Scenario 1: Baseline (reactive) ---
# Triggered at government declaration; excludes CATI; 1-dose vax; 6-month immunity
sc1_reactive <- chlaa_scenario(
    name = "S1: Baseline (reactive)",
    modify = list(
        # Interventions start at government declaration
        chlor_start = gov_declare, chlor_end = gov_declare + 180,
        hyg_start = gov_declare, hyg_end = gov_declare + 180,
        lat_start = gov_declare, lat_end = gov_declare + 365,
        orc_start = gov_declare, orc_end = gov_declare + 180,
        ctc_start = gov_declare, ctc_end = gov_declare + 180,
        # No CATI
        cati_start = sim_days + 1, cati_end = sim_days + 1,
        cati_effect = 0,
        # Vaccination: 1 dose, 52% coverage, starts at declaration
        vax1_start = gov_declare, vax1_end = gov_declare + 90,
        vax1_total_doses = doses_52pct,
        # No dose 2
        vax2_start = sim_days + 1, vax2_end = sim_days + 1,
        vax2_total_doses = 0
    )
)

# --- Scenario 2: Anticipatory action ---
# Triggered at 15 cases/day (day 297); all interventions; 1-dose; 6-month immunity
sc2_aa <- chlaa_scenario(
    name = "S2: Anticipatory action",
    modify = list(
        chlor_start = aa_trigger, chlor_end = aa_trigger + 180,
        hyg_start = aa_trigger, hyg_end = aa_trigger + 180,
        lat_start = aa_trigger, lat_end = aa_trigger + 365,
        orc_start = aa_trigger, orc_end = aa_trigger + 180,
        ctc_start = aa_trigger, ctc_end = aa_trigger + 180,
        cati_start = aa_trigger, cati_end = aa_trigger + 180,
        vax1_start = aa_trigger, vax1_end = aa_trigger + 90,
        vax1_total_doses = doses_52pct,
        vax2_start = sim_days + 1, vax2_end = sim_days + 1,
        vax2_total_doses = 0
    )
)

# --- Scenario 3: AA + 1-dose OCV (vax delayed 1 month) ---
sc3_aa_1dose <- chlaa_scenario(
    name = "S3: AA + 1-dose OCV",
    modify = list(
        chlor_start = aa_trigger, chlor_end = aa_trigger + 180,
        hyg_start = aa_trigger, hyg_end = aa_trigger + 180,
        lat_start = aa_trigger, lat_end = aa_trigger + 365,
        orc_start = aa_trigger, orc_end = aa_trigger + 180,
        ctc_start = aa_trigger, ctc_end = aa_trigger + 180,
        cati_start = aa_trigger, cati_end = aa_trigger + 180,
        # Vaccination delayed by 1 month from trigger
        vax1_start = aa_trigger + vax_delay,
        vax1_end = aa_trigger + vax_delay + 90,
        vax1_total_doses = doses_52pct,
        vax2_start = sim_days + 1, vax2_end = sim_days + 1,
        vax2_total_doses = 0
    )
)

# --- Scenario 4: AA + 2-dose OCV (vax delayed 1 month, 26% coverage, 3-year immunity) ---
sc4_aa_2dose <- chlaa_scenario(
    name = "S4: AA + 2-dose OCV",
    modify = list(
        chlor_start = aa_trigger, chlor_end = aa_trigger + 180,
        hyg_start = aa_trigger, hyg_end = aa_trigger + 180,
        lat_start = aa_trigger, lat_end = aa_trigger + 365,
        orc_start = aa_trigger, orc_end = aa_trigger + 180,
        ctc_start = aa_trigger, ctc_end = aa_trigger + 180,
        cati_start = aa_trigger, cati_end = aa_trigger + 180,
        # Dose 1: starts 1 month after trigger, 26% coverage
        vax1_start = aa_trigger + vax_delay,
        vax1_end = aa_trigger + vax_delay + 90,
        vax1_total_doses = doses_26pct,
        # Dose 2: starts 2 weeks after dose 1
        vax2_start = aa_trigger + vax_delay + 14,
        vax2_end = aa_trigger + vax_delay + 104,
        vax2_total_doses = doses_26pct
    )
)

scenarios <- list(sc1_reactive, sc2_aa, sc3_aa_1dose, sc4_aa_2dose)
cat("Defined", length(scenarios), "scenarios\n")
Defined 4 scenarios

3.3 Run scenarios (deterministic)

Code 
# Run all 4 scenarios with 1 particle to approximate deterministic behaviour
# (chlaa_run_scenarios does not accept a deterministic argument)
scenario_runs_det <- chlaa_run_scenarios(
    pars = base_pars,
    scenarios = scenarios,
    time = time_vec,
    n_particles = 1,
    seed = 1
)

cat("Deterministic scenario runs complete\n")
Deterministic scenario runs complete
Code 
cat("Rows:", nrow(scenario_runs_det), "\n")
Rows: 4388
Code 
cat("Columns:", paste(names(scenario_runs_det), collapse = ", "), "\n")
Columns: scenario, time, particle, S, E, A, M, Sev, Mu, Mt, Sevu, Sevt, Ra, Rs, V1, V2, Du, Dt, C, inc_infections, inc_symptoms, inc_deaths, inc_vax1, inc_vax2, cum_infections, cum_symptoms, cum_deaths, cum_vax1, cum_vax2, cum_orc_treated, cum_ctc_treated

3.4 Scenario comparison — epidemic curves

Code 
chlaa_plot_scenario_overlay(scenario_runs_det, var = "inc_symptoms") +
    labs(
        title = "Daily symptomatic incidence by scenario (deterministic)",
        x = "Day", y = "Symptomatic cases / day"
    ) +
    theme_minimal()

Epidemic curves under 4 intervention scenarios (deterministic)

3.5 Cumulative outcomes

Code 
p_inf <- chlaa_plot_scenario_overlay(scenario_runs_det, var = "cum_infections") +
    labs(title = "Cumulative infections", x = "Day", y = "Infections") +
    theme_minimal()

p_death <- chlaa_plot_scenario_overlay(scenario_runs_det, var = "cum_deaths") +
    labs(title = "Cumulative deaths", x = "Day", y = "Deaths") +
    theme_minimal()

if (requireNamespace("patchwork", quietly = TRUE)) {
    library(patchwork)
    p_inf + p_death
} else {
    print(p_inf)
    print(p_death)
}

Cumulative infections and deaths by scenario

Cumulative infections and deaths by scenario

3.6 Scenario comparison table

Code 
comparison <- chlaa_compare_scenarios(
    scenario_runs_det,
    baseline = "S1: Baseline (reactive)"
)

knitr::kable(comparison, caption = "Scenario outcomes vs baseline (deterministic)")
Scenario outcomes vs baseline (deterministic)
scenario infections cases_symptomatic deaths doses orc_treated ctc_treated infections_averted cases_averted deaths_averted cost dalys cost_diff dalys_averted icer_cost_per_daly_averted icer_cost_per_death_averted mean_cost_vax mean_cost_care mean_cost_wash
S1: Baseline (reactive) 1784253 445449 65129 0 3819 3303 0 0 0 7133430 1955476 0 0.000 NA NA 0 302430 6831000
S2: Anticipatory action 1784002 444967 64436 0 3959 3413 251 482 693 12003630 1934684 4870200 20791.738 234.2373 7027.706 0 312630 11691000
S3: AA + 1-dose OCV 1786004 445132 64349 0 3996 3377 -1751 317 780 12001120 1932075 4867690 23401.143 208.0108 6240.628 0 310120 11691000
S4: AA + 2-dose OCV 1784491 445425 64921 0 3957 3397 -238 24 208 12002330 1949236 4868900 6240.087 780.2616 23408.173 0 311330 11691000

4. Stochastic simulations

4.1 No-intervention stochastic run

Code 
# Run stochastic baseline with no interventions to see variability
pars_no_int <- chlaa_parameters(
    N = 540000, Sev0 = 2,
    trans_prob = 0.055,
    # Push all interventions beyond simulation window (start = end to satisfy validation)
    chlor_start = sim_days + 1, chlor_end = sim_days + 1,
    hyg_start = sim_days + 1, hyg_end = sim_days + 1,
    lat_start = sim_days + 1, lat_end = sim_days + 1,
    orc_start = sim_days + 1, orc_end = sim_days + 1,
    ctc_start = sim_days + 1, ctc_end = sim_days + 1,
    cati_start = sim_days + 1, cati_end = sim_days + 1,
    vax1_start = sim_days + 1, vax1_end = sim_days + 1,
    vax2_start = sim_days + 1, vax2_end = sim_days + 1
)

sim_stoch <- chlaa_simulate(
    pars = pars_no_int,
    time = time_vec,
    n_particles = 50,
    deterministic = FALSE,
    seed = 42
)

chlaa_plot_trajectories(sim_stoch, vars = "inc_symptoms") +
    labs(
        title = "Stochastic ensemble: daily symptomatic incidence (no intervention)",
        x = "Day", y = "Cases"
    ) +
    theme_minimal()

Stochastic ensemble — no intervention (50 particles)

4.2 Stochastic scenario runs

Code 
scenario_runs_stoch <- chlaa_run_scenarios(
    pars = base_pars,
    scenarios = scenarios,
    time = time_vec,
    n_particles = 50,
    seed = 42
)

chlaa_plot_scenario_overlay(scenario_runs_stoch, var = "inc_symptoms") +
    labs(
        title = "Daily symptomatic incidence by scenario (stochastic, 50 particles)",
        x = "Day", y = "Symptomatic cases / day"
    ) +
    theme_minimal()

Stochastic scenario comparison (50 particles)

4.3 Stochastic comparison table

Code 
comparison_stoch <- chlaa_compare_scenarios(
    scenario_runs_stoch,
    baseline = "S1: Baseline (reactive)"
)

knitr::kable(comparison_stoch,
    caption = "Scenario outcomes vs baseline (stochastic, median over particles)"
)
Scenario outcomes vs baseline (stochastic, median over particles)
scenario infections cases_symptomatic deaths doses orc_treated ctc_treated infections_averted cases_averted deaths_averted cost dalys cost_diff dalys_averted icer_cost_per_daly_averted icer_cost_per_death_averted mean_cost_vax mean_cost_care mean_cost_wash
S1: Baseline (reactive) 1784546 444967.8 64478.00 0 3862.82 3308.02 0.00 0.00 0.00 7134270 1935944 0 0.0000 NA NA 0 303269.8 6831000
S2: Anticipatory action 1784944 445007.2 64462.84 0 3924.86 3357.50 -398.22 -39.40 15.16 11998849 1935490 4864579 454.6579 10699.425 320882.51 0 307848.6 11691000
S3: AA + 1-dose OCV 1785285 445096.0 64468.54 0 3928.28 3359.20 -738.58 -128.16 9.46 11999019 1935661 4864749 283.3379 17169.424 514244.08 0 308018.8 11691000
S4: AA + 2-dose OCV 1785204 444872.8 64410.58 0 3919.22 3363.78 -658.08 94.96 67.42 11999295 1933921 4865025 2022.9424 2404.925 72159.96 0 308294.6 11691000

5. Fitting the model to Nyiragongo outbreak data

5.1 Load and prepare data

Code 
# The chlaa package ships example data for Nyiragongo
nyira_file <- system.file("extdata/outbreak/nyiragongo_historical_cases_clean.csv",
    package = "chlaa"
)

if (nchar(nyira_file) > 0) {
    nyira_raw <- read.csv(nyira_file)
    cat("Nyiragongo data loaded:", nrow(nyira_raw), "rows\n")
    str(nyira_raw)
} else {
    cat("No example data found in chlaa package.\n")
    cat("Attempting to extract Nyiragongo from Combined_data...\n")
    combined <- read.csv(file.path(data_dir, "Combined_data.csv"))
    nyira_raw <- combined %>%
        filter(grepl("nyiragongo", hz, ignore.case = TRUE)) %>%
        arrange(year, numsem)
    cat("Extracted", nrow(nyira_raw), "rows for Nyiragongo\n")
}
Nyiragongo data loaded: 731 rows
'data.frame':   731 obs. of  3 variables:
 $ date  : chr  "2022-08-01" "2022-08-02" "2022-08-03" "2022-08-04" ...
 $ cases : int  6 6 6 7 7 7 7 2 1 1 ...
 $ source: chr  "Milano et al. (2025) Fig. 3 historical data (digitized and cleaned)" "Milano et al. (2025) Fig. 3 historical data (digitized and cleaned)" "Milano et al. (2025) Fig. 3 historical data (digitized and cleaned)" "Milano et al. (2025) Fig. 3 historical data (digitized and cleaned)" ...
Code 
# chlaa_prepare_data expects numeric time + integer cases
# The example data has columns: date, cases, source
# Convert date to numeric days since the start of the simulation (Jan 1, 2022)
nyira_raw$date <- as.Date(nyira_raw$date)
sim_origin <- as.Date("2022-01-01")
nyira_raw$time <- as.numeric(nyira_raw$date - sim_origin)

fit_data <- chlaa_prepare_data(
    nyira_raw,
    time_col = "time",
    cases_col = "cases",
    expected_step = 1,
    fill_missing = TRUE
)

cat("Prepared data for fitting:", nrow(fit_data), "time points\n")
Prepared data for fitting: 731 time points
Code 
cat("Time range:", min(fit_data$time), "to", max(fit_data$time), "days\n")
Time range: 212 to 942 days

5.2 Prior specification and fitting

Code 
# pMCMC fitting — computationally intensive, so eval: false by default
# Set eval: true when ready to run (consider submitting as HPC job)

fit <- chlaa_fit_pmcmc(
    data = fit_data,
    pars = base_pars,
    n_particles = 200,
    n_steps = 5000
)

# Save fit object for later use
saveRDS(fit, file.path(data_dir, "chlaa_fit_nyiragongo.rds"))

5.3 Fitting diagnostics

Code 
# Trace plots
chlaa_plot_trace(fit)

# Likelihood trace
chlaa_plot_likelihood_trace(fit)

# Posterior summary
post_summary <- chlaa_posterior_summary(fit)
knitr::kable(post_summary, caption = "Posterior parameter estimates")

# Fit report
chlaa_fit_report(fit)

5.4 Posterior predictive check

Code 
chlaa_plot_ppc(fit) +
    labs(title = "Posterior predictive check: observed vs predicted cases") +
    theme_minimal()

6. Scenario forecasts from fitted model

Code 
# Update parameters with posterior means
pars_fitted <- chlaa_update_from_fit(fit, base_pars)

# Run scenarios with fitted parameters
scenario_runs_fitted <- chlaa_run_scenarios(
    pars = pars_fitted,
    scenarios = scenarios,
    time = time_vec,
    n_particles = 500,
    seed = 42
)

# Plot
chlaa_plot_scenario_overlay(scenario_runs_fitted, var = "inc_symptoms") +
    labs(title = "Scenario forecasts from fitted model (500 particles)") +
    theme_minimal()

# Comparison with economics
comparison_fitted <- chlaa_compare_scenarios(
    scenario_runs_fitted,
    baseline = "S1: Baseline (reactive)",
    include_econ = TRUE
)

knitr::kable(comparison_fitted,
    caption = "Scenario outcomes from fitted model (with cost-effectiveness)"
)

7. Health economic evaluation

Code 
# Cost-effectiveness analysis
econ <- chlaa_health_econ(scenario_runs_fitted, econ = chlaa_econ_defaults())

knitr::kable(econ, caption = "Health economic outcomes by scenario")

# Cost-effectiveness plane
chlaa_plot_ce_plane(scenario_runs_fitted, baseline = "S1: Baseline (reactive)")

# Cost-effectiveness acceptability curve
chlaa_ceac(scenario_runs_fitted, baseline = "S1: Baseline (reactive)")

8. Extending to multi-health-zone modelling

8.1 Assess Combined_data suitability

Code 
combined <- read.csv(file.path(data_dir, "Combined_data.csv"))

cat("Combined_data dimensions:", nrow(combined), "rows x", ncol(combined), "cols\n\n")
Combined_data dimensions: 419098 rows x 13 cols
Code 
# Health zones with case data
hz_with_cases <- combined %>%
    group_by(hz) %>%
    summarise(
        total_cases = sum(cases, na.rm = TRUE),
        total_deaths = sum(deaths, na.rm = TRUE),
        n_weeks = n(),
        years = paste(sort(unique(year)), collapse = ", "),
        has_population = any(!is.na(population)),
        .groups = "drop"
    ) %>%
    arrange(desc(total_cases))

cat(
    "Health zones with any cases:", sum(hz_with_cases$total_cases > 0),
    "of", nrow(hz_with_cases), "\n"
)
Health zones with any cases: 384 of 428
Code 
cat(
    "Health zones with population data:", sum(hz_with_cases$has_population),
    "of", nrow(hz_with_cases), "\n\n"
)
Health zones with population data: 400 of 428
Code 
knitr::kable(head(hz_with_cases, 20),
    caption = "Top 20 health zones by total cholera cases"
)
Top 20 health zones by total cholera cases
hz total_cases total_deaths n_weeks years has_population
goma 35493 71 1041 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024, 2025 TRUE
fizi 29448 60 1042 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024, 2025 TRUE
uvira 23922 84 1042 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024, 2025 TRUE
kalemie 23571 208 1042 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024, 2025 TRUE
karisimbi 23021 69 1041 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024, 2025 TRUE
kirotshe 22231 162 1042 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024, 2025 TRUE
nyiragongo 20604 30 1041 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024, 2025 TRUE
nyemba 19773 159 1042 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024, 2025 TRUE
minova 17240 101 1042 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024, 2025 TRUE
mweso 11882 69 1042 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024, 2025 TRUE
kinkondja 11725 320 1041 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024, 2025 TRUE
kadutu 10461 67 1042 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024, 2025 TRUE
rutshuru 10226 32 1042 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024, 2025 TRUE
katana 9809 35 1042 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024, 2025 TRUE
moba 9775 207 1042 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024, 2025 TRUE
masisi 9229 85 1042 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024, 2025 TRUE
ruzizi 8360 53 1042 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024, 2025 TRUE
bukama 8300 144 1042 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024, 2025 TRUE
malemba 6999 192 561 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024, 2025 FALSE
kenya 6937 108 1040 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024, 2025 TRUE

8.2 Data gaps for multi-HZ modelling

Code 
# Identify data gaps that would need addressing for multi-HZ modelling
cat("=== Key data gaps for multi-HZ modelling ===\n\n")
=== Key data gaps for multi-HZ modelling ===
Code 
cat("1. Population denominators:\n")
1. Population denominators:
Code 
cat("   HZs missing population:", sum(!hz_with_cases$has_population), "\n\n")
   HZs missing population: 28
Code 
cat("2. Temporal coverage per HZ:\n")
2. Temporal coverage per HZ:
Code 
hz_temporal <- combined %>%
    group_by(hz) %>%
    summarise(
        min_year = min(year, na.rm = TRUE),
        max_year = max(year, na.rm = TRUE),
        n_years = n_distinct(year),
        .groups = "drop"
    )
cat(
    "   HZs with full 2017-2024 coverage:",
    sum(hz_temporal$n_years == 8), "of", nrow(hz_temporal), "\n\n"
)
   HZs with full 2017-2024 coverage: 3 of 428
Code 
cat("3. Spatial connectivity data:\n")
3. Spatial connectivity data:
Code 
cat("   Not available in current dataset — would need adjacency matrix\n")
   Not available in current dataset — would need adjacency matrix
Code 
cat("   or mobility data for spatial coupling between HZs.\n")
   or mobility data for spatial coupling between HZs.

9. Next steps

NoteImmediate priorities
  1. Validate deterministic reproduction (Section 3): compare epidemic curve shapes and cumulative outcomes against published Fig. 2-3
  2. Tune trans_prob: paper uses 0.055 but chlaa default is 0.127 — calibrate to match observed outbreak dynamics
  3. Run pMCMC fitting (Section 5): submit as HPC job with n_particles = 200, n_steps = 10000

Medium-term

  1. Sensitivity analysis: replicate paper’s local (15%) and global (100k simulations) sensitivity analyses using chlaa_counterfactual_grid()
  2. Multi-HZ extension: fit independently to top-burden health zones from Combined_data
  3. Seasonal forcing: paper acknowledges this limitation — could add seasonal transmission multiplier

Longer-term

  1. Spatial coupling: connect health zones via mobility-weighted force of infection
  2. Real-time forecasting: use chlaa_forecast_from_fit() with rolling data windows
  3. Policy brief: use chlaa_compare_scenarios() + health economics to produce decision-support outputs