In-sample checks

Whether the fitted joint model reproduces the national data it was fitted to. It also sets the single-stream fits against the joint, breaks the sampler diagnostics down by parameter and re-fits the joint under alternative assumptions. The same checks by province are on the province in-sample checks page and by health zone on the health-zone in-sample checks page. How the model predicts data it has not seen is on the forecast evaluation page.

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 and prior draws this page reads, loaded from the cache here.
chn_joint = load_fit("joint")
prior_chn = joint_prior_draws()
# `sens_no_patches` is the headline with `n_patches = 1`, the check on the
# spatial structure.
chn_no_patches = load_fit("sens_no_patches")
chn_exports = load_fit("exports")
chn_deaths = load_fit("deaths")
chn_cases = load_fit("cases")
chn_confirmed = load_fit("confirmed")
chn_confirmed_deaths = load_fit("confirmed_deaths")
chn_treatment = load_fit("treatment")
chn_onsets = load_fit("onsets")
frozen_lastweek = load_fit("frozen_validation")
if RUN_SENSITIVITY
    chn_joint_community_delay = load_fit("sens_community_delay")
    chn_joint_exp_growth_clock = load_fit("sens_exp_growth_clock")
end
posterior_C_joint = vec(Array(chn_joint[:C_T]))
posterior_C_exports = vec(Array(chn_exports[:C_T]))
posterior_C_deaths = vec(Array(chn_deaths[:C_T]))
posterior_C_cases = vec(Array(chn_cases[:C_T]))
posterior_C_confirmed = vec(Array(chn_confirmed[:C_T]))
posterior_C_treatment = vec(Array(chn_treatment[:C_T]))
posterior_C_onsets = vec(Array(chn_onsets[:C_T]));

Summary

Whether each stream is reproduced, from the checks further down this page. Bias runs from −1 to 1 and is negative when the model under-predicts, and coverage is the fraction of vintages inside the predictive interval.

Streams: 11 of 17 fitted streams have 90% coverage of at least 0.8.
Least well reproduced: Specimens analysed (cumulative) (bias -0.14, 90% coverage 0.5); In-care deaths/day (bias 0.14, 90% coverage 0.81); New suspected deaths/day (bias -0.11, 90% coverage 0.94).
Exports: Uganda exports observed 3 against a predictive median of 1 (90% interval 0–3), and export deaths observed 1 against a predictive median of 0 (90% interval 0–0).

Streams still reporting

New suspects/day: bias 0.03, 90% coverage 0.96 over 109 vintages.
Patients in isolation: bias 0.01, 90% coverage 0.94 over 70 vintages.
Confirmed cases: bias 0.07, 90% coverage 0.91 over 127 vintages.
Confirmed deaths: bias 0.04, 90% coverage 0.89 over 129 vintages.
Recovered (confirmed): bias 0.03, 90% coverage 0.92 over 108 vintages.
Specimens analysed (24h): bias -0.05, 90% coverage 0.97 over 110 vintages.
Onset reports (net correction/snapshot): bias -0.1, 90% coverage 0.62 over 53 vintages.

Streams no longer reporting

Suspected cases: bias -0.02, 90% coverage 0.67 over 9 vintages.
Suspected deaths: bias -0.03, 90% coverage 0.78 over 9 vintages.
New suspected deaths/day: bias -0.11, 90% coverage 0.94 over 32 vintages.
Specimens analysed (cumulative): bias -0.14, 90% coverage 0.5 over 6 vintages.
Admissions/day: bias 0.05, 90% coverage 0.76 over 42 vintages.
In-care deaths/day: bias 0.14, 90% coverage 0.81 over 42 vintages.
Rule-outs/day: bias -0.0, 90% coverage 0.56 over 41 vintages.
Absconded/day: bias 0.09, 90% coverage 0.86 over 42 vintages.
Confirmed in care: bias -0.04, 90% coverage 1.0 over 38 vintages.
Suspects in care: bias -0.07, 90% coverage 0.92 over 38 vintages.

Prior predictive check

Whether the prior, before any data are fitted, brackets the observed counts.

Summarise the joint prior

julia

prior_C_table = summary_table(prior_chn, [:C_T]; digits = 0);

Show prior summary table

1×7 DataFrame

Row	Quantity	Lower 90%	Lower 60%	Lower 30%	Upper 30%	Upper 60%	Upper 90%
	String	Float64	Float64	Float64	Float64	Float64	Float64
1	C_T	1999.0	17714.0	92431.0	2.9923e6	1.53113e7	2.59195e7

Prior pair plot

julia

prior_pair_fig = plot_pair(
    prior_chn,
    [
        :C_T, :R_T, :r, :T, :CFR, :k,
        :p_drc, :p_uganda,
    ]
);

Posterior predictive checks

Whether data replicated from the fitted model reproduce each observed stream.

Streams still reporting

Whether the streams reported within a week of the cut-off are reproduced.

Cumulative

Whether the replicated running totals, or daily counts for a daily stream, track the observed ones.

Joint posterior predictive plot

julia

# The joint posterior predictive, shared with the province page (see
# `joint_posterior_predictive` in `docs/pages/_setup.jl`).
pp_joint = joint_posterior_predictive();

# `predict` stores each stream's per-vintage increments as one
# vector-valued variable (`<stream>_increments.increments`); the slice is
# an iter×chain matrix of per-draw increment vectors, exactly the
# `replicates` shape `plot_vintage_conditional_ppc` grounds on each
# vintage's observed previous cumulative for the one-step-ahead
# predictive. Look it up by its VarName with FlexiChains' `Prefixed`, which
# matches a (submodel-prefixed) key by its varname tail: `Prefixed(@varname(
# reported_increments.increments))` finds `cases_state.reported_increments.
# increments` without hard-coding the `cases_state.` prefix, and matches by
# the varname tail rather than a loose substring, so it cannot be fooled by a
# scalar `expected_*_T` deterministic. `FlexiChains` is a package
# dependency (imported, not exported), so it is reached through the package
# namespace.
const _Prefixed = BVDOutbreakSize.FlexiChains.Prefixed;
_vintage_replicates(pp, vn) = collect(pp[_Prefixed(vn)]);

# Grid day-index → INSP situation-report date label.
_vintage_dates(days) = string.(obs.seeding .+ Day.(days .- 1));

