Health-zone forecasts

This page splits the one-week-ahead confirmed-case forecast across the health zones, as defined in the health-zone model Methods section. The national forecast is on the forecasts page and the split by province on the province forecasts page. How these forecasts have scored against what each zone went on to report is in the forecast by health zone evaluation. The zone estimates behind the split are on the health-zone estimates page.

Load packages, data and fitted chains

# Shared setup: packages, observations and the fit registry. See
# `docs/pages/_setup.jl`.
using BVDOutbreakSize
include(joinpath(pkgdir(BVDOutbreakSize), "docs", "pages", "_setup.jl"))

validation_forecast_from (generic function with 1 method)

# The zone fit this page reads, loaded from the cache here.
chn_local = load_fit("local");

Summary

The bullets rank the zones, and the figure and table below give each zone's projection.

Split the one-week-ahead forecast across the zones

# `zone_stage_inputs` is defined in the shared setup, so the zone
# estimates page draws the same forecast from the same fixed inputs.
zone_inputs = zone_stage_inputs(; forecast = true);
zone_patch = zone_inputs.patch_of_zone;
zone_week_end = obs.cutoff + Day(7);
ZONE_THRESHOLDS = (1, 5, 10, 20)
zone_fc = zone_forecast(chn_local, zone_inputs);
zone_fc_draws = zone_forecast_draws(zone_fc, zone_inputs);
zone_fc_probs = zone_forecast_probabilities(
    zone_fc, zone_inputs; thresholds = ZONE_THRESHOLDS, draws = zone_fc_draws
);
zone_forecast_summary = zone_forecast_table(
    zone_fc, zone_inputs; thresholds = ZONE_THRESHOLDS
);
# Cases allocated to each zone over the past one, two and four weeks, and
# the last vintage on which each zone's count rose.
zone_recent = Dict(
    w => zone_recent_cases(zone_inputs; window = w) for w in (7, 14, 28)
);
zone_last_case = zone_last_case_dates(zone_inputs);
zone_forecast_bullets = let n_zones = length(zone_inputs.zone_labels)
    pct(x) = string(round(Int, 100 * x), "%")
    lead = first(
        sortperm([median(v) for v in zone_fc_draws.zones]; rev = true), 3
    )
    ranked = join(
        [
            "$(zone_inputs.zone_labels[z]) " *
                "$(median_interval_text(zone_fc_draws.zones[z]))"
                for z in lead
        ], "; "
    )
    quiet = findall(iszero, zone_recent[14])
    quiet_text = isempty(quiet) ?
        "every zone has had a case allocated over the past two weeks." :
        let z = quiet[argmax(zone_fc_probs[quiet, 1])]
            "$(length(quiet)) zones have had no case allocated over the " *
            "past two weeks. $(zone_inputs.zone_labels[z]) is the " *
            "most likely of them to report one, at " *
            "$(pct(zone_fc_probs[z, 1]))."
    end
    join(
        [
            "Projected counts are for the week to $(zone_week_end).", "",
            "- **Most new confirmed cases:** $(ranked).",
            "- **Chance of a case:** " *
                "$(count(>=(0.5), zone_fc_probs[:, 1])) of $(n_zones) " *
                "zones have at least an even chance of reporting a " *
                "confirmed case, and " *
                "$(count(>=(0.5), zone_fc_probs[:, 3])) of reporting at " *
                "least ten.",
            "- **Quiet zones:** $(quiet_text)",
        ], "\n"
    )
end;

Projected counts are for the week to 2026-10-27.

• Most new confirmed cases: Mandima about 111 (90% credible interval 47 to 210); Beni about 105 (90% credible interval 51 to 197); Bunia about 70 (90% credible interval 27 to 146).

• Chance of a case: 34 of 65 zones have at least an even chance of reporting a confirmed case, and 14 of reporting at least ten.

• Quiet zones: 21 zones have had no case allocated over the past two weeks. Viadana is the most likely of them to report one, at 44%.

One-week-ahead forecast by health zone

The patch totals in the table are the projections on the province forecasts page.

One-week-ahead zone forecast

