BVDOutbreakSize

Appearance

Health-zone estimates ​

Each patch of the headline joint fit split across its health zones: the reproduction number, the share of the patch, the one-week forecast and the probability of at least a few cases, zone by zone. The health-zone model on the methods page gives the maths. This page carries its results, the checks of the fit and the interactive map. The change in each zone's estimate over the past week has its own section below. The one-week zone forecast is on the health-zone forecasts page and its scores in the health-zone forecast evaluation.

Load packages, data and fitted chains
julia
# 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)
julia
# The fits this page reads, loaded from the cache here: the headline joint,
# the zone fit melded from it, and the frozen zone fit the week-on-week
# comparison reads.
chn_joint = load_fit("joint");
chn_local = load_fit("local");
frozen_local = load_fit("local_frozen_validation");

# The zone stage's fixed inputs with the joint's forecast, and the frozen
# fit's inputs for the comparison below, both from the shared setup, so the
# health-zone forecast page draws the same forecast from the same inputs.
zone_inputs = zone_stage_inputs(; forecast = true);
zone_patch = zone_inputs.patch_of_zone;
frozen_zone_inputs = frozen_zone_stage_inputs();

Last updated: 30 September 2026.

Data as of: 26 September 2026.

Estimates by zone ​

The maps below show, for the health-zone model, the reproduction number at the cut-off, the bounds of the 90% interval on the forecast confirmed cases over the coming week, and the confirmed cases to date, zone by zone. The forecast is mapped as its two bounds rather than a single number, so a zone's colour reads as a range. On the reproduction-number map a zone whose 90% interval straddles one is washed towards white. Zones with no confirmed case, or too few infections for a reproduction number, are grey. The four maps share the zone boundaries. A zone's forecast can be read against its reproduction number and its cases to date. The interactive map below adds the probability of at least cases at a chosen and a filter for zones with or without a case over the past one, two or four weeks.

Health-zone post-processing
julia
# The geojson keys a zone without the province prefix the manifest carries.
zone_map_keys = [
    String(last(split(k, "."; limit = 2)))
        for k in zone_inputs.zone_keys
];
# Daily zone reproduction numbers over the zone grid, each zone draw on
# its own draw of the patch trajectory so the patch uncertainty is
# carried, and the cut-off values the map and ranking read. A zone's
# reproduction number is reported from the day its cumulative infections
# reach the floor in the median draw.
zone_rt_traj = reconstruct_zone_rt(chn_local, zone_inputs);
zone_RT_finite = [filter(isfinite, m[:, obs.n]) for m in zone_rt_traj];
zone_RT_reported = findall(!isempty, zone_RT_finite);
zone_share_T = let vs = vec(collect(chn_local[:share_T_zone]))
    [Float64[v[z] for v in vs] for z in eachindex(zone_map_keys)]
end;
_zq(v, p) = quantile(v, p)
# The one-week zone forecast drawn from the zone model, and the
# probability of at least K cases per zone from the same draws.
zone_fc = zone_forecast(chn_local, zone_inputs);
zone_fc_draws = zone_forecast_draws(zone_fc, zone_inputs);
ZONE_THRESHOLDS = (1, 5, 10, 20)
zone_fc_probs = zone_forecast_probabilities(
    zone_fc, zone_inputs; thresholds = ZONE_THRESHOLDS, draws = zone_fc_draws
);
zone_overview = zone_overview_table(chn_local, zone_inputs);
# 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_map_fig = plot_zone_map_panels(
    [
        (;
            values = [median(zone_RT_finite[z]) for z in zone_RT_reported],
            zones = zone_map_keys[zone_RT_reported],
            lower = [_zq(zone_RT_finite[z], 0.05) for z in zone_RT_reported],
            upper = [_zq(zone_RT_finite[z], 0.95) for z in zone_RT_reported],
            diverging_at = 1.0, scale = log10,
            title = "Reproduction number at the cut-off",
            colorbar_label = "R",
        ),
        (;
            values = [_zq(v, 0.05) for v in zone_fc_draws.zones],
            zones = zone_map_keys, scale = CairoMakie.Makie.pseudolog10,
            title = "Confirmed cases over the coming week (lower 90%)",
            colorbar_label = "cases",
        ),
        (;
            values = [_zq(v, 0.95) for v in zone_fc_draws.zones],
            zones = zone_map_keys, scale = CairoMakie.Makie.pseudolog10,
            title = "Confirmed cases over the coming week (upper 90%)",
            colorbar_label = "cases",
        ),
        (;
            values = Float64.(zone_inputs.cumulative), zones = zone_map_keys,
            scale = CairoMakie.Makie.pseudolog10,
            title = "Confirmed cases to date", colorbar_label = "cases",
        ),
    ];
    ncols = 2, title = "Health zones at the cut-off"
);

The panels below trace the reproduction number of the twelve zones with most confirmed cases, each against its patch's own implied reproduction number in grey. Where a zone's line departs from the grey patch line, the gap is the zone's fitted deviation from its patch.