reported_panel = (;
    id = :suspected_cases,
    title = "Suspected cases",
    dates = _vintage_dates(obs.reported_history.days),
    replicates = _vintage_replicates(
        pp_joint, @varname(reported_increments.increments)
    ),
    observed = obs.reported_history.counts, colour = :steelblue,
);
# Daily new-suspect inflow: a per-day count (not cumulative), so the panel
# is drawn with `cumulative = false` — each replicate is its own daily
# count against the observed daily count rather than a running total. Its
# days pick up where the cumulative suspected panel freezes.
suspected_daily_panel = (;
    id = :suspected_daily,
    title = "New suspects/day",
    dates = _vintage_dates(obs.suspected_daily_history.days),
    replicates = _vintage_replicates(
        pp_joint, @varname(suspected_daily.increments)
    ),
    observed = obs.suspected_daily_history.counts,
    colour = :slateblue, cumulative = false,
);
# Isolation/treatment-bed occupancy: a census stock, so the panel is drawn
# with `cumulative = false` — each replicate is the modelled bed count on a
# report day against the observed "Patients en isolement" count. It is a
# level, not a count of new events, so it carries its own `ylabel` rather
# than the "Daily count" and "New per vintage" defaults, which would read as
# an accumulating total on a series that rises through the outbreak. The count
# is the suspect inflow carried through a length-of-stay survival, so its
# level and lag reflect the admission proportion and the stays. The censored-
# occupancy likelihood stores its per-day predictive draws under the submodel
# `obs` variable (not `increments`), so the replicates are read from that key.
# The treatment model scores the total occupancy only on the days without a
# published confirmed/suspect split (a per-day total-or-split switch): on the
# split days the two sub-stock census panels carry the fit instead, so the
# `isolation.obs` predictive holds only the non-split days. Drop the split
# days from the panel's dates and observed counts to match that length.
_iso_split_days = Set(Int.(obs.treatment_confirmed_incare_history.days))
_iso_keep = [!(Int(d) in _iso_split_days) for d in obs.isolation_history.days]
isolation_panel = (;
    id = :isolation_beds,
    title = "Patients in isolation",
    dates = _vintage_dates(obs.isolation_history.days[_iso_keep]),
    replicates = _vintage_replicates(
        pp_joint, @varname(isolation.obs)
    ),
    observed = obs.isolation_history.counts[_iso_keep],
    colour = :darkorange, cumulative = false,
    ylabel = "Beds occupied",
);
deaths_panel = (;
    id = :suspected_deaths,
    title = "Suspected deaths",
    dates = _vintage_dates(obs.deaths_history.days),
    replicates = _vintage_replicates(
        pp_joint, @varname(death_increments.increments)
    ),
    observed = obs.deaths_history.counts, colour = :firebrick,
);
# Daily new suspected deaths: a per-day count (not cumulative), so the panel
# is drawn with `cumulative = false` — each replicate is its own daily count
# against the observed daily count rather than a running total. Its days
# pick up where the cumulative suspected-death panel freezes, the deaths
# analogue of the new-suspects-per-day panel.
suspected_daily_deaths_panel = (;
    id = :suspected_daily_deaths,
    title = "New suspected deaths/day",
    dates = _vintage_dates(obs.suspected_daily_deaths_history.days),
    replicates = _vintage_replicates(
        pp_joint, @varname(suspected_daily_deaths.increments)
    ),
    observed = obs.suspected_daily_deaths_history.counts,
    colour = :indianred, cumulative = false,
);
# Specimens analysed is the single modelled laboratory volume (the
# report-to-analysed delay and tested-fraction throughput), fit to the
# cumulative analysed series, so it gets the same cumulative conditional
# check as the suspected streams. This is the testing volume the
# confirmed-positivity denominator is built from.
tests_analysed_panel = (;
    id = :tests_analysed,
    title = "Specimens analysed (cumulative)",
    dates = _vintage_dates(obs.lab_history.days),
    replicates = _vintage_replicates(
        pp_joint, @varname(analysed_increments.increments)
    ),
    observed = obs.lab_history.counts, colour = :seagreen,
);
# Post-cutoff 24h analysed volume: once the cumulative series stops, INSP
# reports a 24h analysed count on some days. These are fitted as per-day
# volumes (not cumulative), so the panel is a standalone daily check
# (`cumulative = false`): the modelled daily analysed volume against the
# observed 24h count on each reported day.
tests_analysed_daily_panel = (;
    id = :tests_analysed_daily,
    title = "Specimens analysed (24h)",
    dates = _vintage_dates(obs.lab_daily_history.days),
    replicates = _vintage_replicates(
        pp_joint, @varname(analysed_daily_increments.increments)
    ),
    observed = obs.lab_daily_history.counts, colour = :teal,
    cumulative = false,
);

# Confirmed cases are scored over two groups of laboratory windows: the
# early confirmed vintages (no per-vintage analysed denominator, scored
# as counts against the modelled laboratory volume) and the observed
# windows (a Binomial of the observed analysed denominator). Both groups
# produce per-window replicate increments in `predict`, so concatenating
# them oldest-first gives the per-vintage cumulative confirmed-case
# trajectory, grounded on the observed cumulative confirmed at each window
# end-day. The 24-25 May analysis stall merges into 26 May, so the window
# grid is slightly coarser than the raw confirmed history.
_conf_windows = BVDOutbreakSize.confirmed_positivity_windows(
    obs.confirmed_history, obs.lab_history, obs.lab_daily_history
);
# Oldest-first: early (no denominator) → observed (analysed Binomial) →
# late (post-28 May; trusted 24h-analysed days are Binomial windows, the
# rest unanchored windows scored against the modelled volume).
_conf_window_days = vcat(
    _conf_windows.early_days, _conf_windows.obs_days,
    _conf_windows.late_days
);
function _confirmed_at(day)
    i = searchsortedlast(obs.confirmed_history.days, day)
    return i == 0 ? 0 : Int(obs.confirmed_history.counts[i])
end;
_conf_early = _vintage_replicates(
    pp_joint, @varname(early_increments.increments)
);
_conf_obs = collect(
    first(
        pp_joint[k]
            for k in keys(pp_joint)
            if occursin("confirmed_state.confirmed_positives.positives", string(k))
    )
);
_conf_late = _vintage_replicates(
    pp_joint, @varname(late_increments.increments)
);
confirmed_panel = (;
    id = :confirmed_cases,
    title = "Confirmed cases",
    dates = _vintage_dates(_conf_window_days),
    replicates = [
        vcat(collect(e), collect(p), collect(l))
            for (e, p, l) in zip(vec(_conf_early), vec(_conf_obs), vec(_conf_late))
    ],
    observed = [_confirmed_at(d) for d in _conf_window_days],
    colour = :goldenrod,
);

# Confirmed deaths are a per-vintage stream, scored as increments of the
# modelled confirmed-death trajectory up to the cut-off, so they get the
# same cumulative conditional check.
confirmed_deaths_panel = (;
    id = :confirmed_deaths,
    title = "Confirmed deaths",
    dates = _vintage_dates(obs.confirmed_deaths_history.days),
    replicates = _vintage_replicates(
        pp_joint, @varname(cdeath_increments.increments)
    ),
    observed = obs.confirmed_deaths_history.counts, colour = :purple,
);

# Recovered among confirmed ("cumul guéris") is a cumulative per-vintage
# stream fitted through the increments of the modelled recovered trajectory
# (the confirmation-to-recovery convolution of the daily confirmed cases) up
# to the cut-off, so it gets the same cumulative conditional check.
recovered_panel = (;
    id = :recovered,
    title = "Recovered (confirmed)",
    dates = _vintage_dates(obs.recovered_history.days),
    replicates = _vintage_replicates(
        pp_joint, @varname(recovered_increments.increments)
    ),
    observed = obs.recovered_history.counts, colour = :mediumseagreen,
);

# Tableau 6 treatment-centre daily flows (the new patient-movement data
# sources): admissions and the discharge reasons (in-care deaths, rule-outs,
# absconded). Per-day counts, so drawn with `cumulative = false` — each
# replicate is the modelled daily flow on a report day against the observed
# Tableau 6 count.
admissions_panel = (;
    id = :treatment_admissions,
    title = "Admissions/day",
    dates = _vintage_dates(obs.treatment_admissions_history.days),
    replicates = _vintage_replicates(
        pp_joint, @varname(admissions.obs)
    ),
    observed = obs.treatment_admissions_history.counts,
    colour = :teal, cumulative = false,
);
incare_deaths_panel = (;
    id = :treatment_deaths,
    title = "In-care deaths/day",
    dates = _vintage_dates(obs.treatment_deaths_history.days),
    replicates = _vintage_replicates(
        pp_joint, @varname(incare_deaths.increments)
    ),
    observed = obs.treatment_deaths_history.counts,
    colour = :darkred, cumulative = false,
);
ruleouts_panel = (;
    id = :treatment_ruleouts,
    title = "Rule-outs/day",
    dates = _vintage_dates(obs.treatment_ruleout_history.days),
    replicates = _vintage_replicates(
        pp_joint, @varname(ruleouts.increments)
    ),
    observed = obs.treatment_ruleout_history.counts,
    colour = :goldenrod, cumulative = false,
);
absconded_panel = (;
    id = :treatment_absconded,
    title = "Absconded/day",
    dates = _vintage_dates(obs.treatment_absconded_history.days),
    replicates = _vintage_replicates(
        pp_joint, @varname(absconded.increments)
    ),
    observed = obs.treatment_absconded_history.counts,
    colour = :slategray, cumulative = false,
);

