Health-zone in-sample checks
This page checks whether the health-zone model reproduces the per-zone data it was fitted to. The national checks are on the national in-sample checks page and the provincial ones on the province in-sample checks page. The zone data enter the model as compositions of their patch's totals, so the checks are on each zone's share of its patch. The last section sets each zone's reproduction number from the current fit against the same fit a week earlier.
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
# carries the patch trajectories the zone stage conditions on, and the
# frozen zone fit is the same fit a week earlier.
chn_joint = load_fit("joint");
chn_local = load_fit("local");
frozen_local = load_fit("local_frozen_validation");Summary
The overall bullets come first, then a block per province, from the calibration below. Bias and coverage are defined as on the national in-sample checks page, with a nominal coverage of 0.9.
Least well reproduced: Haut-Uele (cases) (bias 0.03, 90% coverage 0.95); Nord-Kivu (cases) (bias -0.01, 90% coverage 0.98); Tshopo (cases) (bias 0.01, 90% coverage 0.97).
Coverage: 16 of 16 composition rows have 90% coverage of at least 0.8; the lowest is Haut-Uele (cases) at 0.95.
Ituri
Confirmed cases: bias -0.0 and 90% coverage 0.97 over 3388 cells.
Confirmed deaths: bias 0.0 and 90% coverage 0.98 over 3304 cells.
Nord-Kivu
Confirmed cases: bias -0.01 and 90% coverage 0.98 over 1972 cells.
Confirmed deaths: bias -0.0 and 90% coverage 0.98 over 1887 cells.
Haut-Uele
Confirmed cases: bias 0.03 and 90% coverage 0.95 over 567 cells.
Confirmed deaths: bias 0.01 and 90% coverage 0.98 over 469 cells.
Sud-Kivu
Confirmed cases: bias 0.0 and 90% coverage 1.0 over 1 cells.
Confirmed deaths: bias 0.0 and 90% coverage 1.0 over 1 cells.
Tshopo
Confirmed cases: bias 0.01 and 90% coverage 0.97 over 266 cells.
Confirmed deaths: bias 0.0 and 90% coverage 0.99 over 168 cells.
Bas-Uele
Confirmed cases: bias 0.01 and 90% coverage 1.0 over 30 cells.
Confirmed deaths: bias -0.0 and 90% coverage 1.0 over 12 cells.
Sud-Ubangi
Confirmed cases: bias 0.0 and 90% coverage 1.0 over 2 cells.
Confirmed deaths: bias 0.0 and 90% coverage 1.0 over 1 cells.
Nord-Ubangi
Confirmed cases: bias 0.0 and 90% coverage 1.0 over 2 cells.
Confirmed deaths: bias 0.0 and 90% coverage 1.0 over 1 cells.
Zone compositions
The zone compositions sit one level below the province compositions. The cases run on the infection-to-report delay and the deaths on the infection-to-confirmed-death delay, each with its own overdispersion, so the two are checked separately.
Load the zone fit's inputs and prior
julia
# `zone_stage_inputs` and `zone_prior_draws` are defined in the shared setup,
# so the zone estimates and forecast pages read the same fixed inputs.
zone_inputs = zone_stage_inputs();
zone_patch = zone_inputs.patch_of_zone;
prior_chn_zone = zone_prior_draws(zone_inputs);Confirmed cases
Zone case composition predictive checks
julia
# Per-draw expected, predictive and observed zone shares at every vintage,
# from the posterior and from the zone prior, per vintage and cumulative.
zone_ppc = zone_composition_draws(chn_local, zone_inputs);
zone_ppc_prior = zone_composition_draws(prior_chn_zone, zone_inputs);
zone_ppc_cum = zone_composition_draws(
chn_local, zone_inputs;
cumulative = true
);
zone_ppc_prior_cum = zone_composition_draws(
prior_chn_zone, zone_inputs;
cumulative = true
);
zone_ppc_fig = plot_zone_shares(
zone_ppc.expected, zone_ppc.observed,
zone_inputs.dates, zone_inputs.zone_labels;
zone_patch, patch_labels = zone_inputs.patch_labels,
pred_draws = zone_ppc.predictive_share,
prior_draws = zone_ppc_prior.predictive_share, top = 15, ncols = 5
);
zone_ppc_cum_fig = plot_zone_shares(
zone_ppc_cum.expected,
zone_ppc_cum.observed, zone_inputs.dates, zone_inputs.zone_labels;
zone_patch, patch_labels = zone_inputs.patch_labels,
pred_draws = zone_ppc_cum.predictive_share,
prior_draws = zone_ppc_prior_cum.predictive_share, top = 15, ncols = 5,
title = "Zone share of the patch's cumulative confirmed cases, " *
"modelled against observed"
);
Confirmed deaths
Zone death composition predictive check
julia
zone_death_ppc_cum = zone_composition_draws(
chn_local, zone_inputs;
stream = :deaths, cumulative = true
);
zone_death_ppc_prior_cum = zone_composition_draws(
prior_chn_zone, zone_inputs;
stream = :deaths, cumulative = true
);
zone_death_ppc_cum_fig = plot_zone_shares(
zone_death_ppc_cum.expected, zone_death_ppc_cum.observed,
zone_inputs.dates, zone_inputs.zone_labels;
zone_patch, patch_labels = zone_inputs.patch_labels,
pred_draws = zone_death_ppc_cum.predictive_share,
prior_draws = zone_death_ppc_prior_cum.predictive_share,
top = 15, ncols = 5,
title = "Zone share of the patch's cumulative confirmed deaths, " *
"modelled against observed"
);
Zone composition calibration
Both compositions are scored the way the stream calibration scores the national streams, over every zone and vintage of each patch.
Zone composition calibration
julia
# The two compositions in one table, each patch row named by the
# composition it belongs to so the calibration plot reads them apart.
zone_calibration_table = let
function tagged(stream, what)
tbl = zone_composition_calibration(
chn_local, zone_inputs; stream = stream
)
tbl[!, "Stream"] = string.(tbl[!, "Stream"], " ($what)")
return tbl
end
vcat(tagged(:cases, "cases"), tagged(:deaths, "deaths"))
end;
zone_calibration_fig = plot_stream_calibration(zone_calibration_table);
| Stream | Vintages | Bias | 50% coverage | 90% coverage |
|---|---|---|---|---|
| Ituri (cases) | 3388 | 0 | 0.83 | 0.97 |
| Nord-Kivu (cases) | 1972 | -0.01 | 0.87 | 0.98 |
| Haut-Uele (cases) | 567 | 0.03 | 0.81 | 0.95 |
| Sud-Kivu (cases) | 1 | 0 | 1 | 1 |
| Tshopo (cases) | 266 | 0.01 | 0.89 | 0.97 |
| Bas-Uele (cases) | 30 | 0.01 | 0.9 | 1 |
| Sud-Ubangi (cases) | 2 | 0 | 1 | 1 |
| Nord-Ubangi (cases) | 2 | 0 | 1 | 1 |
| All zones (cases) | 6228 | 0 | 0.85 | 0.97 |
| Ituri (deaths) | 3304 | 0 | 0.9 | 0.98 |
| Nord-Kivu (deaths) | 1887 | 0 | 0.89 | 0.98 |
| Haut-Uele (deaths) | 469 | 0.01 | 0.9 | 0.98 |
| Sud-Kivu (deaths) | 1 | 0 | 1 | 1 |
| Tshopo (deaths) | 168 | 0 | 0.91 | 0.99 |
| Bas-Uele (deaths) | 12 | 0 | 0.92 | 1 |
| Sud-Ubangi (deaths) | 1 | 0 | 1 | 1 |
| Nord-Ubangi (deaths) | 1 | 0 | 1 | 1 |
| All zones (deaths) | 5843 | 0 | 0.9 | 0.98 |
Change over the past week
The comparison reads each zone's reproduction number on the frozen cut-off day from the frozen and live fits, alongside the live estimate at the current cut-off. The two fits condition on different weeks of data and on different parent fits. A zone carries its own walk only once it has reported 30 confirmed cases, so the comparison is restricted to zones walking in both fits.
Zone cut-off summaries shared by the comparisons
julia
# The frozen fit's fixed inputs, from the shared setup.
frozen_zone_inputs = frozen_zone_stage_inputs();
# 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)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,
ranking = zone_inputs.cumulative[zs],
ranking_label = "confirmed cases to date", top = 12,
reference_rt = frozen, reference_label = "Frozen fit",
title = "Zone reproduction number from the live fit, " *
"with the frozen fit behind"
)
end
| Zone | Province | Cases | R (frozen fit at its cut-off) | R (live fit on the same day) | R (live fit at its cut-off) |
|---|---|---|---|---|---|
| Bunia | Ituri | 1914 | 0.85 (0.33–1.65) | 0.99 (0.54–1.51) | 0.97 (0.4–1.72) |
| Rwampara | Ituri | 1173 | 0.8 (0.32–1.65) | 1.0 (0.57–1.62) | 1.01 (0.41–1.92) |
| Nizi | Ituri | 832 | 0.88 (0.35–1.87) | 0.85 (0.43–1.48) | 0.87 (0.31–1.84) |
| Mongbwalu | Ituri | 717 | 0.63 (0.22–1.79) | 0.71 (0.31–1.68) | 0.8 (0.27–2.12) |
| Katwa | Nord-Kivu | 648 | 0.92 (0.51–1.6) | 0.86 (0.54–1.27) | 1.01 (0.57–1.68) |
| Beni | Nord-Kivu | 553 | 1.07 (0.67–1.59) | 1.12 (0.84–1.47) | 1.11 (0.75–1.55) |
| Mangala | Ituri | 382 | 0.74 (0.26–1.95) | 0.76 (0.35–1.49) | 0.83 (0.32–2.06) |
| Mandima | Ituri | 362 | 0.72 (0.24–1.67) | 1.11 (0.59–1.73) | 0.98 (0.34–1.79) |
| Nia-Nia | Ituri | 334 | 0.73 (0.21–2.18) | 0.65 (0.28–1.37) | 0.72 (0.26–1.97) |
| Butembo | Nord-Kivu | 319 | 0.91 (0.49–1.63) | 0.89 (0.55–1.33) | 1.04 (0.58–1.7) |
Saving zone in-sample outputs
Write the summary bullets
julia
# The bullets under the summary heading at the top of the page. They read
# the calibration table built further down, so they are written here and
# read back when the site is assembled.
evaluation_insample_zone_summary = let
fmt(x) = ismissing(x) || !isfinite(x) ? "n/a" :
string(round(x; digits = 2))
# The per-patch rows only: the table's "All zones" rows pool them and
# would be counted twice.
fitted = filter(
r -> isfinite(r["Bias"]) && !startswith(r["Stream"], "All zones"),
zone_calibration_table
)
worst = first(
sort(fitted, "Bias"; by = abs, rev = true), min(3, size(fitted, 1))
)
low = first(sort(fitted, "90% coverage"), min(1, size(fitted, 1)))
overall = [
string(
"- **Least well reproduced:** ",
join(
[
string(
r["Stream"], " (bias ", fmt(r["Bias"]),
", 90% coverage ", fmt(r["90% coverage"]), ")"
)
for r in eachrow(worst)
], "; "
), "."
),
string(
"- **Coverage:** ", count(>=(0.8), fitted[!, "90% coverage"]),
" of ", size(fitted, 1),
" composition rows have 90% coverage of at least 0.8",
size(low, 1) == 0 ? "." :
string(
"; the lowest is ", low[1, "Stream"], " at ",
fmt(low[1, "90% coverage"]), "."
)
),
]
rows(tbl) = Dict(r["Stream"] => r for r in eachrow(tbl))
cal = rows(zone_calibration_table)
function calibration(p, what, label)
r = get(cal, string(PROVINCE_LABELS[p], " (", what, ")"), nothing)
r === nothing && return string("- ", label, ": not scored.")
return string(
"- ", label, ": bias ", fmt(r["Bias"]), " and 90% coverage ",
fmt(r["90% coverage"]), " over ", r["Vintages"], " cells."
)
end
detail = [
join(
[
string("**", PROVINCE_LABELS[p], "**"), "",
calibration(p, "cases", "Confirmed cases"),
calibration(p, "deaths", "Confirmed deaths"),
], "\n"
)
for p in 1:N_PATCHES
]
join(vcat([join(overall, "\n")], detail), "\n\n")
end
dashboard_dir = joinpath(
pkgdir(BVDOutbreakSize), "docs", "src", "summary_assets"
)
mkpath(dashboard_dir)
write(
joinpath(dashboard_dir, "evaluation_insample_zone.md"),
evaluation_insample_zone_summary
);