Zone reproduction-number trajectories
julia
# Each patch's implied reproduction number from the joint draws with the
# same generation interval the zone stage fixes, so the grey reference is
# the quantity the zone values average to.
zone_grid = zone_inputs.t0:obs.n
_patch_infection_draws = vec(collect(chn_joint[:infections_patch]));
patch_implied_rt = [
    let m = Matrix{Float64}(
            undef,
            length(_patch_infection_draws), obs.n
        )
        for (i, v) in enumerate(_patch_infection_draws)
            I = reshape(Float64.(v), N_PATCHES, obs.n)
            m[i, :] .= implied_national_Rt(I[p, :], zone_inputs.g)
        end
        m
    end
        for p in 1:N_PATCHES
];
zone_rt_fig = plot_rt_zones(
    [replace(m[:, zone_grid], NaN => missing) for m in zone_rt_traj],
    zone_inputs.zone_labels, zone_patch;
    patch_labels = zone_inputs.patch_labels,
    dates = grid_date.(zone_grid), as_of_date = obs.cutoff,
    cumulative = zone_inputs.cumulative, top = 12,
    patch_rt = [m[:, zone_grid] for m in patch_implied_rt]
);

The ranking below orders the zones by the posterior probability that their reproduction number exceeds one. A zone whose reproduction number is from its province, not modelled separately, is drawn hollow in grey.

Zone ranking
julia
zone_ranking_fig = plot_zone_ranking(
    zone_overview;
    patch_labels = zone_inputs.patch_labels
);
# The overview as displayed: the interval strings, without the numeric
# columns the figure reads.
zone_overview_display = let d = zone_overview[
        :,
        [:zone, :patch, :cases, :share, :R_T, :p_R_above_1, :delta_T],
    ]
    d[!, "R modelled separately"] = [
        w ? "yes" : "no" for w in zone_overview.walking
    ]
    d
end;

The table gives the twenty highest-ranked zones: the confirmed cases to date, the zone's share of its patch's infections at the cut-off in percent, its reproduction number, the probability that it exceeds one, its log-transmission deviation at the cut-off and whether its reproduction number is modelled separately. A zone whose reproduction number is not modelled separately takes it from its province. The share, the reproduction number and the deviation are each a median with a 90% interval. Every zone is listed in the fold below it.