# Tableau 6 occupancy split (`dont confirmes` / `dont suspects`): the two
# in-care prevalence sub-stocks. Per-day census counts, so drawn with
# `cumulative = false` — each replicate is the modelled confirmed-in-care or
# suspect-in-care bed count on a report day against the observed sub-stock.
# On these split days the total-occupancy panel is not scored, so the two
# sub-stock panels carry the window instead.
confirmed_incare_panel = (;
    id = :treatment_beds,
    title = "Confirmed in care",
    dates = _vintage_dates(obs.treatment_confirmed_incare_history.days),
    replicates = _vintage_replicates(
        pp_joint, @varname(confirmed_incare_obs.increments)
    ),
    observed = obs.treatment_confirmed_incare_history.counts,
    colour = :darkgoldenrod, cumulative = false,
    ylabel = "Beds occupied",
);
suspect_incare_panel = (;
    id = :suspect_beds,
    title = "Suspects in care",
    dates = _vintage_dates(obs.treatment_suspect_incare_history.days),
    replicates = _vintage_replicates(
        pp_joint, @varname(suspect_incare_obs.increments)
    ),
    observed = obs.treatment_suspect_incare_history.counts,
    colour = :chocolate, cumulative = false,
    ylabel = "Beds occupied",
);

# Symptom-onset reporting triangle: one cell per (onset day, report day)
# pair. The cells sharing a report day are summed into one net correction
# per snapshot to fit the per-vintage panel shape. The onset-by-report grid
# is the snapshot nowcast figure further down. A net correction is not a
# running total, so `cumulative = false`.
_onset_ppc_report_days = sort(unique(obs.onset_curve_history.report_days))
_onset_ppc_groups = [
    findall(==(r), obs.onset_curve_history.report_days)
        for r in _onset_ppc_report_days
]
_onset_ppc_replicates_raw = _vintage_replicates(
    pp_joint, @varname(onset_report_state.increments)
)
onset_panel = (;
    id = :onset_reports,
    title = "Onset reports (net correction/snapshot)",
    dates = _vintage_dates(_onset_ppc_report_days),
    replicates = [
        [sum(collect(rep)[g]) for g in _onset_ppc_groups]
            for rep in vec(_onset_ppc_replicates_raw)
    ],
    observed = [
        sum(obs.onset_curve_history.increments[g])
            for g in _onset_ppc_groups
    ],
    colour = :mediumpurple, cumulative = false,
);

# Each panel runs to its own last vintage, so a stream that keeps
# reporting shows the full series the model is fitting rather than the
# window the streams that stopped earlier cover. The per-stream
# calibration table reads this ordered list too, so it stays whole and
# the two stream groups are filtered out of it.
vintage_panels = [
    reported_panel, suspected_daily_panel, isolation_panel, confirmed_panel,
    deaths_panel, suspected_daily_deaths_panel, confirmed_deaths_panel,
    recovered_panel, tests_analysed_panel, tests_analysed_daily_panel,
    admissions_panel, incare_deaths_panel, ruleouts_panel, absconded_panel,
    confirmed_incare_panel, suspect_incare_panel, onset_panel,
];
# The incidence view drops the treatment-centre flow and occupancy-split
# panels, whose per-day counts are already their own incidence.
vintage_incidence_panels = [
    reported_panel, suspected_daily_panel, isolation_panel, confirmed_panel,
    deaths_panel, suspected_daily_deaths_panel, confirmed_deaths_panel,
    recovered_panel, tests_analysed_panel, tests_analysed_daily_panel,
    onset_panel,
];
# Whether a panel's stream was still being reported at the cut-off, from
# the shared registry rule (last vintage within a week of the cut-off)
# rather than a per-page list of dates that goes stale.
_still_reporting(p) = stream_reporting(obs, p.id);
reporting_panels = filter(_still_reporting, vintage_panels);
stopped_panels = filter(!_still_reporting, vintage_panels);
reporting_incidence_panels = filter(
    _still_reporting, vintage_incidence_panels
);
stopped_incidence_panels = filter(
    !_still_reporting, vintage_incidence_panels
);
joint_vintage_ppc_fig = plot_vintage_conditional_ppc(reporting_panels);

Per-vintage incidence

Whether the count reported between consecutive situation reports is reproduced.

Per-vintage incidence posterior predictive plot

julia

joint_vintage_incidence_fig = plot_vintage_incidence_ppc(
    reporting_incidence_panels
);

Streams no longer reporting

Whether the streams that stopped before the cut-off are reproduced over the dates they cover.

Cumulative

Joint posterior predictive plot

julia

joint_vintage_ppc_stopped_fig = plot_vintage_conditional_ppc(stopped_panels);

Per-vintage incidence

Per-vintage incidence posterior predictive plot

julia

joint_vintage_incidence_stopped_fig = plot_vintage_incidence_ppc(
    stopped_incidence_panels
);

Stream calibration

Whether each stream's per-vintage predictions are calibrated against the observed counts.

julia

stream_calibration_table = stream_calibration(vintage_panels);

Per-stream calibration plot

julia

stream_calibration_fig = plot_stream_calibration(stream_calibration_table);

Per-stream calibration table

17×5 DataFrame

Row	Stream	Vintages	Bias	50% coverage	90% coverage
	String	Int64	Float64	Float64	Float64
1	Suspected cases	9	-0.02	0.33	0.67
2	New suspects/day	109	0.03	0.69	0.96
3	Patients in isolation	70	0.01	0.6	0.94
4	Confirmed cases	127	0.07	0.63	0.91
5	Suspected deaths	9	-0.03	0.56	0.78
6	New suspected deaths/day	32	-0.11	0.59	0.94
7	Confirmed deaths	129	0.04	0.56	0.89
8	Recovered (confirmed)	108	0.03	0.54	0.92
9	Specimens analysed (cumulative)	6	-0.14	0.17	0.5
10	Specimens analysed (24h)	110	-0.05	0.7	0.97
11	Admissions/day	42	0.05	0.45	0.76
12	In-care deaths/day	42	0.14	0.55	0.81
13	Rule-outs/day	41	-0.0	0.27	0.56
14	Absconded/day	42	0.09	0.43	0.86
15	Confirmed in care	38	-0.04	0.71	1.0
16	Suspects in care	38	-0.07	0.68	0.92
17	Onset reports (net correction/snapshot)	53	-0.1	0.25	0.62

Exports

Whether the modelled Uganda export and export-death totals match the observed ones.

Scalar posterior predictive plot

julia

# The dated counts are nested under their submodel prefix as a single
# per-day count vector `<prefix>.counts`; look it up by its VarName with
# `Prefixed` (matching the key by its `<obs>.counts` tail) so the
# deterministic `expected_*_T` quantities cannot be picked up by a loose
# substring, then sum each replicate's per-day vector into the total.
function _dated_total(pp, vn)
    return [sum(v) for v in vec(Array(pp[_Prefixed(vn)]))]
end;

pp_exports = _dated_total(pp_joint, @varname(export_obs.counts));
pp_exports_deaths = _dated_total(
    pp_joint, @varname(death_obs.counts)
);

joint_ppc_fig = plot_posterior_predictive(
    pp_exports, nothing,
    obs.exported_cases, nothing;
    pp_exports_deaths = pp_exports_deaths,
    obs_exports_deaths = obs.exports_deaths
);

Onset snapshot nowcasts

Whether the fitted reporting delay nowcasts each digitised onset snapshot to the latest figure covering its dates (see the symptom-onset reporting delay Methods section).

Nowcasts of the digitised reporting-triangle snapshots

julia

# Each snapshot's own printed counts, read from the source blocks since the
# fitted stream holds only the corrections between snapshots.
_onset_readings = onset_snapshot_readings()
_onset_snap_by_day = Dict(
    obs.n - value(obs.cutoff - b.report_date) => b
        for b in _onset_readings.snaps
)
_onset_cells_by_report = Dict{Int, Vector{Int}}()
for (i, r) in enumerate(obs.onset_curve_history.report_days)
    push!(get!(_onset_cells_by_report, r, Int[]), i)