zone_forecast_fig = plot_zone_forecast(
    zone_fc_draws.zones,
    zone_inputs.zone_labels, zone_patch;
    patch_labels = zone_inputs.patch_labels, top = 15
);
# The fifteen largest zone forecasts, then every patch total.
zone_forecast_display = let t = zone_forecast_summary
    zones = sort(t[t.zone .!= "Patch total", :], :median; rev = true)
    vcat(first(zones, 15), t[t.zone .== "Patch total", :])
end;

The fifteen zones with the largest forecasts of confirmed cases over the week after the cut-off.

zone	patch	median	lower_90	lower_60	lower_30	upper_30	upper_60	upper_90	p_ge_1	p_ge_5	p_ge_10	p_ge_20
Mandima	Ituri	111	47	74	92	132	158	210	1	1	1	1
Beni	Nord-Kivu	105	51	74	90	123	148	197	1	1	1	1
Bunia	Ituri	70	27	45	57	84	106	146	1	1	1	0.98
Rwampara	Ituri	35	9	19	27	46	60	95	1	0.99	0.95	0.8
Mambasa	Ituri	21	2	9	14	28	40	67	0.99	0.9	0.78	0.53
Katwa	Nord-Kivu	21	5	11	16	28	36	57	1	0.95	0.84	0.56
Mabalako	Nord-Kivu	20	4	10	15	27	38	60	0.99	0.93	0.81	0.52
Nizi	Ituri	17	1	7	11	23	34	58	0.98	0.87	0.71	0.44
Wamba	Haut-Uele	15	3	8	11	19	25	39	1	0.92	0.71	0.33
Butembo	Nord-Kivu	14	2	6	10	19	26	43	0.99	0.86	0.67	0.34
Oicha	Nord-Kivu	13	1	5	9	19	26	44	0.98	0.83	0.62	0.31
Kalunguta	Nord-Kivu	12	2	5	9	17	25	40	0.98	0.82	0.61	0.29
Lita	Ituri	10	0	3	6	15	22	41	0.93	0.72	0.52	0.26
Pawa	Haut-Uele	10	2	5	7	13	17	31	0.99	0.82	0.52	0.17
Komanda	Ituri	7	0	2	4	11	19	39	0.89	0.62	0.41	0.19
Patch total	Ituri	350	205	269	316	390	434	533
Patch total	Nord-Kivu	243	135	185	214	278	318	390
Patch total	Haut-Uele	33	10	19	26	40	53	80
Patch total	Sud-Kivu	0	0	0	0	0	0	2
Patch total	Tshopo	15	4	8	11	20	26	44
Patch total	Bas-Uele	1	0	0	1	2	3	7
Patch total	Sud-Ubangi	0	0	0	0	0	1	2
Patch total	Nord-Ubangi	1	0	0	0	1	3	6

The fifteen zones with the largest forecasts, then each patch total. The last four columns are the probability that the zone reports at least 1, 5, 10 and 20 confirmed cases over the week.

Zones quiet for two weeks

Quiet zones ranked by the chance of a case

zone_quiet_display = let
    quiet = findall(iszero, zone_recent[14])
    order = quiet[sortperm(zone_fc_probs[quiet, 1]; rev = true)]
    DataFrame(
        zone = zone_inputs.zone_labels[order],
        patch = zone_inputs.patch_labels[zone_patch[order]],
        cases = zone_inputs.cumulative[order],
        cases_past_4_weeks = zone_recent[28][order],
        last_case = [
            ismissing(d) ? "" : string(d) for d in zone_last_case[order]
        ],
        p_ge_1 = round.(zone_fc_probs[order, 1]; digits = 2),
        p_ge_5 = round.(zone_fc_probs[order, 2]; digits = 2),
    )
end;