zonepatchcasesshareR_Tp_R_above_1delta_TR modelled separately
MandimaIturi11850.7 (25.2–73.9)1.85 (1.21–2.63)0.981.38 (0.72–2.09)yes
BeniNord-Kivu32961.6 (43.6–76.2)1.03 (0.69–1.44)0.560.5 (0.08–0.97)yes
MabalakoNord-Kivu316.0 (2.4–14.5)1.06 (0.59–1.85)0.560.49 (0.05–1.03)yes
PawaHaut-Uele7748.6 (27.9–69.3)1.02 (0.61–1.54)0.530.3 (-0.02–0.66)yes
KalungutaNord-Kivu516.0 (2.8–13.6)0.93 (0.5–1.59)0.420.34 (-0.03–0.74)yes
WambaHaut-Uele10929.1 (13.9–51.6)0.94 (0.52–1.53)0.420.21 (-0.13–0.65)yes
Boma MangbetuHaut-Uele3910.0 (4.6–21.5)0.94 (0.55–1.51)0.410.15 (-0.16–0.5)yes
Nia-NiaIturi2848.7 (3.2–21.5)0.77 (0.34–1.41)0.260.47 (-0.14–1.15)yes
KyondoNord-Kivu331.6 (0.7–3.8)0.73 (0.4–1.28)0.18-0.02 (-0.5–0.43)yes
KomandaIturi1666.3 (2.5–15.0)0.63 (0.32–1.18)0.120.28 (-0.28–0.85)yes
LolwaIturi330.4 (0.1–1.5)0.59 (0.3–1.14)0.10.14 (-0.42–0.77)yes
DamasIturi400.2 (0.1–0.8)0.53 (0.26–1.03)0.06-0.11 (-0.71–0.49)yes
MongbwaluIturi6975.0 (2.0–12.1)0.55 (0.26–1.04)0.060.14 (-0.39–0.68)yes
NyankundeIturi1270.2 (0.1–0.6)0.58 (0.32–1.0)0.05-0.14 (-0.72–0.41)yes
IsiroHaut-Uele1007.2 (2.8–16.4)0.55 (0.28–1.0)0.05-0.39 (-0.95–0.03)yes
KiloIturi350.1 (0.0–0.4)0.55 (0.29–0.98)0.04-0.13 (-0.61–0.37)yes
TchomiaIturi760.4 (0.1–1.3)0.52 (0.25–0.95)0.04-0.04 (-0.54–0.43)yes
MusieneneNord-Kivu1212.0 (0.9–4.5)0.52 (0.27–0.94)0.04-0.33 (-0.81–0.1)yes
BambuIturi2011.3 (0.5–3.4)0.47 (0.23–0.9)0.03-0.05 (-0.46–0.37)yes
RwamparaIturi10783.4 (1.4–8.4)0.5 (0.24–0.9)0.030.0 (-0.4–0.41)yes
All health zones
zonepatchcasesshareR_Tp_R_above_1delta_TR modelled separately
MandimaIturi11850.7 (25.2–73.9)1.85 (1.21–2.63)0.981.38 (0.72–2.09)yes
BeniNord-Kivu32961.6 (43.6–76.2)1.03 (0.69–1.44)0.560.5 (0.08–0.97)yes
MabalakoNord-Kivu316.0 (2.4–14.5)1.06 (0.59–1.85)0.560.49 (0.05–1.03)yes
PawaHaut-Uele7748.6 (27.9–69.3)1.02 (0.61–1.54)0.530.3 (-0.02–0.66)yes
KalungutaNord-Kivu516.0 (2.8–13.6)0.93 (0.5–1.59)0.420.34 (-0.03–0.74)yes
WambaHaut-Uele10929.1 (13.9–51.6)0.94 (0.52–1.53)0.420.21 (-0.13–0.65)yes
Boma MangbetuHaut-Uele3910.0 (4.6–21.5)0.94 (0.55–1.51)0.410.15 (-0.16–0.5)yes
Nia-NiaIturi2848.7 (3.2–21.5)0.77 (0.34–1.41)0.260.47 (-0.14–1.15)yes
KyondoNord-Kivu331.6 (0.7–3.8)0.73 (0.4–1.28)0.18-0.02 (-0.5–0.43)yes
KomandaIturi1666.3 (2.5–15.0)0.63 (0.32–1.18)0.120.28 (-0.28–0.85)yes
LolwaIturi330.4 (0.1–1.5)0.59 (0.3–1.14)0.10.14 (-0.42–0.77)yes
DamasIturi400.2 (0.1–0.8)0.53 (0.26–1.03)0.06-0.11 (-0.71–0.49)yes
MongbwaluIturi6975.0 (2.0–12.1)0.55 (0.26–1.04)0.060.14 (-0.39–0.68)yes
NyankundeIturi1270.2 (0.1–0.6)0.58 (0.32–1.0)0.05-0.14 (-0.72–0.41)yes
IsiroHaut-Uele1007.2 (2.8–16.4)0.55 (0.28–1.0)0.05-0.39 (-0.95–0.03)yes
KiloIturi350.1 (0.0–0.4)0.55 (0.29–0.98)0.04-0.13 (-0.61–0.37)yes
TchomiaIturi760.4 (0.1–1.3)0.52 (0.25–0.95)0.04-0.04 (-0.54–0.43)yes
MusieneneNord-Kivu1212.0 (0.9–4.5)0.52 (0.27–0.94)0.04-0.33 (-0.81–0.1)yes
BambuIturi2011.3 (0.5–3.4)0.47 (0.23–0.9)0.03-0.05 (-0.46–0.37)yes
RwamparaIturi10783.4 (1.4–8.4)0.5 (0.24–0.9)0.030.0 (-0.4–0.41)yes
ButemboNord-Kivu2564.4 (2.1–8.5)0.51 (0.27–0.91)0.03-0.31 (-0.64–0.01)yes
FatakiIturi820.3 (0.1–0.8)0.46 (0.24–0.87)0.02-0.31 (-0.9–0.24)yes
LitaIturi2832.3 (0.9–5.7)0.46 (0.22–0.89)0.02-0.05 (-0.46–0.39)yes
KatwaNord-Kivu5557.1 (3.6–13.2)0.48 (0.26–0.84)0.02-0.39 (-0.72–-0.09)yes
BuniaIturi16978.6 (3.4–17.4)0.44 (0.2–0.8)0.01-0.1 (-0.48–0.27)yes
MangalaIturi3462.3 (0.9–5.4)0.4 (0.18–0.79)0.01-0.17 (-0.61–0.22)yes
NiziIturi7711.9 (0.8–4.6)0.38 (0.17–0.74)0.01-0.25 (-0.63–0.16)yes
MangoboOther provinces59.6 (5.6–17.1)1.26 (0.83–1.82)0.820.04 (0.0–0.12)no
Makiso-KisanganiOther provinces2343.4 (30.8–57.4)1.25 (0.84–1.81)0.810.08 (0.03–0.19)no
KabondoOther provinces57.2 (3.5–13.7)1.24 (0.82–1.76)0.80.03 (-0.02–0.1)no
BafwasendeOther provinces76.9 (2.9–15.1)1.22 (0.81–1.73)0.780.04 (-0.02–0.14)no
ViadanaOther provinces66.6 (2.8–14.7)1.2 (0.81–1.7)0.770.02 (-0.05–0.1)no
MutwangaNord-Kivu120.9 (0.5–1.6)0.86 (0.55–1.27)0.26-0.05 (-0.18–-0.0)no
KaynaNord-Kivu10.3 (0.2–0.6)0.84 (0.55–1.24)0.24-0.13 (-0.37–-0.03)no
OichaNord-Kivu241.6 (1.0–2.6)0.82 (0.53–1.23)0.210.0 (-0.06–0.05)no
BienaNord-Kivu180.8 (0.4–1.4)0.76 (0.47–1.15)0.150.02 (-0.03–0.07)no
LuberoNord-Kivu211.1 (0.6–1.8)0.77 (0.48–1.15)0.15-0.01 (-0.1–0.03)no
MaserekaNord-Kivu281.3 (0.7–2.2)0.76 (0.47–1.14)0.13-0.0 (-0.06–0.04)no
VuhoviNord-Kivu291.5 (0.8–2.5)0.75 (0.47–1.13)0.130.02 (-0.01–0.08)no
AdjaIturi110.1 (0.0–0.1)0.76 (0.48–1.1)0.12-0.1 (-0.33–-0.02)no
AriwaraIturi80.2 (0.1–0.3)0.77 (0.48–1.1)0.12-0.09 (-0.31–-0.01)no
AruIturi80.2 (0.1–0.3)0.76 (0.48–1.1)0.12-0.12 (-0.35–-0.03)no
MahagiIturi40.2 (0.1–0.2)0.75 (0.48–1.08)0.11-0.12 (-0.37–-0.03)no
ManguredjipaNord-Kivu130.7 (0.3–1.4)0.72 (0.44–1.09)0.10.06 (0.01–0.16)no
AungbaIturi120.1 (0.1–0.2)0.74 (0.46–1.07)0.09-0.11 (-0.31–-0.03)no
BogaIturi20.1 (0.1–0.2)0.71 (0.43–1.04)0.08-0.02 (-0.12–0.04)no
KambalaIturi20.1 (0.1–0.2)0.71 (0.46–1.03)0.07-0.1 (-0.27–-0.02)no
RimbaIturi120.2 (0.1–0.3)0.72 (0.45–1.03)0.07-0.1 (-0.32–-0.02)no
LogoIturi130.2 (0.1–0.3)0.7 (0.43–1.01)0.06-0.08 (-0.26–-0.02)no
GetyIturi70.2 (0.1–0.4)0.69 (0.43–1.0)0.05-0.04 (-0.16–0.01)no
MambasaIturi270.3 (0.2–0.6)0.67 (0.42–1.0)0.050.02 (-0.02–0.1)no
DrodroIturi150.3 (0.2–0.6)0.62 (0.37–0.93)0.02-0.02 (-0.12–0.02)no
GomaNord-Kivu10.2 (0.1–0.3)NaN-0.11 (-0.34–-0.02)no
DunguHaut-Uele10.5 (0.3–0.9)NaN-0.12 (-0.31–-0.04)no
GombariHaut-Uele30.5 (0.2–1.1)NaN-0.05 (-0.17–-0.0)no
RunguHaut-Uele20.9 (0.5–1.5)NaN-0.07 (-0.23–-0.01)no
Miti-MurhesaOther provinces32.8 (1.0–6.3)NaN-0.08 (-0.25–-0.01)no
LubungaOther provinces13.9 (1.8–7.7)NaN0.01 (-0.08–0.07)no
TshopoOther provinces12.5 (1.1–5.5)NaN-0.03 (-0.17–0.03)no
Wanie-RukulaOther provinces12.4 (1.0–5.6)NaN-0.03 (-0.16–0.04)no
ButaOther provinces12.6 (1.0–5.4)NaN-0.05 (-0.24–0.02)no
GangaOther provinces34.3 (1.7–9.6)NaN-0.0 (-0.11–0.07)no
BuluOther provinces22.2 (0.6–7.2)NaN0.03 (-0.08–0.12)no