end
_onset_hazard = fitted_onset_hazard(fit_model("joint"), chn_joint)
_onset_daily_draws = onset_daily_draws(chn_joint)
_onset_replicated = onset_bar_replicator(
    chn_joint, Random.MersenneTwister(20260729)
)

# Each snapshot is nowcast to the delay of the figure each of its onset
# dates was last printed on, so the band and the latest reading are the
# same quantity.
_onset_panels = map(sort(collect(keys(_onset_cells_by_report)))) do R
    snap = _onset_snap_by_day[R]
    us = sort(obs.onset_curve_history.onset_days[_onset_cells_by_report[R]])
    observed = Float64[get(snap.onsets, grid_date(u), 0) for u in us]
    nowcast = onset_nowcast_draws(
        us, observed, [R - u for u in us],
        _onset_daily_draws, _onset_hazard; grid_start = _onset_grid_start,
        target_delays = [_onset_readings.last_report_day[u] - u for u in us]
    )
    (;
        title = string(snap.report_date), dates = grid_date.(us), observed,
        nowcast = [_onset_replicated(d) for d in nowcast],
        latest = [_onset_readings.last_printed[u] for u in us],
    )
end

onset_fit_fig = plot_onset_nowcast_grid(_onset_panels);

Posterior correlations and stream totals

Which headline quantities trade off against each other, and whether the stream totals match the observed ones.

Posterior correlation heatmap

julia

correlation_fig = plot_correlation_heatmap(
    chn_joint,
    [
        :C_T, :R_T, :T, :CFR, :p_drc, :p_uganda, :lambda_bg, :tau_test,
        :expected_reports_T, :expected_deaths_T, :expected_confirmed_T,
    ];
    labels = Dict(
        :C_T => raw"C_T", :R_T => raw"R_T", :T => raw"T",
        :CFR => raw"\mathrm{CFR}", :p_drc => raw"p_\mathrm{drc}",
        :p_uganda => raw"p_\mathrm{ug}", :lambda_bg => raw"\lambda_\mathrm{bg}",
        :tau_test => raw"\tau_\mathrm{test}",
        :expected_reports_T => raw"\mathrm{susp.\ cases}",
        :expected_deaths_T => raw"\mathrm{susp.\ deaths}",
        :expected_confirmed_T => raw"\mathrm{conf.\ cases}"
    )
);

Stream totals against observed

julia

# Per-draw modelled total of each stream, summed over its own reporting
# vintages (the confirmed total adds the unscored first-vintage baseline),
# reusing the posterior-predictive replicates built for the vintage panels.
_stream_total(reps) = [sum(Float64.(collect(r))) for r in vec(reps)]
_conf_baseline = isempty(obs.confirmed_history.counts) ? 0 :
    Int(obs.confirmed_history.counts[1])
stream_totals = (;
    suspected_cases = _stream_total(reported_panel.replicates),
    suspected_deaths = _stream_total(deaths_panel.replicates),
    confirmed_cases = _stream_total(confirmed_panel.replicates) .+ _conf_baseline,
    confirmed_deaths = _stream_total(confirmed_deaths_panel.replicates),
    analysed = _stream_total(tests_analysed_panel.replicates),
);
stream_observed = (;
    suspected_cases = Float64(obs.reported_history.counts[end]),
    suspected_deaths = Float64(obs.deaths_history.counts[end]),
    confirmed_cases = Float64(obs.confirmed_cases),
    confirmed_deaths = Float64(obs.confirmed_deaths_history.counts[end]),
    analysed = Float64(obs.lab_history.counts[end]),
);
stream_pairs_fig = plot_stream_pairs(stream_totals, stream_observed);

Outbreak size estimated by each data stream

The table below puts the posteriors over the infection count side by side, the single-stream fits and the joint, to show what each stream implies alone and what the joint adds.

Per-stream infection-count table

julia

streams_C_table = streams_table(
    "exports" => posterior_C_exports,
    "deaths (DRC)" => posterior_C_deaths,
    "cases (DRC)" => posterior_C_cases,
    "confirmed (DRC)" => posterior_C_confirmed,
    "isolation (DRC)" => posterior_C_treatment,
    "onsets (DRC)" => posterior_C_onsets,
    "joint" => posterior_C_joint
);

Stream	Lower 90%	Lower 60%	Lower 30%	Upper 30%	Upper 60%	Upper 90%
exports	4261	24650	82751	753193	3259686	22461433
deaths (DRC)	32687	65113	98150	206392	325116	926798
cases (DRC)	43456	47549	51692	64628	79283	134807
confirmed (DRC)	58645	69146	78311	105958	130201	228859
isolation (DRC)	12299	14597	16234	20308	24685	35718
onsets (DRC)	15697	23625	31255	50812	65694	115560
joint	12371	14885	16574	19815	21957	26103

The first figure shows each single-stream fit's cumulative-infection trajectory projected to the cut-off, with a dotted rule in each stream's colour marking where its data stops and the ribbon beyond it becomes a forward projection. The count axis is cropped to twice the joint fit's 90% upper bound, as the density figure below is, and a stream whose band runs past the crop is marked with an open triangle where it leaves the axis.

Per-stream projected-trajectory plot

julia

# Per-draw cumulative-infection trajectory carried by each single-stream
# fit out to the cut-off on day `n`, so streams whose data ends earlier are
# still projected to today.
function _cuminf(chn)
    mat = chn[:cumulative_infections]
    return [collect(v) for v in vec(collect(mat))]
end
# Grid day a stream's data last reports, used for the dotted rule. The
# suspected case and death histories freeze at 26 May; exports and confirmed
# run to the cut-off.
_last_day(days) = isempty(days) ? nothing : maximum(days)

stream_traj_fig = plot_stream_trajectories(
    [
        (;
            label = "exports", trajs = _cuminf(chn_exports),
            last_day = _last_day(
                vcat(
                    obs.export_case_days,
                    obs.export_death_days
                )
            ), colour = :seagreen,
        ),
        (;
            label = "deaths (DRC)", trajs = _cuminf(chn_deaths),
            last_day = _last_day(obs.deaths_history.days),
            colour = :firebrick,
        ),
        (;
            label = "cases (DRC)", trajs = _cuminf(chn_cases),
            last_day = _last_day(obs.reported_history.days),
            colour = :steelblue,
        ),
        (;
            label = "confirmed (DRC)", trajs = _cuminf(chn_confirmed),
            last_day = _last_day(obs.confirmed_history.days),
            colour = :goldenrod,
        ),
        (;
            label = "isolation (DRC)", trajs = _cuminf(chn_treatment),
            last_day = _last_day(obs.isolation_history.days),
            colour = :darkorange,
        ),
        (;
            label = "onsets (DRC)", trajs = _cuminf(chn_onsets),
            last_day = _last_day(obs.onset_curve_history.report_days),
            colour = :mediumpurple,
        ),
    ];
    n = obs.n, seeding = obs.seeding,
    # Twice the joint fit's 90% upper, the crop the cut-off density figure
    # below already uses. The exports-only fit bounds the infection count
    # so weakly that its 90% upper reaches the source population, which on
    # a free axis puts every other stream on the baseline.
    ymax = 2.0 * quantile(posterior_C_joint, 0.95)
);

The second figure is the posterior density of each fit's cumulative infection count at the cut-off. The x-axis is scaled to a multiple of the joint-fit 90% upper bound so the bulk of the streams stays visible rather than being flattened by the wide, ill-defined confirmed-only tail.

Cut-off infection-count density plot

julia

# Scale the x-axis to twice the joint-fit 90% upper bound, so the joint and
# the streams that track it read clearly while the confirmed-only tail runs
# off the axis rather than dominating it.
density_xmax = 2.0 * quantile(posterior_C_joint, 0.95)

cumulative_density_fig = plot_cumulative_cases(
    "exports" => posterior_C_exports,
    "deaths (DRC)" => posterior_C_deaths,
    "cases (DRC)" => posterior_C_cases,
    "confirmed (DRC)" => posterior_C_confirmed,
    "isolation (DRC)" => posterior_C_treatment,
    "onsets (DRC)" => posterior_C_onsets,
    "joint" => posterior_C_joint;
    scenarios = [], xmax = density_xmax
);