zone	patch	cases	cases_past_4_weeks	last_case	p_ge_1	p_ge_5
Viadana	Bas-Uele	6	3	2026-09-26	0.44	0.02
Bulu	Sud-Ubangi	2	0	2026-09-22	0.27	0.01
Lubunga	Tshopo	1	0	2026-07-11	0.26	0.01
Wanie-Rukula	Tshopo	1	0	2026-08-01	0.24	0.01
Buta	Bas-Uele	1	0	2026-08-12	0.24	0
Fataki	Ituri	84	3	2026-10-06	0.21	0.07
Miti-Murhesa	Sud-Kivu	3	0	2026-06-01	0.19	0.01
Logo	Ituri	14	1	2026-10-05	0.18	0.06
Kilo	Ituri	38	4	2026-10-04	0.18	0.06
Drodro	Ituri	15	0	2026-09-04	0.15	0.05
Alimbongo	Nord-Kivu	5	5	2026-10-06	0.15	0.04
Boga	Ituri	2	0	2026-08-24	0.12	0.04
Kayna	Nord-Kivu	3	2	2026-10-02	0.1	0.03
Gombari	Haut-Uele	3	0	2026-08-22	0.1	0
Ariwara	Ituri	8	0	2026-08-17	0.08	0.02
Aungba	Ituri	12	0	2026-08-20	0.07	0.02
Dungu	Haut-Uele	1	0	2026-09-18	0.06	0
Mahagi	Ituri	4	0	2026-07-27	0.06	0.02
Adja	Ituri	11	1	2026-10-06	0.05	0.01
Kambala	Ituri	2	0	2026-06-30	0.04	0.01
Goma	Nord-Kivu	1	0	2026-07-02	0.04	0.01

Zones with no allocated case over the past two weeks, ranked by the probability that they report at least one, with their cases over the past four weeks and the date of their last case.

Saving zone forecast assets

The summary dashboard shows the summary bullets.

Write the zone forecast asset

dashboard_dir = joinpath(
    pkgdir(BVDOutbreakSize), "docs", "src", "summary_assets"
)
mkpath(dashboard_dir)
open(joinpath(dashboard_dir, "zone_forecast.md"), "w") do io
    print(io, zone_forecast_bullets)
end

Past zone forecasts against what was observed

Each release archives its zone forecast, and later releases score it in the forecast by health zone across releases evaluation.

Load the archived zone forecasts and their outcomes

# Written by `scripts/score_releases.jl` from each release's
# `zone_forecast.csv`, in the national overlay's schema. A missing file
# reads as an empty table, which the figure reports as nothing scored.
zone_overlay_df = _release_data(
    joinpath("zone", "forecast_overlay.csv"),
    (;
        release = String, made_date = Date, stream = String, horizon = Int,
        target_date = Date, fit = String, observed = Float64,
        median = Float64, lo30 = Float64, hi30 = Float64, lo60 = Float64,
        hi60 = Float64, lo90 = Float64, hi90 = Float64,
    )
)
zone_past_keys = zone_inputs.zone_keys[
    first(sortperm([median(v) for v in zone_fc_draws.zones]; rev = true), 10),
]
zone_past_fig = plot_forecast_overlay(
    zone_score_rows(
        scored_overlay(zone_overlay_df), zone_past_keys,
        zone_inputs.zone_labels[indexin(zone_past_keys, zone_inputs.zone_keys)]
    );
    empty_message = "No zone forecast has been scored yet. Scores appear " *
        "from the first release that archives the zone forecast onwards."
);

Past forecasts for the ten zones with the largest forecasts in this build, by the cut-off each was made from. Each shows its median and 90% predictive interval beside the persistence baseline and the count the zone went on to report. A week is left out when the zone tables carry no vintage on its last day, when it holds a harmonisation-break day, or when its province moved unallocated cases into named zones during it.

Archiving the zone forecast

Each release carries the zone forecast draws, for scoring once the week has been reported.

Write the zone forecast archive

output_dir = get(
    ENV, "BVD_OUTPUT_DIR",
    joinpath(pkgdir(BVDOutbreakSize), "output")
)
mkpath(output_dir)
CSV.write(
    joinpath(output_dir, "zone_forecast.csv"),
    zone_forecast_archive(
        zone_fc, zone_inputs; made_date = zone_inputs.cutoff, thin = 5
    )
);

---

The full analysis code, data and model definitions are in the epiforecasts/BVDOutbreakSize repository.