Composition checks ​

Whether the model reproduces each zone's observed share of its patch's confirmed cases and deaths is on the in-sample checks page.

Data currency ​

The zone blocks are the per-zone confirmed-case and confirmed-death tables of the situation reports. Each is read here against the cut-off, so a block that stops being picked up shows as a date before it rather than as a flat series.

Zone block currency
julia
# Every zone block whose last vintage falls before the cut-off. The grace
# the national table allows is not applied here: one vintage behind is
# already worth reading.
zone_currency = let
    status = stream_report_status(obs; stratum = :zone)
    behind = status[[ismissing(d) || d > 0 for d in status.days_since], :]
    isempty(behind) ?
        Markdown.parse(
            "Every zone block reports to the cut-off, $(obs.cutoff)."
        ) :
        MarkdownTable(
            DataFrame(
                "Block" => behind.label,
                "Last reported" => [
                    ismissing(d) ? "never" : string(d)
                    for d in behind.last_date
                ],
                "Days before cut-off" => [
                    ismissing(d) ? "-" : string(d)
                    for d in behind.days_since
                ]
            )
        )
end;

Every zone block reports to the cut-off, 2026-09-26.

Health-zone fit diagnostics ​

The table gives the sampler diagnostics of the two zone fits: the worst R-hat and smallest effective sample sizes over every stored quantity but the zone reproduction number, and the divergences. Per chain it gives the fraction of iterations at the tree-depth cap, the energy fraction of missing information and the adapted step size. The per-zone R-hat and effective sample sizes of the cut-off reproduction number, share and deviation are in the fold. The rows whose reproduction number is modelled separately are the ones to read.

Zone fit diagnostics
julia
_zone_per_chain(v) = join(string.(round.(v; sigdigits = 3)), " / ")
function _zone_sampler_row(label, chn)
    d = zone_sampler_diagnostics(chn)
    return (
        fit = label, max_rhat = round(d.max_rhat; digits = 3),
        min_ess_bulk = round(d.min_ess_bulk; digits = 0),
        min_ess_tail = round(d.min_ess_tail; digits = 0),
        divergences = d.n_divergent,
        depth_cap = _zone_per_chain(d.depth_cap_fraction),
        ebfmi = _zone_per_chain(d.ebfmi),
        step_size = _zone_per_chain(d.step_size),
    )
end
zone_sampler_table = DataFrame(
    [
        _zone_sampler_row("health zones", chn_local),
        _zone_sampler_row("health zones (frozen)", frozen_local.chn),
    ]
);
zone_diagnostics = let d = zone_diagnostics_table(chn_local, zone_inputs)
    for c in names(d)[3:end]
        d[!, c] = round.(d[!, c]; digits = startswith(c, "rhat") ? 3 : 0)
    end
    DataFrame(
        [
            n == "walking" ?
                "R modelled separately" => [w ? "yes" : "no" for w in d[!, n]] :
                n => d[!, n]
                for n in names(d)
        ]
    )