Reproduction number estimated by each data stream

The reproduction number each stream implies on its own, one panel per stream with the joint fit overlaid in grey as the reference.

Per-stream implied-Rt plot

julia

# The per-stream fits walk Rt from day 1 (the default `rt_start`), while the
# joint walks from `RT_WALK_LEAD` days before the first situation report; the
# shared `display_start` is the joint renewal start so every stream reads over
# the same established window. `ramp` matches the joint Rt figure.
_rt_walk_start_joint = clamp(_BREAKPOINT - RT_WALK_LEAD, _rt_start_plot, obs.n);
stream_rt_fig = plot_rt_streams(
    [
        (;
            label = "exports", chn = chn_exports, rt_start = 1,
            rt_walk_start = 1, colour = :seagreen,
        ),
        (;
            label = "deaths (DRC)", chn = chn_deaths, rt_start = 1,
            rt_walk_start = 1, colour = :firebrick,
        ),
        (;
            label = "cases (DRC)", chn = chn_cases, rt_start = 1,
            rt_walk_start = 1, colour = :steelblue,
        ),
        (;
            label = "confirmed (DRC)", chn = chn_confirmed, rt_start = 1,
            rt_walk_start = 1, colour = :goldenrod,
        ),
        (;
            label = "isolation (DRC)", chn = chn_treatment, rt_start = 1,
            rt_walk_start = 1, colour = :darkorange,
        ),
        (;
            label = "onsets (DRC)", chn = chn_onsets, rt_start = 1,
            rt_walk_start = 1, colour = :mediumpurple,
        ),
    ];
    joint = (;
        label = "joint", chn = chn_joint, rt_start = _rt_start_plot,
        rt_walk_start = _rt_walk_start_joint,
    ),
    n = obs.n, breakpoint = _BREAKPOINT,
    as_of_date = string(obs.cutoff), seeding = obs.seeding,
    display_start = _rt_start_plot, ramp = RT_INTERVENTION_RAMP
);

Sensitivity to assumptions

The joint model re-fitted with one assumption changed at a time.

Delay sensitivity

The death stream dates the outbreak from how far deaths lag symptom onset, so the assumed onset-to-death delay sets the implied infection count. The baseline uses the hospital-pathway delay from the Isiro 2012 line-list reanalysis (onset to admission then admission to death, implied mean about 12 d). We re-fit the joint model under the community-pathway delay from the same reanalysis: the delay for deaths that occur in the community without a recorded admission. This delay is shorter (implied mean about 8 d). Both pathways come from the line list, so this varies the actual delay assumption rather than an arbitrary scenario. The re-fit is the headline model, with its provinces and in-care split, and changes only the delay. It uses the headline's sampler settings.

The infection count to date shifts with the assumed delay, and the table and overlaid densities below show how far.

Re-fit the joint under the community-pathway onset-to-death delay

julia

# The sensitivity re-fits (community-delay variant) are
# defined in the fit registry (`docs/fits/registry.jl`) and loaded through the cache
# (when enabled) in the setup block above.
posterior_C_community_delay = RUN_SENSITIVITY ?
    vec(Array(chn_joint_community_delay[:C_T])) : nothing;

Delay-sensitivity infection-count table

julia

delay_sensitivity_table = RUN_SENSITIVITY ?
    streams_table(
        "baseline (hospital pathway)" => posterior_C_joint,
        "community pathway" => posterior_C_community_delay
    ) :
    Markdown.md"_Delay sensitivity analysis not shown in this build._";

Stream	Lower 90%	Lower 60%	Lower 30%	Upper 30%	Upper 60%	Upper 90%
baseline (hospital pathway)	12371	14885	16574	19815	21957	26103
community pathway	12494	14706	16277	19096	20991	24755

Delay-sensitivity infection-count density plot

julia

delay_sensitivity_fig = RUN_SENSITIVITY ?
    plot_cumulative_cases(
        "baseline (hospital pathway)" => posterior_C_joint,
        "community pathway" => posterior_C_community_delay; scenarios = []
    ) :
    Markdown.md"_Delay sensitivity analysis not shown in this build._";

Tree-prior sensitivity

The outbreak-age estimate depends on the coalescent tree prior assumed in the BEAST X analysis. The baseline uses the more flexible Skygrid non-parametric model, which dates the common ancestor to 15 March 2026 ( HPD 09 Feb – 12 Apr). The report also fits an Exponential growth tree prior, which dates the common ancestor about a week earlier to 08 March 2026 ( HPD 01 Feb – 05 Apr) (Mbala-Kingebeni and others, 2026). Both priors give similar evolutionary rates ( subs/site/year). We re-fit the joint model under the Exponential growth TMRCA and compare the infection count to date and the outbreak age. As for the delay, the re-fit is the headline model with only the common-ancestor date changed.

Re-fit the joint under the Exponential growth tree prior

julia

# The Exponential-growth re-fit (and its `tmrca_days` offset) is defined in the fit
# registry (`docs/fits/registry.jl`) and loaded through the cache (when enabled) in the
# setup block above.
posterior_C_exp_growth = RUN_SENSITIVITY ?
    vec(Array(chn_joint_exp_growth_clock[:C_T])) : nothing
T_skygrid = vec(Array(chn_joint[:T]))
T_exp_growth = RUN_SENSITIVITY ? vec(Array(chn_joint_exp_growth_clock[:T])) : nothing;

The infection count to date under the two tree priors, side by side. A slightly earlier common ancestor (Exponential growth) permits a marginally older outbreak, though the difference is small because the evolutionary rates are nearly identical.

Tree-prior infection-count table

julia

clock_sensitivity_C_table = RUN_SENSITIVITY ?
    streams_table(
        "Skygrid (baseline)" => posterior_C_joint,
        "Exponential growth" => posterior_C_exp_growth
    ) :
    Markdown.md"_Tree-prior sensitivity analysis not shown in this build._";

Stream	Lower 90%	Lower 60%	Lower 30%	Upper 30%	Upper 60%	Upper 90%
Skygrid (baseline)	12371	14885	16574	19815	21957	26103
Exponential growth	12682	15024	16520	20010	22478	26931

Tree-prior infection-count density plot

julia

clock_sensitivity_C_fig = RUN_SENSITIVITY ?
    plot_cumulative_cases(
        "Skygrid (baseline)" => posterior_C_joint,
        "Exponential growth" => posterior_C_exp_growth; scenarios = []
    ) :
    Markdown.md"_Tree-prior sensitivity analysis not shown in this build._";

The outbreak age, the number of days from seeding to the cut-off, under the two tree priors.

Tree-prior outbreak-age table

julia

clock_sensitivity_T_table = RUN_SENSITIVITY ?
    streams_table(
        "Skygrid (baseline)" => T_skygrid,
        "Exponential growth" => T_exp_growth; digits = 0
    ) :
    Markdown.md"_Tree-prior sensitivity analysis not shown in this build._";

Stream	Lower 90%	Lower 60%	Lower 30%	Upper 30%	Upper 60%	Upper 90%
Skygrid (baseline)	184	189	192	198	202	211
Exponential growth	191	195	198	205	210	220

Tree-prior outbreak-age density plot

julia

clock_sensitivity_T_fig = RUN_SENSITIVITY ?
    plot_density_overlay(
        "Skygrid (baseline)" => T_skygrid,
        "Exponential growth" => T_exp_growth;
        xlabel = "Outbreak age (days before cut-off)",
        title = "Posterior outbreak age by tree prior", lower = 0
    ) :
    Markdown.md"_Tree-prior sensitivity analysis not shown in this build._";

Fit diagnostics by parameter

One parameter or the whole model

Per-parameter diagnostics for every fit

julia

# R-hat and both effective sample sizes over several thousand parameters
# are not free to compute, so each fit's per-parameter frame is built once
# here and handed to every table and figure in this section.
diagnostic_fits = [
    "joint" => chn_joint,
    "joint, no patches" => chn_no_patches,
    "exports" => chn_exports,
    "deaths (DRC)" => chn_deaths,
    "cases (DRC)" => chn_cases,
    "confirmed (DRC)" => chn_confirmed,
    "confirmed deaths (DRC)" => chn_confirmed_deaths,
    "isolation (DRC)" => chn_treatment,
    "onsets (DRC)" => chn_onsets,
    "frozen (1wk back)" => frozen_lastweek.chn,
    (
        RUN_SENSITIVITY ?
            [
                "delay sensitivity" => chn_joint_community_delay,
                "clock sensitivity (ExpGrowth)" => chn_joint_exp_growth_clock,
            ] :
            []
    )...,
]
diagnostic_frames = [
    label => parameter_diagnostics(chn)
        for (label, chn) in diagnostic_fits
]
diagnostic_frame = Dict(diagnostic_frames)
joint_diagnostics = diagnostic_frame["joint"]
diagnostic_spread = MarkdownTable(
    diagnostic_spread_table(diagnostic_frames...; labels = display_names)
);

fit	parameters	rhat_above_1.01	rhat_above_1.1	percent_above_1.1	ess_bulk_below_100	lowest_ess_parameter
joint	6703	1232	1	0	140	onset_report_state.σ_h0
joint, no patches	2093	139	0	0	3	onset_report_state.σ_h0
exports	295	3	0	0	0	asc_state.τ_logit
deaths (DRC)	511	2	0	0	0	intervention effect
cases (DRC)	515	0	0	0	0	intervention effect
confirmed (DRC)	434	1	0	0	0	intervention effect
confirmed deaths (DRC)	443	1	0	0	0	asc_state.τ_logit
isolation (DRC)	349	6	0	0	0	treatment_state.adm_delay_state.delay_sd
onsets (DRC)	504	10	0	0	1	Rt step size
frozen (1wk back)	6441	3986	771	12	2118	onset_report_state.σ_h0
delay sensitivity	6709	2529	2	0	36	onset_report_state.σ_h0
clock sensitivity (ExpGrowth)	6801	1685	0	0	146	bartlett_diag

R-hat spread figure

julia

rhat_spread_fig = plot_rhat_spread(
    "joint" => joint_diagnostics,
    "cases (DRC)" => diagnostic_frame["cases (DRC)"],
    "deaths (DRC)" => diagnostic_frame["deaths (DRC)"],
    "confirmed (DRC)" => diagnostic_frame["confirmed (DRC)"],
    "exports" => diagnostic_frame["exports"],
    "frozen (1wk back)" => diagnostic_frame["frozen (1wk back)"]
);

Which parameters mix worst

Worst-mixing parameters of the joint fit

julia

joint_worst_parameters = MarkdownTable(
    worst_parameters_table(
        joint_diagnostics; n = 15,
        labels = display_names
    )
);

parameter	rhat	ess_bulk	ess_tail
onset_report_state.σ_h0	1.103	28	57
onset_report_state.z_h0[6]	1.076	39	124
onset_report_state.z_h0[5]	1.076	40	267
onset_report_state.z_h0[8]	1.07	42	309
onset_report_state.z_h0[7]	1.064	43	200
onset_report_state.z_h0[10]	1.073	45	321
onset_report_state.η0	1.05	46	58
onset_report_state.z_h0[4]	1.061	47	322
onset_report_state.z_h0[2]	1.048	59	279
onset_report_state.z_h0[3]	1.052	62	298
infections_patch[710]	1.032	63	188
infections_patch[714]	1.029	67	181
infections_patch[706]	1.029	69	170
bg_split_state.τ_bg	1.044	69	139
bg_split_state.z_bg[1]	1.037	71	150

The same diagnostics grouped by parameter rather than by element.

Worst-mixing parameters, grouped

julia

joint_worst_groups = MarkdownTable(
    family_diagnostics_table(
        joint_diagnostics; n = 12,
        labels = display_names
    )
);

parameter	elements	max_rhat	min_ess_bulk	above_1.1
onset_report_state.σ_h0	1	1.103	28	1
onset_report_state.z_h0	27	1.076	39	0
onset_report_state.η0	1	1.05	46	0
infections_patch	773	1.032	63	0
bg_split_state.τ_bg	1	1.044	69	0
bg_split_state.z_bg	3	1.04	71	0
cases_state.bg_state.steps	22	1.037	77	0
treatment_state.cap_share_state.τ_cap	1	1.028	81	0
treatment_state.cap_share_state.z_cap	3	1.024	81	0
onset_report_state.σ_a	1	1.012	88	0
delta_patch_start	4	1.04	89	0
delta_knots	96	1.04	89	0

Mixing over time varying paramters

julia

joint_index_fig = plot_parameter_index_diagnostics(
    joint_diagnostics;
    n_groups = 3, labels = display_names
);

Where the divergent transitions sit

Sampler behaviour by chain

julia

joint_chain_table = MarkdownTable(sampler_by_chain_table(chn_joint));

chain	draws	divergences	percent_divergent	step_size	deepest_tree
1	1000	0	0	0.00182	10
2	1000	14	1.4	0.00157	10

Divergence location table

julia

joint_divergence_table = MarkdownTable(
    divergence_location_table(chn_joint; n = 12, labels = display_names)
);

parameter	all_draws	divergent_draws	separation
recovered_state.delay_state.delay_sd	1.56–12.1	1.1–12.8	-1.03
deaths_state.asc_state.p_death	0.786–0.954	0.751–0.926	-0.91
confirmed_deaths_state.expected_confirmed_deaths	3440.0–3900.0	3660.0–3980.0	0.88
deaths_state.asc_state.logit_p_death	1.3–3.03	1.1–2.53	-0.86
lab_composition_state.ρ	0.133–0.254	0.128–0.236	-0.73
onset-to-detection shape	0.789–1.82	0.785–1.49	-0.66
treatment_state.expected_confirmed_incare	240.0–346.0	269.0–348.0	0.65
tau_death	0.449–0.826	0.475–0.963	0.63
treatment_state.incare_confirm_log	-0.33–0.578	-0.481–0.345	-0.62
intervention effect	-0.635–-0.032	-0.559–-0.0219	-0.61
treatment_state.ruleout_los_state.delay_mean	5.84–7.28	5.86–6.83	-0.61
confirmed_deaths_state.spec_state.spec	0.926–0.995	0.944–0.994	-0.6

Divergent draws against the posterior

julia

joint_divergence_fig = plot_divergence_locations(
    chn_joint,
    [:C_T, :R_T, :r, :T, :CFR, :k];
    labels = Dict(
        :C_T => "cumulative infections",
        :R_T => "reproduction number at the cut-off",
        :r => "latest growth rate", :T => "outbreak age",
        :CFR => "case-fatality ratio",
        :k => "surveillance dispersion"
    )
);

The joint fit against the single-stream fits

Joint against single-stream contrast

julia

stream_contrast = diagnostic_contrast(
    "joint" => joint_diagnostics,
    "exports" => diagnostic_frame["exports"],
    "deaths (DRC)" => diagnostic_frame["deaths (DRC)"],
    "cases (DRC)" => diagnostic_frame["cases (DRC)"],
    "confirmed (DRC)" => diagnostic_frame["confirmed (DRC)"],
    "isolation (DRC)" => diagnostic_frame["isolation (DRC)"],
    "onsets (DRC)" => diagnostic_frame["onsets (DRC)"]
)
stream_contrast_table = MarkdownTable(
    diagnostic_contrast_table(
        stream_contrast; n = 15,
        labels = display_names
    )
);