end;
fitmax_rhatmin_ess_bulkmin_ess_taildivergencesdepth_capebfmistep_size
health zones1.02714610900.0 / 0.00.857 / 0.9020.0099 / 0.0104
health zones (frozen)1.01921425100.0 / 0.00.966 / 0.9170.00963 / 0.0113
Per-zone R-hat and effective sample sizes
zoneR modelled separatelyrhat_R_Tess_bulk_R_Tess_tail_R_Trhat_share_Tess_bulk_share_Tess_tail_share_Trhat_delta_Tess_bulk_delta_Tess_tail_delta_T
AdjanoNaNNaNNaN1125212681.00313641121
AriwaranoNaNNaNNaN1.001122610951.00113891434
ArunoNaNNaNNaN1112812070.99910981275
AungbanoNaNNaNNaN1.001120811791.0011063996
Bambuyes1230412271214014301.00129401341
BoganoNaNNaNNaN1.001143913451.00118761013
Buniayes1196110051.0011444938126981030
Damasyes1.001270212081.00117941508132701176
Drodrono1.008195212891.00114401379111761160
Fatakiyes1.002200313861.001175914971.00226351337
GetynoNaNNaNNaN1.00112281194114861166
KambalanoNaNNaNNaN1.0011286120319391111
Kiloyes1.003186412141.001181611981.00234051293
Komandayes1245810211203014071.00129511449
Litayes1214010971.001210412471.00527131224
Logono1.008182113281.003123713931.0018571041
Lolwayes1274911901.01195814911.00135101372
MahaginoNaNNaNNaN1.001114712421.0019661184
Mambasano1.003175013211.00313738741.00114381073
Mandimayes1198911951.00119151367121721216
Mangalayes1215114401.001180013411.00326311128
Mongbwaluyes0.999223611331.002203812371.00229821169
Nia-Niayes1.003277910380.999258814450.99926331167
Niziyes1238613491.004133813390.99934211083
Nyankundeyes1.00224881337120141224137301203
RimbanoNaNNaNNaN1.00212489841.003905888
Rwamparayes0.999218111831.00588410201.00325971078
Tchomiayes1.003193610861.00512291346129261155
Beniyes1.001322312711130513820.99923801307
Bienano1262111041.0011219127711360826
Butemboyes1324811081.001215811001.00133301178
GomanoNaNNaNNaN1.0016454951.00114511130
Kalungutayes1.00129649321.001176414850.99925851186
Katwayes0.999289912681.003190412431.00126531187
KaynanoNaNNaNNaN16657391.00112281248
Kyondoyes1.00131899701.001219617171.00231661304
Luberono125769911.00110158311.00312731092
Mabalakoyes1231012670.999248913191.0071478642
ManguredjipanoNaNNaNNaN1.002146514311.00310271344
Maserekano1320211371.0027219791.0021394966
Musieneneyes1.001278313390.9991853147912468920
MutwanganoNaNNaNNaN0.99971310080.99912741116
Oichano122479601.00290012671.00215941300
Vuhovino1.00128518851.00178211281.00211751170
Boma Mangbetuyes1.001383510781.001225713041.00129491119
DungunoNaNNaNNaN1.001151015001.00111471275
GombarinoNaNNaNNaN11675128111627860
Isiroyes1271512351.001211110571.00129761080
Pawayes1.006291513831.001201812491.00524301283
RungunoNaNNaNNaN1.00115671358112321167
Wambayes1390311871.00320891312129941215
Miti-MurhesanoNaNNaNNaN1155114251.00213971180
BafwasendenoNaNNaNNaN1.003145710241.0021397990
KabondonoNaNNaNNaN0.99921321429116751055
LubunganoNaNNaNNaN118051169117641201
Makiso-Kisanganino1.002391810661.002187712511.00310241156
MangobonoNaNNaNNaN0.999205615621.00214131057
TshoponoNaNNaNNaN1186712651.00116691230
Wanie-RukulanoNaNNaNNaN1202512901.00120301126
ButanoNaNNaNNaN1176215251.00219641330
GanganoNaNNaNNaN119191363116991134
ViadananoNaNNaNNaN1.001173614230.99911861303
BulunoNaNNaNNaN0.999151512130.99912011254

The figure below sets the reproduction number implied by the zone stage's own patch trajectories against the one implied by the headline joint fit, nationally and for each patch. Both are computed from infections with the generation interval the zone stage fixes. Agreement says the melding stage has kept the joint's patch trajectories rather than moved them to fit the zone data.

Reproduction number from the zone stage and the joint fit
julia
# The implied reproduction number of each row of `draws` (draws × days).
function _implied_rt_matrix(draws::AbstractMatrix)
    m = similar(draws, Float64)
    for i in axes(draws, 1)
        m[i, :] .= implied_national_Rt(view(draws, i, :), zone_inputs.g)
    end
    return m
end
# The zone stage's deformed patch trajectories on one side and the joint's
# draws on the other, nationally then per patch. The joint's per-patch
# values are the grey references of the trajectory panels above.
zone_stage_rt = let I = zone_patch_infections(chn_local, zone_inputs)
    vcat([_implied_rt_matrix(sum(I))], _implied_rt_matrix.(I))
end;
joint_stage_rt = let draws = _patch_infection_draws
    national = Float64[
        sum(reshape(Float64.(v), N_PATCHES, obs.n)[:, t])
            for v in draws, t in 1:obs.n
    ]
    vcat([_implied_rt_matrix(national)], patch_implied_rt)