fit	parameter	ess_bulk	ess_bulk_reference	ess_ratio
exports	deaths_state.cfr_state.CFR	2103	118	0.06
exports	rt_state.z[13]	2774	211	0.08
isolation (DRC)	treatment_state.cap_state.z[16]	2022	179	0.09
exports	rt_state.z[7]	2116	189	0.09
exports	rt_state.z[12]	1934	174	0.09
exports	asc_state.p_drc	1305	117	0.09
exports	asc_state.z_drc	2588	243	0.09
exports	rt_state.log_R0	1085	105	0.1
isolation (DRC)	treatment_state.occupancy_step[2]	1620	157	0.1
cases (DRC)	asc_state.p_drc	1205	117	0.1
cases (DRC)	susceptible_fraction[210]	1150	117	0.1
cases (DRC)	susceptible_fraction[211]	1149	117	0.1
cases (DRC)	susceptible_fraction[209]	1150	117	0.1
cases (DRC)	susceptible_fraction[206]	1152	118	0.1
cases (DRC)	susceptible_fraction[205]	1153	118	0.1

Joint against single-stream figure

julia

stream_contrast_fig = plot_diagnostic_contrast(
    stream_contrast;
    xlabel = "Bulk effective sample size, single-stream fit",
    ylabel = "Bulk effective sample size, joint fit",
    title = "Mixing in the joint against each stream fitted alone"
);

The joint fit against the same fit a week earlier

Live against frozen contrast

julia

frozen_contrast = diagnostic_contrast(
    "joint" => joint_diagnostics,
    "one week earlier" => diagnostic_frame["frozen (1wk back)"]
)
frozen_contrast_table = MarkdownTable(
    diagnostic_contrast_table(
        frozen_contrast; n = 15,
        labels = display_names
    )
);

fit	parameter	ess_bulk	ess_bulk_reference	ess_ratio
one week earlier	cases_state.bg_state.steps[10]	210	77	0.37
one week earlier	recovered_state.rec_state.recovery_offset	272	117	0.43
one week earlier	province_capacity_share[3]	1602	739	0.46
one week earlier	province_capacity_share[3]	1570	739	0.47
one week earlier	province_capacity_share[3]	1561	739	0.47
one week earlier	province_capacity_share[3]	1482	739	0.5
one week earlier	province_capacity_share[3]	1464	739	0.5
one week earlier	confirmed_state.δ0	340	176	0.52
one week earlier	province_capacity_share[3]	1302	739	0.57
one week earlier	onset_report_state.z_h0[11]	139	82	0.59
one week earlier	infections_patch[710]	108	63	0.59
one week earlier	province_capacity_share[3]	1200	739	0.62
one week earlier	infections_patch[714]	107	67	0.63
one week earlier	infections_patch[706]	106	69	0.65
one week earlier	onset_report_state.η0	71	46	0.65

Live against frozen figure

julia

frozen_contrast_fig = plot_diagnostic_contrast(
    frozen_contrast;
    xlabel = "Bulk effective sample size, fit a week earlier",
    ylabel = "Bulk effective sample size, live fit",
    title = "Mixing in the live fit against the same fit a week earlier"
);

Parameter recovery

Whether the model recovers known values when fitted to data it simulated itself. Each seed is one prior draw of the model run past the cut-off, kept when its outbreak size is within a factor of five of the one observed, and fitted with the headline joint's sampler settings. The top panel shows each seed's posterior median with its 50% and 90% intervals divided by that seed's true value, so a recovered quantity straddles the line at one. The growth rate is shown as the ratio of daily growth factors, . The intervention effect, a change in , is shown the same way, as the ratio of the multipliers it implies. Below, each quantity is on its own scale: the prior in grey, each seed's posterior in its colour and its true value as a dashed line in the same colour. The prior is the fitted model's, before the factor-of-five selection of the truths. The forecasts are scored against the simulated future and a persistence baseline, where a relative CRPS below one beats the baseline.

Recovery figure

julia

recovery = recovery_results()
recovery_national_quantities = [
    "C_T", "T", "R_T", "r", "CFR", "p_drc", "tau_test", "lambda_bg",
    "growth_state.G", "rt_state.sigma_rw", "onset_report_state.τ",
    "rt_state.intervention_effect",
]
recovery_labels = Dict(
    "growth_state.G" => "G", "rt_state.sigma_rw" => "Rt step size",
    "onset_report_state.τ" => "onset read SD",
    "rt_state.intervention_effect" => "intervention effect",
)
recovery_fig = isempty(recovery.params) ? nothing : plot_recovery(
        recovery.params, recovery.draws, recovery.prior;
        quantities = recovery_national_quantities,
        labels = merge(
            recovery_labels,
            Dict(
                "r" => "r (as exp(r))",
                "rt_state.intervention_effect" => "intervention effect (as exp)",
            )
        ),
        panel_labels = merge(recovery_labels, Dict("r" => "r")),
        log_x = ["C_T", "lambda_bg"],
        difference = ["r", "rt_state.intervention_effect"]
    );

The error of each seed's posterior median relative to the truth, and the z-score of the truth, summarised across seeds.

Summary table across seeds

julia

recovery_summary_national = isempty(recovery.params) ? DataFrame() :
    recovery_summary_table(recovery.params; province = false);

recovery_summary_national_display = isempty(recovery_summary_national) ?
    Markdown.parse("No parameter-recovery run is available for this build.") :
    MarkdownTable(recovery_summary_national);

Quantity	Truth, lowest	Truth, highest	Relative error, median	Relative error, range	z, median	Truth in 90%	Outside 99%
CFR	0.157	0.263	0.07	-0.14 to 0.2	0.41	3/3	0
C_T	3290	12300	0.03	-0.32 to 0.2	0.27	2/3	0
R_T	0.765	1.42	0.02	-0.19 to 0.12	0.29	2/3	0
T	194	228	-0.01	-0.06 to 0.06	-0.25	3/3	0
growth_state.G	14.2	15.5	0	-0.06 to 0.03	0.04	3/3	0
lambda_bg	0.34	8.18	0.16	-0.15 to 0.29	0.45	3/3	0
onset_report_state.τ	0.7	1.04	0.01	-0.04 to 0.02	0.43	3/3	0
p_drc	0.557	0.67	-0.05	-0.06 to 0.45	-0.14	2/3	0
r	-0.018	0.0259	0.18	-0.63 to 0.4	0.4	2/3	0
rt_state.intervention_effect	-0.69	-0.0961	0.42	-2.3 to 0.63	1.35	3/3	0
rt_state.sigma_rw	0.0569	0.105	-0.08	-0.17 to -0.06	-0.12	3/3	0
tau_test	0.536	0.763	0.3	-0.15 to 0.33	1.53	2/3	1

Each seed's fit and recovered values

julia

recovery_seeds = isempty(recovery.params) ? DataFrame() :
    recovery_seed_verdicts(recovery.params);
recovery_national = isempty(recovery.params) ? DataFrame() :
    recovery.params[
        .!occursin.("[", recovery.params.quantity), [
            :seed, :quantity, :truth, :median, :lower_90, :upper_90, :covered_90,
        ],
    ];

recovery_seeds_display = isempty(recovery_seeds) ?
    Markdown.parse("No parameter-recovery run is available for this build.") :
    MarkdownTable(recovery_seeds);
recovery_national_display = isempty(recovery_national) ? Markdown.parse("") :
    MarkdownTable(recovery_national);

seed	status	coverage_90	outside	fit_minutes	max_rhat	min_ess_bulk	divergences
1	unconverged	0.938	tau_test	210.1	1.79	5.323	138
2	pass	0.906		190.3	1.013	424.327	43
3	pass	0.906		142.3	1.009	371.723	4

seed	quantity	truth	median	lower_90	upper_90	covered_90
1	CFR	0.263	0.28	0.212	0.363	true
1	C_T	3285.37	3376.631	2518.957	4690.918	true
1	R_T	0.765	0.855	0.649	1.038	true
1	T	198.243	195.28	183.179	210.755	true
1	growth_state.G	14.205	14.235	12.78	15.725	true
1	lambda_bg	0.34	0.439	0.237	0.991	true
1	onset_report_state.τ	0.7	0.711	0.685	0.748	true
1	p_drc	0.67	0.633	0.448	0.846	true
1	r	-0.018	-0.0108	-0.0285	0.00258	true
1	rt_state.intervention_effect	-0.69	-0.252	-0.762	-0.0894	true
1	rt_state.sigma_rw	0.105	0.0969	0.0543	0.16	true
1	tau_test	0.536	0.696	0.59	0.828	false
2	CFR	0.157	0.188	0.138	0.244	true
2	C_T	7489.771	5060.262	4363.27	6680.522	false
2	R_T	1.254	1.281	1.089	1.551	true
2	T	227.753	214.297	199.442	232.492	true
2	growth_state.G	15.541	14.607	13.033	16.167	true
2	lambda_bg	6.021	6.989	3.22	13.325	true
2	onset_report_state.τ	1.041	1.057	1.001	1.113	true
2	p_drc	0.557	0.809	0.605	0.946	false
2	r	0.015	0.0177	0.00602	0.0323	true
2	rt_state.intervention_effect	-0.68	-0.398	-0.75	-0.0792	true
2	rt_state.sigma_rw	0.0815	0.0678	0.0377	0.121	true
2	tau_test	0.562	0.745	0.54	0.925	true
3	CFR	0.192	0.164	0.102	0.245	true
3	C_T	12312.614	14737.193	9687.327	25539.51	true
3	R_T	1.416	1.147	0.926	1.37	false
3	T	193.694	205.023	190.265	226.487	true
3	growth_state.G	14.154	14.612	13.067	16.158	true
3	lambda_bg	8.179	6.956	3.815	12.074	true
3	onset_report_state.τ	0.743	0.714	0.668	0.761	true
3	p_drc	0.563	0.53	0.306	0.815	true
3	r	0.0259	0.00971	-0.00518	0.0223	false
3	rt_state.intervention_effect	-0.0961	-0.317	-0.672	-0.059	true
3	rt_state.sigma_rw	0.0569	0.0534	0.0116	0.118	true
3	tau_test	0.763	0.647	0.468	0.853	true

Forecast scores from the recovery fits

julia

recovery_forecasts = isempty(recovery.forecasts) ? DataFrame() :
    recovery.forecasts[
        :, [
            :seed, :horizon, :quantity, :truth, :baseline, :crps, :baseline_crps,
            :relative_crps, :covered_90,
        ],
    ];

recovery_forecasts_display = isempty(recovery_forecasts) ?
    Markdown.parse("No recovery forecast is available for this build.") :
    MarkdownTable(recovery_forecasts);

seed	horizon	quantity	truth	baseline	crps	baseline_crps	relative_crps	covered_90
1	7	confirmed_deaths_new	33	9	7.913	24	0.33	true
1	7	confirmed_new	45	2535	4.135	2490	0.00166	true
1	14	confirmed_deaths_new	49	37	5.533	12	0.461	true
1	14	confirmed_new	76	4945	10.909	4869	0.00224	true
2	7	confirmed_deaths_new	16	22	3.561	6	0.593	true
2	7	confirmed_new	98	2235	53.806	2137	0.0252	false
2	14	confirmed_deaths_new	42	41	3.208	1	3.208	true
2	14	confirmed_new	277	4376	58.22	4099	0.0142	true
3	7	confirmed_deaths_new	152	117	10.602	35	0.303	true
3	7	confirmed_new	512	2095	50.642	1583	0.032	true
3	14	confirmed_deaths_new	291	229	11.319	62	0.183	true
3	14	confirmed_new	1069	4060	108.385	2991	0.0362	true

Saving in-sample outputs

The per-stream infection-count table is written to the shared output directory.

Write in-sample outputs

julia

output_dir = get(
    ENV, "BVD_OUTPUT_DIR",
    joinpath(pkgdir(BVDOutbreakSize), "output")
)
mkpath(output_dir)
CSV.write(
    joinpath(output_dir, "cumulative_cases_by_stream.csv"),
    streams_C_table
)

"/home/runner/work/BVDOutbreakSize/BVDOutbreakSize/output/cumulative_cases_by_stream.csv"

Write the summary bullets

julia

dashboard_dir = joinpath(
    pkgdir(BVDOutbreakSize), "docs", "src", "summary_assets"
)
mkpath(dashboard_dir)

# The bullets under the summary heading at the top of the page. They read
# tables built further down, so they are written here and read back when
# the site is assembled.
evaluation_insample_national_summary = let
    fmt(x) = ismissing(x) || !isfinite(x) ? "n/a" :
        string(round(x; digits = 2))
    cal = filter(r -> isfinite(r["90% coverage"]), stream_calibration_table)
    n_cov = count(>=(0.8), cal[!, "90% coverage"])
    calibrated(r) = string(
        r["Stream"], " (bias ", fmt(r["Bias"]), ", 90% coverage ",
        fmt(r["90% coverage"]), ")"
    )
    worst = first(
        sort(cal, "Bias"; by = abs, rev = true), min(3, size(cal, 1))
    )
    pred(x, observed) = string(
        "observed ", observed, " against a predictive median of ",
        round(Int, quantile(x, 0.5)), " (90% interval ",
        round(Int, quantile(x, 0.05)), "–",
        round(Int, quantile(x, 0.95)), ")"
    )
    overall = [
        string(
            "- **Streams:** ", n_cov, " of ", size(cal, 1),
            " fitted streams have 90% coverage of at least 0.8."
        ),
        string(
            "- **Least well reproduced:** ",
            join(calibrated.(eachrow(worst)), "; "), "."
        ),
        string(
            "- **Exports:** Uganda exports ",
            pred(pp_exports, obs.exported_cases), ", and export deaths ",
            pred(pp_exports_deaths, obs.exports_deaths), "."
        ),
    ]
    # Per stream, split the same way as the posterior predictive checks.
    reporting = Set(p.title for p in reporting_panels)
    function block(lead, keep)
        rows = filter(r -> keep(r["Stream"] in reporting), cal)
        size(rows, 1) == 0 && return nothing
        return join(
            vcat(
                [string("**", lead, "**"), ""],
                [
                    string(
                        "- ", r["Stream"], ": bias ", fmt(r["Bias"]),
                        ", 90% coverage ", fmt(r["90% coverage"]), " over ",
                        r["Vintages"], " vintages."
                    )
                        for r in eachrow(rows)
                ]
            ), "\n"
        )
    end
    blocks = filter(
        !isnothing,
        [
            block("Streams still reporting", identity),
            block("Streams no longer reporting", !),
        ]
    )
    join(vcat([join(overall, "\n")], blocks), "\n\n")
end
write(
    joinpath(dashboard_dir, "evaluation_insample_national.md"),
    evaluation_insample_national_summary
);

In-sample checks ​

Summary ​

Prior predictive check ​

Posterior predictive checks ​

Streams still reporting ​

Cumulative ​

Per-vintage incidence ​

Streams no longer reporting ​

Cumulative ​

Per-vintage incidence ​

Stream calibration ​

Exports ​

Onset snapshot nowcasts ​

Posterior correlations and stream totals ​

Outbreak size estimated by each data stream ​

Reproduction number estimated by each data stream ​

Sensitivity to assumptions ​

Delay sensitivity ​

Tree-prior sensitivity ​

Fit diagnostics by parameter ​

One parameter or the whole model ​

Which parameters mix worst ​

Where the divergent transitions sit ​

The joint fit against the single-stream fits ​

The joint fit against the same fit a week earlier ​

Parameter recovery ​

Saving in-sample outputs ​

In-sample checks

Summary

Prior predictive check

Posterior predictive checks

Streams still reporting

Cumulative

Per-vintage incidence

Streams no longer reporting

Cumulative

Per-vintage incidence

Stream calibration

Exports

Onset snapshot nowcasts

Posterior correlations and stream totals

Outbreak size estimated by each data stream

Reproduction number estimated by each data stream

Sensitivity to assumptions

Delay sensitivity

Tree-prior sensitivity

Fit diagnostics by parameter

One parameter or the whole model

Which parameters mix worst

Where the divergent transitions sit

The joint fit against the single-stream fits

The joint fit against the same fit a week earlier

Parameter recovery

Saving in-sample outputs