end;
zone_meld_rt_fig = plot_rt_zones(
    [m[:, zone_grid] for m in zone_stage_rt],
    vcat(["National"], zone_inputs.patch_labels),
    vcat([N_PATCHES + 1], 1:N_PATCHES);
    patch_labels = vcat(zone_inputs.patch_labels, ["National"]),
    patch_colours = [:firebrick, :steelblue, :seagreen, :darkorange, :black],
    dates = grid_date.(zone_grid), as_of_date = obs.cutoff,
    top = N_PATCHES + 1, ncols = 3,
    reference_rt = [m[:, zone_grid] for m in joint_stage_rt],
    reference_label = "Headline joint fit",
    title = "Reproduction number from the zone stage and the joint fit"
);

Health-zone parameters against their priors ​

The table sets the posterior of each zone hyperparameter against its prior, with the ratio of their standard deviations. A ratio near one says the zone data add little to the prior. A ratio above one says the posterior is wider than the prior, which happens when the data move a parameter into the prior's wider tail. The drift scale is one per patch. The pair plot overlays the prior on the posterior of the scalar hyperparameters.

Draw from the zone model's prior
julia
prior_chn_zone = zone_prior_draws(zone_inputs);
Compute the zone prior and posterior table
julia
# Each scalar zone hyperparameter with its table label and pair-plot axis
# label, kept where both chains carry it: the mixing and correlation
# blocks are sampled only when their inputs are on.
zone_hyper = [
    h for h in (
            (:region_sd_zone, "Level spread σ_level", "σ_level"),
            (:region_halflife_zone, "Deviation half-life (days)", "half-life (days)"),
            (:correlation_reference_zone, "Correlation ρ_corr", "ρ_corr"),
            (
                :zone_ascertainment_sd, "Ascertainment spread σ_ascertainment",
                "σ_ascertainment",
            ),
            (:zone_severity_sd, "Severity spread σ_severity", "σ_severity"),
            (:mixing_within_zone, "Within-patch mixing ε_within", "ε_within"),
            (:mixing_departure_zone, "Mixing departure τ_mix", "τ_mix"),
            (:composition_rho_zone, "Case composition ρ", "ρ"),
            (:composition_rho_death_zone, "Death composition ρ_death", "ρ_death"),
        )
        if BVDOutbreakSize._has_key(chn_local, h[1]) &&
        BVDOutbreakSize._has_key(prior_chn_zone, h[1])
]
zone_hyper_keys = first.(zone_hyper)
_hyper_draws(chn, k) = Float64.(vec(collect(chn[k])))
_drift_draws(chn, p) = Float64[
    v[p] for v in vec(collect(chn[:region_drift_sd_zone]))
]
function _prior_posterior_row(label, post, prior)
    f(x) = string(round(x; sigdigits = 3))
    ci(v) = string(
        f(median(v)), " (", f(quantile(v, 0.05)), "–",
        f(quantile(v, 0.95)), ")"
    )
    return (
        parameter = label, posterior = ci(post), prior = ci(prior),
        sd_ratio = round(std(post) / std(prior); digits = 2),
    )
end
zone_prior_table = DataFrame(
    vcat(
        [
            _prior_posterior_row(
                label, _hyper_draws(chn_local, k),
                _hyper_draws(prior_chn_zone, k)
            )
                for (k, label, _) in zone_hyper
        ],
        [
            _prior_posterior_row(
                "Drift scale σ_δ ($(zone_inputs.patch_labels[p]))",
                _drift_draws(chn_local, p), _drift_draws(prior_chn_zone, p)
            )
                for p in eachindex(zone_inputs.patch_labels)
        ]
    )
);
# The pair plot's axes take the symbols alone.
zone_hyper_pair_fig = plot_pair(
    chn_local, zone_hyper_keys;
    prior = prior_chn_zone,
    labels = Dict(k => sym for (k, _, sym) in zone_hyper)
);
parameterposteriorpriorsd_ratio
Level spread σ_level0.575 (0.369–0.826)0.216 (0.0189–0.59)0.79
Deviation half-life (days)60.9 (36.6–114.0)42.2 (14.6–106.0)0.88
Correlation ρ_corr0.00518 (0.000286–0.0383)0.347 (0.0358–0.876)0.06
Ascertainment spread σ_ascertainment0.257 (0.0692–0.381)0.126 (0.0252–0.669)0.3
Severity spread σ_severity0.0988 (0.0092–0.256)0.069 (0.00693–0.182)1.39
Within-patch mixing ε_within0.0188 (0.00859–0.035)0.0326 (0.00232–0.136)0.19
Mixing departure τ_mix0.48 (0.0367–1.31)0.366 (0.0309–1.01)1.26
Case composition ρ0.0308 (0.0271–0.0345)0.0329 (0.023–0.0457)0.33
Death composition ρ_death0.0528 (0.0457–0.0607)0.067 (0.0493–0.0873)0.39
Drift scale σ_δ (Ituri)0.244 (0.182–0.319)0.118 (0.067–0.211)0.93
Drift scale σ_δ (Nord-Kivu)0.229 (0.162–0.338)0.122 (0.071–0.218)1.2
Drift scale σ_δ (Haut-Uele)0.209 (0.133–0.318)0.122 (0.0684–0.215)1.24
Drift scale σ_δ (Other provinces)0.122 (0.0692–0.22)0.12 (0.0663–0.221)0.97
Zone hyperparameter pair plot (prior overlaid)

Change over the past week ​

The health-zone model is melded onto the headline fit's patch infections one way, so the zone data do not update the national and province estimates. The comparison below reads the reproduction number of every zone walking in both fits at the frozen and live cut-offs, matched by key. The dot plot shows the fifteen zones the frozen fit ranks highest, the trajectories the twelve with most confirmed cases, and the table the ten with most confirmed cases.

Zone cut-off summaries shared by the comparisons
julia
# One row per zone in `zs` with the median and the 50% and 90% intervals
# of its draws (one vector per zone), in the schema the zone dot plots
# read, plus the zone's key to match variants by and its confirmed cases
# to rank by. A zone with no finite draw is left out.
function _zone_cutoff_summary(draws, inputs, zs = eachindex(draws))
    keep = [i for (i, v) in enumerate(draws) if any(isfinite, v)]
    z = zs[keep]
    tbl = zone_summary_table(
        [filter(isfinite, draws[i]) for i in keep],
        inputs.zone_labels[z], inputs.patch_of_zone[z]
    )
    tbl.key = inputs.zone_keys[z]
    tbl.cases = inputs.cumulative[z]
    return tbl
end

# The `top` zones by confirmed cases in the first variant, one column per
# variant for the reproduction number. Each cell is a median with its 90%
# interval. Variants are matched by zone key.
function _zone_comparison_table(variants; top::Integer = 10)
    function cell(t, key, d, scale)
        i = findfirst(==(key), t.key)
        i === nothing && return ""
        f(x) = string(round(scale * x; digits = d))
        return string(f(t.median[i]), " (", f(t.lo90[i]), "–", f(t.hi90[i]), ")")
    end
    base = first(last(first(variants)))
    order = sortperm(base.cases; rev = true)[1:min(top, size(base, 1))]
    df = DataFrame(
        "Zone" => base.label[order],
        "Province" => [PROVINCE_LABELS[p] for p in base.patch[order]],
        "Cases" => base.cases[order]
    )
    for (q, name, d, scale) in ((:R, "R", 2, 1),),
            (label, s) in variants

        haskey(s, q) || continue
        df[!, "$name ($label)"] = [
            cell(s[q], k, d, scale)
                for k in base.key[order]
        ]
    end
    return df
end
_zone_comparison_table (generic function with 1 method)

The frozen zone fit and the live one condition on different weeks of data and on different parent fits. The comparison reads each zone's reproduction number on the frozen cut-off day from both fits, alongside the live estimate at the current cut-off. A zone carries its own walk only once it has reported 30 confirmed cases, and the comparison is restricted to zones walking in both fits. The trajectory panels draw the frozen fit behind the live one for the twelve such zones with most confirmed cases.

Zone reproduction numbers from the frozen and live fits
julia
# Daily reproduction numbers rebuilt from both fits, and the zones walking
# in both as live index => frozen index, matched by key. The summaries
# take their labels and cases from the live inputs.
zone_rt_live = reconstruct_zone_rt(chn_local, zone_inputs)
zone_rt_frozen = reconstruct_zone_rt(frozen_local.chn, frozen_zone_inputs)
zone_week_pairs = [
    z => j
        for (z, k) in enumerate(zone_inputs.zone_keys)
        for j in (findfirst(==(k), frozen_zone_inputs.zone_keys),)
        if j !== nothing && zone_inputs.walking[z] &&
        frozen_zone_inputs.walking[j]
]
zone_week_variants = let zs = first.(zone_week_pairs), js = last.(zone_week_pairs),
        n_f = frozen_zone_inputs.n

    [
        "frozen fit at its cut-off" => (;
            R = _zone_cutoff_summary(
                [zone_rt_frozen[j][:, n_f] for j in js],
                zone_inputs, zs
            ),
        ),
        "live fit on the same day" => (;
            R = _zone_cutoff_summary(
                [zone_rt_live[z][:, n_f] for z in zs],
                zone_inputs, zs
            ),
        ),
        "live fit at its cut-off" => (;
            R = _zone_cutoff_summary(
                [zone_rt_live[z][:, end] for z in zs],
                zone_inputs, zs
            ),
        ),
    ]
end
zone_week_table = _zone_comparison_table(zone_week_variants);
zone_week_fig = plot_zone_comparison(
    [l => s.R for (l, s) in zone_week_variants];
    xlabel = "Reproduction number", reference_line = 1.0,
    title = "Zone reproduction number from the frozen and live fits"
);

# The frozen trajectories padded onto the live grid, undefined past the
# frozen cut-off, behind the live ones over the zone grid. The plot reads
# an undefined day as `missing`, where the reconstruction writes `NaN`.
zone_week_rt_fig = let grid = zone_inputs.t0:obs.n, n_f = frozen_zone_inputs.n
    asmissing(m) = replace(m, NaN => missing)
    frozen = map(zone_week_pairs) do (z, j)
        m = fill(NaN, size(zone_rt_frozen[j], 1), obs.n)
        m[:, 1:n_f] .= zone_rt_frozen[j]
        asmissing(m[:, grid])
    end
    zs = first.(zone_week_pairs)
    plot_rt_zones(
        [asmissing(zone_rt_live[z][:, grid]) for z in zs],
        zone_inputs.zone_labels[zs], zone_inputs.patch_of_zone[zs];
        patch_labels = zone_inputs.patch_labels,
        dates = grid_date.(grid), as_of_date = obs.cutoff,
        cumulative = zone_inputs.cumulative[zs], top = 12,
        reference_rt = frozen, reference_label = "Frozen fit",
        title = "Zone reproduction number from the live fit, " *
            "with the frozen fit behind"
    )
end

ZoneProvinceCasesR (frozen fit at its cut-off)R (live fit on the same day)R (live fit at its cut-off)
BuniaIturi16970.79 (0.49–1.18)0.54 (0.33–0.81)0.44 (0.2–0.8)
RwamparaIturi10780.71 (0.43–1.13)0.6 (0.36–0.93)0.5 (0.24–0.9)
NiziIturi7710.64 (0.38–1.04)0.46 (0.27–0.72)0.38 (0.17–0.74)
MongbwaluIturi6971.12 (0.61–1.81)0.69 (0.41–1.09)0.55 (0.26–1.04)
KatwaNord-Kivu5550.71 (0.44–1.1)0.46 (0.28–0.72)0.48 (0.26–0.84)
MangalaIturi3460.71 (0.39–1.18)0.49 (0.28–0.81)0.4 (0.18–0.79)
BeniNord-Kivu3291.2 (0.79–1.71)1.13 (0.79–1.56)1.03 (0.69–1.44)
Nia-NiaIturi2841.03 (0.5–1.94)0.96 (0.55–1.56)0.77 (0.34–1.41)
LitaIturi2830.82 (0.49–1.3)0.58 (0.34–0.92)0.46 (0.22–0.89)
ButemboNord-Kivu2560.8 (0.49–1.22)0.5 (0.3–0.81)0.51 (0.27–0.91)

Health-zone map ​

Each affected health zone coloured by its current reproduction number, with the chance that number exceeds one, the seven-day confirmed-case forecast, the chance of at least a chosen number of cases, the confirmed cases to date and the zone's share of its patch's infections available from the switcher. A filter shows the zones with or without a case over the past one, two or four weeks, and the reproduction number and the forecast can be read at their median or at either bound of the 90% interval. A zone whose reproduction number is from its province, not modelled separately, is hatched, and a filter shows either kind alone. A reproduction number whose 90% interval spans one is paler, and province outlines are drawn over the zones. Hover over a zone for its estimate and 90% credible interval, click it for every number, or open the table view to sort by any column. The map header gives the data cut-off and a link to download the estimates as a CSV file. The map needs a browser. It appears only on the documentation site.

Open the map in a new tab.

Saving zone assets ​

The interactive map above reads the per-zone estimates written here. The zone forecast figure is written by the health-zone forecasts page and the frozen zone forecast and its scores by the health-zone forecast evaluation page.

Write the zone dashboard and release assets
julia
dashboard_dir = joinpath(
    pkgdir(BVDOutbreakSize), "docs", "src", "summary_assets"
)
mkpath(dashboard_dir)
# The health-zone maps and the one-week zone forecast for the summary
# dashboard, and the per-zone estimates the interactive map reads: one row per zone keyed as the
# geojson keys it, with the cases and deaths to date, the reproduction
# number and the chance it exceeds one, the one-week forecast and the
# share of the patch's infections, whether the zone's reproduction number
# is modelled separately (`walking`) and the data cut-off. A zone below the
# reporting floor carries no reproduction number.
CairoMakie.save(joinpath(dashboard_dir, "zone_rt_map.png"), zone_map_fig)
_zone_deaths = [
    let h = obs.zone_death_history
        haskey(h, prov) && haskey(h[prov], z) &&
            !isempty(h[prov][z].counts) ? Int(h[prov][z].counts[end]) :
            0
    end
        for (prov, z) in zip(
            zone_inputs.zone_province,
            zone_inputs.zone_names
        )
]
_zone_rt_stat(f) = [isempty(r) ? missing : f(r) for r in zone_RT_finite]
zone_estimates = DataFrame(
    zone = zone_map_keys,
    label = zone_inputs.zone_labels, province = zone_inputs.zone_province,
    patch = zone_inputs.patch_labels[zone_patch],
    cases = zone_inputs.cumulative, deaths = _zone_deaths,
    R_T_median = _zone_rt_stat(median),
    p_rt_above_one = [
        isempty(r) ? missing : mean(r .> 1) for r in zone_RT_finite
    ],
    R_T_lower = _zone_rt_stat(r -> quantile(r, 0.05)),
    R_T_upper = _zone_rt_stat(r -> quantile(r, 0.95)),
    forecast_median = [median(v) for v in zone_fc_draws.zones],
    forecast_lower = [quantile(v, 0.05) for v in zone_fc_draws.zones],
    forecast_upper = [quantile(v, 0.95) for v in zone_fc_draws.zones],
    share_median = [median(v) for v in zone_share_T],
    share_lower = [quantile(v, 0.05) for v in zone_share_T],
    share_upper = [quantile(v, 0.95) for v in zone_share_T],
    p_ge_1 = zone_fc_probs[:, 1], p_ge_5 = zone_fc_probs[:, 2],
    p_ge_10 = zone_fc_probs[:, 3], p_ge_20 = zone_fc_probs[:, 4],
    cases_last_7 = zone_recent[7], cases_last_14 = zone_recent[14],
    cases_last_28 = zone_recent[28],
    last_case_date = [ismissing(d) ? "" : string(d) for d in zone_last_case],
    walking = Int.(zone_inputs.walking),
    as_of = fill(string(obs.cutoff), length(zone_map_keys))
)
CSV.write(joinpath(dashboard_dir, "zone_estimates.csv"), zone_estimates)
# The same frame into the release outputs.
output_dir = get(
    ENV, "BVD_OUTPUT_DIR",
    joinpath(pkgdir(BVDOutbreakSize), "output")
)
mkpath(output_dir)
CSV.write(joinpath(output_dir, "zone_estimates.csv"), zone_estimates)
"/home/runner/work/BVDOutbreakSize/BVDOutbreakSize/output/zone_estimates.csv"

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