BVDOutbreakSize

Appearance

Sensitivity and comparison analyses ​

This page continues the main analysis from the one-week-ahead forecast onward: forecast validation, forecast scoring across releases, the reproduction number by release, the outbreak size each data stream implies alone, how the estimate has evolved across releases, comparisons with McCabe et al. and Chamla et al., and the delay and tree-prior sensitivity re-fits. It renders from the same fitted chains as the main analysis, loaded through the shared setup, so no model is re-fit here beyond the frozen re-fits and the optional sensitivity re-fits below.

Load packages, data and fitted chains
julia
# Shared setup: packages, observations, the fit registry and every model fit
# (loaded from the content-addressed cache). See docs/examples/_setup.jl.
using BVDOutbreakSize
include(joinpath(pkgdir(BVDOutbreakSize), "docs", "examples", "_setup.jl"))
45

Forecast validation ​

How last week's forecast held up against the data since observed, using the frozen re-fit and one-week projection defined in forecast-versus-frozen evaluation. Only the streams the situation reports are still updating are validated here. A stream that has stopped being reported carries a cumulative total that repeats its last reported value, so there is no observation for the past week to score against. The scoring tables further down withhold such a window for the same reason. The projections for those streams are shown separately below. The frozen fit also conditions on the isolation beds, so the projected bed occupancy is scored against the beds held a week later. The bed validation is weak at a one-week-back freeze. The reported occupancy rate starts only on 9 June, so the capacity has no implied-capacity anchor and rides its random walk back to the freeze date. This widens the projected bed interval. Like the scores further down, the confirmed new-count rows here take out any retrospective harmonisation step the week contained. Such a step reattaches records notified earlier, so it is not something the forecast was predicting. A validation week containing no harmonisation day carries no correction, and the rows then read as the raw window counts. The cumulative rows are scored against the published total, harmonisation included.

Fit one week back and validate the one-week-ahead forecast
julia
# frozen_lastweek and frozen_lastweek_streams are computed in the setup
# block above.
# `obs_recovered` is passed so the frozen fit's forecast carries a
# `recovered_new` column (materialised only when the recovered origin is
# given), letting the recovered stream be scored against the observed count
# below like the other streams.
# The onset grid is the one the FROZEN fit saw, not the live one, so the
# validation forecast carries an `onset reports` row scored on the triangle
# the frozen fit was actually fitted to.
_val_onset_days = frozen_lastweek.o.onset_curve_history.onset_days
_val_grid_start = isempty(_val_onset_days) ? nothing :
                  minimum(_val_onset_days)
_val_grid_end = isnothing(_val_grid_start) ? nothing :
                max(maximum(frozen_lastweek.o.onset_curve_history.report_days),
    _val_grid_start)
validation_forecast = forecast_reported(frozen_lastweek.chn;
    horizon = 7,
    obs_cases = frozen_lastweek.o.reported_cases,
    obs_deaths = frozen_lastweek.o.total_deaths,
    obs_confirmed = frozen_lastweek.o.confirmed_cases,
    obs_confirmed_deaths = frozen_lastweek.o.confirmed_deaths,
    obs_recovered = frozen_lastweek.o.recovered_cases,
    grid_n = frozen_lastweek.o.n,
    onset_grid_start = _val_grid_start, onset_grid_end = _val_grid_end);

# Each frozen individual (single-stream) fit's own one-week-ahead new-count
# forecast at the same cut-off as `frozen_lastweek`, from
# [`forecast_stream`](@ref) (the same per-stream forecaster
# `stream_forecasts.csv` uses), so the validation plots below can show the
# individual fit alongside the joint rather than the joint alone. Recovered
# has no individual fit and is absent here, as it is throughout this report.
# Only the still-reported streams are fitted at the validation cut-off, so
# a stream the situation reports have stopped updating is absent from
# `frozen_lastweek_streams` and carries no individual series here.
function _validation_individual_new(sid, stream::Symbol, obs_field)
    haskey(frozen_lastweek_streams, sid) || return nothing
    f = frozen_lastweek_streams[sid]
    bp = f.o.n - f.o.who_first_sitrep_days
    return Float64.(forecast_stream(f.chn, stream; horizon = 7,
        obs_value = getproperty(f.o, obs_field), n = f.o.n, breakpoint = bp,
        rt_start = 1, rt_walk_start = 1))
end
validation_individual = NamedTuple(
    k => v
for (k, v) in pairs((;
        cases_new = _validation_individual_new(
            "cases", :reported_cases, :reported_cases),
        deaths_new = _validation_individual_new(
            "deaths", :suspected_deaths, :total_deaths),
        confirmed_new = _validation_individual_new(
            "confirmed", :confirmed_cases, :confirmed_cases),
        confirmed_deaths_new = _validation_individual_new(
            "confirmed_deaths", :confirmed_deaths, :confirmed_deaths)))
if !isnothing(v))
# The frozen individual (treatment-only) fit's own bed-occupancy forecast,
# anchored on the beds occupied at ITS OWN cut-off (the frozen fit's own
# `o`, not the current `obs`), matching how the joint frozen forecast is
# itself anchored.
# `nothing` when the beds have stopped being reported, so the treatment fit
# is absent; the bed panel then draws the joint alone.
# A `let` block, not a bare `if`: a top-level `if` shares the script's
# global scope, so its working names would leak into the rest of the page.
validation_individual_isolation = let
    if haskey(frozen_lastweek_streams, "treatment")
        tf = frozen_lastweek_streams["treatment"]
        beds = isempty(tf.o.isolation_history.counts) ? 0.0 :
               Float64(tf.o.isolation_history.counts[end])
        Float64.(forecast_stream(tf.chn, :isolation_beds; horizon = 7,
            obs_value = beds, n = tf.o.n,
            breakpoint = tf.o.n - tf.o.who_first_sitrep_days,
            rt_start = 1, rt_walk_start = 1))
    else
        nothing
    end
end

# The observed beds at the current cut-off (the forecast target), so the
# frozen-fit bed forecast is scored against what the beds actually held.
# Held back once the beds stop being reported, since the last count would
# then be carried forward rather than observed at the target date.
_obs_beds = stream_reporting(obs, :isolation_beds) ?
            obs.isolation_history.counts[end] : missing
# Same observed/baseline keying as the plot below, so the table covers every
# fitted count stream (cumulative and new-count rows) plus the bed level.
# A harmonisation-break day between the frozen cut-off and the current one
# puts records into the confirmed cumulative that were never notified in that
# week, so the new-count truth carries a step the forecast was never
# predicting. Take it out, the same correction `score_releases.jl` applies.
# Grid days are relative to a seeding date fixed by the genetic tmrca, so the
# frozen fit's own `n` and the current `obs.n` index the same grid.
validation_breaks = (
    confirmed_cum = confirmed_break_correction(
        obs, frozen_lastweek.o.n, obs.n),
    confirmed_deaths_cum = confirmed_break_correction(
        obs, frozen_lastweek.o.n, obs.n; deaths = true))

# Observed cumulative at the target date per stream, keyed by the forecast's
# cumulative column; `baseline` is each stream's origin cumulative (the
# frozen cut-off), so the new count is scored against observed minus origin,
# less any harmonisation the window carries (see `validation_breaks`). Both
# the table and the plot below take the still-reported streams
# (`reporting_cum_cols`, from the setup block): a stream the situation
# reports have stopped updating has an origin and a target reading the same
# repeated total, so its cumulative truth is stale and its new-count truth is
# a guaranteed zero.
validation_observed = (cases_cum = obs.reported_cases,
    deaths_cum = obs.total_deaths,
    confirmed_cum = obs.confirmed_cases,
    confirmed_deaths_cum = obs.confirmed_deaths,
    recovered_cum = obs.recovered_cases)
validation_baseline = (cases_cum = frozen_lastweek.o.reported_cases,
    deaths_cum = frozen_lastweek.o.total_deaths,
    confirmed_cum = frozen_lastweek.o.confirmed_cases,
    confirmed_deaths_cum = frozen_lastweek.o.confirmed_deaths,
    recovered_cum = frozen_lastweek.o.recovered_cases)

validation_table = forecast_vs_truth(validation_forecast;
    observed = keep_streams(validation_observed, reporting_cum_cols),
    baseline = keep_streams(validation_baseline, reporting_cum_cols),
    breaks = validation_breaks,
    isolation = _obs_beds);
Forecast-versus-observed validation table
julia
# `MarkdownTable` rather than a bare table expression: a DataFrame is
# `text/html`-showable, Literate prefers that mime, and the `@raw html`
# block it writes crosses Documenter's raw-block regex limit once the
# table grows. `MarkdownTable` is markdown-showable and not
# html-showable, so the table goes out as an ordinary markdown table
# rather than a fixed-width block of printed output. See its docstring
# for the mechanism. The same treatment is applied to every DataFrame
# display in this file and in `analysis.jl`.
StreamQuantityObservedLower 90%Lower 60%Lower 30%Upper 30%Upper 60%Upper 90%Within 90% PI
DRC confirmed casescumulative by T+76686656066486712682769017050yes
DRC confirmed casesnew this week586460548612727801950yes
DRC confirmed deathscumulative by T+73226317632123238328233083367yes
DRC confirmed deathsnew this week276226262288332358417yes
DRC recoveredcumulative by T+71563155916971761187119502113yes
DRC recoverednew this week180176314378488567730yes
DRC isolation bedsoccupancy at T+78196187197728579231038yes

The observation panels histogram the one-week-ahead forecast made from the frozen fit: a cumulative and a new-count panel for each still-reported count stream the forecast carries. The 90% predictive interval is shaded, and the count observed by the current cut-off is a dashed black rule. Each stream draws only when the forecast carries its column and the observation covers the target date, so a fit observing fewer streams shows fewer panels. Where a stream has its own individual (single-stream) fit, that fit's forecast from the same frozen cut-off is overlaid as a dotted density alongside the joint's histogram. Recovered has no individual fit and draws the joint alone.

Forecast-versus-observed plot
julia
validation_fig = plot_forecast_vs_truth(validation_forecast;
    observed = keep_streams(validation_observed, reporting_cum_cols),
    baseline = keep_streams(validation_baseline, reporting_cum_cols),
    breaks = validation_breaks,
    individual = keep_streams(validation_individual, reporting_cum_cols));

The bed panel scores last week's projected occupancy against the beds occupied now (the dashed rule), with the individual (treatment-only) fit's own projection overlaid as a dotted density alongside the joint.

Bed forecast-versus-observed plot
julia
validation_beds_fig = plot_forecast_beds_vs_truth(validation_forecast;
    isolation = _obs_beds, individual = validation_individual_isolation);

The latent quantities are not observed, so they are scored distribution against distribution: what the frozen fit forecast for the past week's new infections, onsets and deaths against what the current fit now estimates for the same window.

Forecast-versus-now latent plot
julia
# Current fit's draws of the new latent counts over the past week, the last
# seven days of each cumulative-trajectory deterministic.
function _now_new(chn, key)
    mat = chn[key]
    trajs = [collect(v) for v in vec(collect(mat))]
    return Float64[t[end] - t[max(1, length(t) - 7)] for t in trajs]
end
now_latent = (;
    infections_new = _now_new(chn_joint, :cumulative_infections),
    onsets_new = _now_new(chn_joint, :cumulative_onsets),
    deaths_latent_new = _now_new(chn_joint, :cumulative_expected_deaths))

validation_latent_fig = plot_forecast_vs_truth_latent(
    validation_forecast; now = now_latent);

Streams no longer reported ​

The situation reports have stopped updating some of the streams the model fits, listed with the date each was last reported below. The panels show what the frozen fit projected for those streams over the same week, without an observed rule, since the count they would be scored against has not moved since the stream stopped. These are the model's projections rather than a validation of them.

Forecast for the streams no longer reported
julia
# The last-reported date per stopped stream, and the frozen fit's own
# projection for them. `plot_forecast` draws a panel per new-count column
# the frame carries, so passing the stopped streams' columns alone gives the
# projection without the fabricated truth rule the validation figure would
# otherwise draw against a repeated total.
validation_stopped_streams = let s = stream_report_status(obs),
    ids = [stream_id(c) for c in stopped_cum_cols]

    keep = [r.stream in ids for r in eachrow(s)]
    DataFrame("Stream" => s[keep, :label],
        "Last reported" => s[keep, :last_date])
end
_stopped_new_cols = [c
                     for c in new_cols(stopped_cum_cols)
                     if c in propertynames(validation_forecast)]
validation_stopped_fig = plot_forecast(
    validation_forecast[!, _stopped_new_cols]);
julia
# See the comment above `validation_table`'s display for why this wraps
# the table in `MarkdownTable` instead of showing it directly.

Forecast scoring across releases ​

Every release's saved one- to four-week-ahead forecast is scored against the data observed since, against a persistence baseline and, where one exists, the stream's own individual fit as well as the joint. The tables in this section are the joint model's, one row per stream. Each stream's individual fit is scored the same way and tabulated in Individual fits against the baseline below, so a fit appears in one table rather than two. See forecast scoring against a persistence baseline for how the scores, the relative skill and the baseline are built. Recovered has no individual fit of its own, so its comparison is the baseline against the joint only. Reported cases and suspected deaths stopped being updated by the situation reports partway through the outbreak, and exports' confirmed-detection series is anchored to an earlier cut-off. Exports therefore contributes no scored forecast, and reported cases and suspected deaths each rest on exactly one matched forecast, a single window rather than a settled sample.

Only a minority of the daily releases examined contribute a row to the table below, each a reconstruction of an earlier model version rather than the current fit. The table below is therefore not a verdict on the current fit. One reconstruction is dropped from scoring entirely: its chain forecasts a near-zero median at every horizon and stream, with the upper predictive tail occasionally reaching five- and six-digit values. This is the signature of a chain that failed to sample properly rather than a genuine forecast, so the scoring script flags and excludes it. Only the newest few releases carry the current model's own individual-stream forecasts, and the backfilled reconstructions carry none at all. The comparison against each stream's individual fit therefore rests on those releases alone, filling in one horizon at a time as their targets resolve. Every row also rests on one to a handful of matched forecasts, shown as its own count rather than rounded away. A ratio here should therefore be read as an early signal rather than a settled result.

The symptom-onset stream is scored on the new reported count each vintage adds rather than on its level, because every vintage rereads the whole figure. Its printed total therefore moves with the scan error as well as with late reporting. It appears only from the release that first carried it. Its intervals are dominated by that scan error rather than by epidemic uncertainty, so read its skill against the baseline rather than its coverage.

Load and summarise the cross-release forecast scores
julia
# scripts/score_releases.jl writes these after the fits. The committed files
# are header-only until a release carries the asset, so the common path
# reads a real file to a zero-row frame; the typed `schema` is the fallback
# for a file that is absent entirely, since CSV.read throws on a missing
# path and would take the whole docs build with it.
function _release_data(name, schema::NamedTuple)
    path = joinpath(pkgdir(BVDOutbreakSize), "data", name)
    isfile(path) && return CSV.read(path, DataFrame)
    return DataFrame([k => T[] for (k, T) in pairs(schema)])
end

forecast_scores_df = _release_data("forecast_scores.csv",
    (; release = String, made_date = Date, stream = String, horizon = Int,
        target_date = Date, fit = String, crps = Float64,
        log_crps = Float64, dispersion = Float64, overprediction = Float64,
        underprediction = Float64, coverage_50 = Float64,
        coverage_90 = Float64,
        bias = Float64, n_samples = Int,
        log_rel_to_baseline = Float64))
forecast_overlay_df = _release_data("forecast_overlay.csv",
    (; release = String, made_date = Date, stream = String, horizon = Int,
        target_date = Date, fit = String, observed = Float64,
        median = Float64, lo30 = Float64, hi30 = Float64, lo60 = Float64,
        hi60 = Float64, lo90 = Float64, hi90 = Float64))
# One row per (stream, fit) pooled over every horizon and release. The
# by-horizon and by-release detail tables carry the same columns at a finer
# grain (see src/scoring.jl). Every fit is kept here, since the
# relative-skill figure below compares the roles against each other. The
# tables rendered in this section select the joint role, and the individual
# fits are tabulated in their own section.
forecast_score_overview_table = forecast_score_overview(forecast_scores_df)
forecast_score_by_horizon_table = forecast_score_by_horizon(forecast_scores_df)
forecast_score_by_release_table = forecast_score_by_release(forecast_scores_df)

joint_score_overview_table = select_fit_role(
    forecast_score_overview_table, "joint")
joint_score_by_horizon_table = select_fit_role(
    forecast_score_by_horizon_table, "joint")
# The trailing `;` on this last assignment matters: without it, this whole
# setup chunk's last statement (the DataFrame it assigns) is Literate's
# implicitly displayed "result" for the chunk, on top of the deliberate
# display further down -- and a bare DataFrame is html-showable, so it
# goes out as a second, undisplayed-in-source `@raw html` block that (for
# a table this size) can itself hit the PCRE limit described above.
joint_score_by_release_table = select_fit_role(
    forecast_score_by_release_table, "joint");

The headline pools every horizon and release into one row per stream for the joint model: the mean CRPS and its decomposition, coverage, bias, and the relative skill against the persistence baseline, on both the natural and the log scale. Each row also carries relative skill against the stream's own individual fit where one exists. Column definitions are in forecast scoring against a persistence baseline.

streamfitncrpsrel_to_baselinelog_crpslog_rel_to_baselinerel_to_individuallog_rel_to_individualdispersionoverpredictionunderpredictioncoverage_50coverage_90bias
confirmed casesjoint701159.77.030.2331.790.80.751097.3661.760.580.8410.21
confirmed deathsjoint681226.2610.820.5813.20.921.661099.465.62121.180.380.99-0.46
isolation bedsjoint6573.631.360.0961.30.880.6944.5821.967.10.60.940.19
onset reportsjoint25166.50.890.4670.680.810.8977.4664.1624.890.40.960.19
recoveredjoint11139.41.680.3761.09missingmissing120.026.3113.0711-0.18

The same relative skill against the baseline, by horizon: one panel per stream, one series per fit role, on a log-scaled skill axis with the reference line at one. This is the one place the two roles are drawn against each other, so it carries each stream's individual fit alongside the joint. A fit that beats the baseline on average but not at every cut-off is visible as a series that crosses the line rather than sitting under it throughout.

julia
forecast_relative_skill_fig = plot_forecast_relative_skill(
    forecast_score_by_horizon_table);

The same columns as a table for the joint model, broken out by horizon, and again broken out by release and averaged across horizons, are behind the two dropdowns below.

Scores by horizon
streamhorizonfitncrpsrel_to_baselinelog_crpslog_rel_to_baselinerel_to_individuallog_rel_to_individualdispersionoverpredictionunderpredictioncoverage_50coverage_90bias
confirmed cases7joint22105.012.850.2353.070.810.7282.0722.070.870.7710.24
confirmed cases14joint19234.993.170.2372.970.80.74189.3444.950.70.8910.18
confirmed cases21joint16471.592.180.2191.330.770.75386.8584.340.40.8810.18
confirmed cases28joint135142.9511.380.2380.960.820.815017.12125.690.140.8510.22
confirmed deaths7joint2283.753.570.5945.751.151.3233.340.6549.760.450.95-0.43
confirmed deaths14joint18203.512.990.5693.861.251.9588.041.93113.540.331-0.49
confirmed deaths21joint15393.792.630.5882.71.011.95211.547.22175.040.21-0.5
confirmed deaths28joint135536.3919.350.5661.770.71.645328.6417.29190.460.541-0.43
isolation beds7joint2064.531.550.0831.440.680.5231.2726.257.010.450.950.29
isolation beds14joint1774.391.530.0931.390.930.6944.8724.045.470.590.940.22
isolation beds21joint1578.481.20.1021.151.010.8151.7918.138.570.730.930.12
isolation beds28joint1381.051.20.1141.251.030.9156.3417.047.660.690.920.09
onset reports7joint11121.830.610.7860.640.880.9655.648.7117.520.2710.18
onset reports14joint8212.151.20.2660.910.70.7278.185.748.350.250.880.26
onset reports21joint6187.541.040.1490.730.630.63116.6763.757.120.8310.12
recovered7joint326.962.030.3191.67missingmissing20.931.284.7611-0.18
recovered14joint389.71.950.3931.45missingmissing68.015.8115.8711-0.23
recovered21joint3297.163.240.4221.15missingmissing269.916.0511.211-0.09
recovered28joint2145.990.640.3670.56missingmissing121.86024.1311-0.26
Scores by release
made_datestreamfitncrpsrel_to_baselinelog_crpslog_rel_to_baselinerel_to_individuallog_rel_to_individualdispersionoverpredictionunderpredictioncoverage_50coverage_90bias
2026-06-07confirmed casesjoint2363.9114.180.91311.98missingmissing341.0622.860110.17
2026-06-10confirmed casesjoint3534.876.690.7914.97missingmissing435.6799.20110.27
2026-07-01confirmed casesjoint415862.7874.740.5832.84missingmissing15377.31485.470010.66
2026-07-06confirmed casesjoint4421.281.310.2690.95missingmissing311.85109.4300.510.49
2026-07-08confirmed casesjoint4390.981.230.2530.97missingmissing239.5151.4700.2510.52
2026-07-23confirmed casesjoint4335.570.860.1690.661.020.82299.2735.680.62110.13
2026-07-25confirmed casesjoint4222.50.980.1340.990.640.63218.044.460110.07
2026-07-26confirmed casesjoint4235.691.070.1411.060.790.67229.655.680.36110.06
2026-07-27confirmed casesjoint4247.971.220.1421.070.930.83240.087.880110.11
2026-07-31confirmed casesjoint4259.641.850.1491.690.640.68235.7323.760.15110.13
2026-08-01confirmed casesjoint4231.261.970.1391.990.690.66226.844.310.11110.05
2026-08-02confirmed casesjoint4278.452.120.1451.750.940.85263.7714.680110.14
2026-08-03confirmed casesjoint4247.581.290.140.960.850.8237.889.70110.08
2026-08-04confirmed casesjoint4221.562.430.1332.390.630.69220.70.860110.03
2026-08-07confirmed casesjoint4213.542.390.142.010.690.73209.7603.7811-0.09
2026-08-11confirmed casesjoint3144.143.410.133.340.610.71137.2806.8611-0.15
2026-08-15confirmed casesjoint3277.435.280.1833.471.280.9200.1477.290110.38
2026-08-17confirmed casesjoint2202.192.310.1871.841.40.89125.1277.070110.4
2026-08-22confirmed casesjoint2157.122.860.1612.630.980.84110.9346.20110.44
2026-08-24confirmed casesjoint177.263.260.1272.920.870.8233.9543.310010.63
2026-08-25confirmed casesjoint174.232.90.1212.380.820.7733.6240.610010.64
2026-08-30confirmed casesjoint136.490.840.0580.720.410.4232.573.920110.2
2026-06-07confirmed deathsjoint141.491.870.9032.55missingmissing36.9404.5511-0.2
2026-06-10confirmed deathsjoint264.932.270.7353.05missingmissing62.70.991.2411-0.03
2026-07-01confirmed deathsjoint41719680.450.5491.21missingmissing17103.1192.890110.28
2026-07-06confirmed deathsjoint4197.20.890.8712.29missingmissing110.35086.850.251-0.49
2026-07-08confirmed deathsjoint4177.590.780.7131.84missingmissing84.29093.30.251-0.61
2026-07-23confirmed deathsjoint4271.491.670.4842.450.391.08172.61098.8811-0.45
2026-07-25confirmed deathsjoint4241.831.980.5043.411.022.15106.460135.3701-0.57
2026-07-26confirmed deathsjoint4268.462.310.5544.161.052.17101.590166.8701-0.64
2026-07-27confirmed deathsjoint4249.472.290.5343.990.82.12109.30140.1701-0.59
2026-07-31confirmed deathsjoint4274.552.310.563.960.881.94108.790165.760.251-0.58
2026-08-01confirmed deathsjoint4313.482.690.6834.661.122.794.480219.0101-0.68
2026-08-02confirmed deathsjoint4258.172.850.5074.381.011.83128.190129.980.251-0.53
2026-08-03confirmed deathsjoint4268.891.960.582.911.22.37123.560145.330.251-0.55
2026-08-04confirmed deathsjoint4298.063.910.6567.331.22.3293.540204.5201-0.67
2026-08-07confirmed deathsjoint4366.294.840.888.051.582.7176.060290.2300.75-0.79
2026-08-11confirmed deathsjoint3256.4612.090.82121.351.542.569.10187.3601-0.75
2026-08-15confirmed deathsjoint3154.313.710.3214.390.910.97128.64025.6711-0.28
2026-08-17confirmed deathsjoint2105.062.110.3072.720.971.1383.42021.6411-0.3
2026-08-22confirmed deathsjoint297.181.330.3032.110.740.6374.43022.7511-0.3
2026-08-24confirmed deathsjoint117.411.570.0551.490.180.1152.410110.22
2026-08-25confirmed deathsjoint114.890.920.0490.820.140.0613.831.060110.16
2026-08-30confirmed deathsjoint119.281.930.0651.760.180.0714.165.120110.29
2026-07-01isolation bedsjoint490.191.480.1381.52missingmissing13.48076.7100-0.99
2026-07-06isolation bedsjoint442.870.620.0660.64missingmissing21.07021.801-0.64
2026-07-08isolation bedsjoint454.751.540.0731.54missingmissing11.6743.080010.8
2026-07-23isolation bedsjoint480.062.720.1012.631.441.1429.3950.670010.58
2026-07-25isolation bedsjoint480.842.190.1052.161.561.2832.6548.1900.510.53
2026-07-26isolation bedsjoint473.831.970.0971.91.20.9636.9636.8700.510.46
2026-07-27isolation bedsjoint471.582.470.0932.351.20.9539.5232.0600.510.44
2026-07-31isolation bedsjoint471.051.680.0881.630.860.6751.6919.3600.7510.34
2026-08-01isolation bedsjoint4721.120.0941.090.810.6164.517.480.02110.18
2026-08-02isolation bedsjoint483.691.630.1021.470.990.7252.6731.0200.7510.39
2026-08-03isolation bedsjoint478.831.850.0991.720.990.7460.4218.410110.3
2026-08-04isolation bedsjoint4100.741.430.1371.361.160.9863.4737.2800.7510.38
2026-08-07isolation bedsjoint483.811.440.1221.610.80.7169.03014.7811-0.24
2026-08-11isolation bedsjoint364.560.390.0910.380.620.5261.8802.6811-0.09
2026-08-15isolation bedsjoint373.180.960.0820.850.590.4264.938.250110.17
2026-08-17isolation bedsjoint273.471.560.0831.330.750.5257.8315.640110.32
2026-08-22isolation bedsjoint263.11.560.0671.380.440.352.8110.290110.24
2026-08-24isolation bedsjoint158.732.430.0652.180.440.336.0322.70110.44
2026-08-25isolation bedsjoint139.50.670.0450.610.290.2239.500110
2026-08-30isolation bedsjoint164.572.690.0732.450.550.4234.1230.450010.52
2026-08-02onset reportsjoint383.60.440.1050.331.080.7877.912.752.9511-0.01
2026-08-03onset reportsjoint3160.111.780.3581.241.251.1976.4983.6200.3310.6
2026-08-04onset reportsjoint3124.60.710.3910.771.21.1374.1248.452.030.6710.38
2026-08-07onset reportsjoint3134.381.150.2271.150.60.4969065.380.331-0.51
2026-08-11onset reportsjoint3204.630.660.3080.340.860.783.450121.180.331-0.54
2026-08-15onset reportsjoint3277.32.290.261.220.570.67107.05170.2400.3310.55
2026-08-17onset reportsjoint2211.092.150.2911.20.560.6678.83132.2600.50.50.54
2026-08-22onset reportsjoint2183.240.610.4610.670.670.7967.44115.80010.72
2026-08-24onset reportsjoint1127.780.550.720.720.890.9457.570.290010.6
2026-08-25onset reportsjoint1183.540.43.9640.750.941.1661.19122.350010.88
2026-08-30onset reportsjoint1108.690.580.5360.350.970.7961.1047.5901-0.55
2026-07-01recoveredjoint3312.696.520.4461.86missingmissing289.5523.140110.27
2026-07-06recoveredjoint485.760.790.4070.85missingmissing68.39017.3711-0.33
2026-07-08recoveredjoint463.080.750.2921missingmissing44.51018.5711-0.38

Forecasts made at each release against the value observed since, one panel per stream and horizon, the observed value in black. The median and 90% interval are coloured by fit role: the persistence baseline, the stream's individual fit and the joint. The x-axis is the date each forecast was made, so an incident stream's observed window pairs unambiguously with the forecast that made it. Each panel's axis is cropped to a small multiple of what that stream actually reached, so one very wide interval cannot squash every other series flat. An interval or median too wide for the panel is clamped at the top and marked with an open triangle rather than silently cut off.

Forecasts-versus-now overlay
julia
forecast_overlay_fig = plot_forecast_overlay(forecast_overlay_df);

Frozen-fit forecast evaluation ​

The current model, frozen at earlier data cut-offs (see Forecast-versus-frozen evaluation), is scored the same way as the cross-release forecasts above, against the same persistence baseline. Only the joint model is scored here, so no individual single-stream fit appears in the tables and figures below. The May cut-offs predate the first reported bed occupancy and the first reported recoveries, so those windows are left unscored rather than scored against a series that had not started. The baseline carries a weaker data-vintage guarantee than the cross-release one, since its snapshot was taken weeks after the frozen cut-off and can hold later revisions to earlier days (see forecast scoring against a persistence baseline).

Load and summarise the frozen-fit forecast scores
julia
frozen_scores_df = _release_data("forecast_scores_frozen.csv",
    (; release = String, made_date = Date, stream = String, horizon = Int,
        target_date = Date, fit = String, crps = Float64,
        log_crps = Float64, dispersion = Float64, overprediction = Float64,
        underprediction = Float64, coverage_50 = Float64,
        coverage_90 = Float64,
        bias = Float64, n_samples = Int,
        log_rel_to_baseline = Float64))
frozen_overlay_df = _release_data("forecast_overlay_frozen.csv",
    (; release = String, made_date = Date, stream = String, horizon = Int,
        target_date = Date, fit = String, observed = Float64,
        median = Float64, lo30 = Float64, hi30 = Float64, lo60 = Float64,
        hi60 = Float64, lo90 = Float64, hi90 = Float64))
# The frozen evaluation never carries an individual single-stream fit
# (it scores only the joint model at past cut-offs), so the individual-fit
# comparison columns are dropped rather than shown as a column of missing.
frozen_score_overview_table = drop_individual_fit_columns(
    forecast_score_overview(frozen_scores_df))
frozen_score_by_horizon_table = drop_individual_fit_columns(
    forecast_score_by_horizon(frozen_scores_df))
frozen_score_by_release_table = drop_individual_fit_columns(
    forecast_score_by_release(frozen_scores_df))

# `fit` is single-valued (`FROZEN_FIT`) by construction in every one of
# these tables, not just for the releases scored so far, so it is dropped
# from the display tables below as a degenerate column. The `..._table`
# frames above keep it and still feed the relative-skill plots, which read
# it to colour each series in the joint role.
frozen_score_overview_display = drop_degenerate_fit_column(
    frozen_score_overview_table)
frozen_score_by_horizon_display = drop_degenerate_fit_column(
    frozen_score_by_horizon_table)
# See the comment above `joint_score_by_release_table`'s assignment for why
# this setup chunk's last statement needs a trailing `;`.
frozen_score_by_release_display = drop_degenerate_fit_column(
    frozen_score_by_release_table);

The headline frozen table pools one row per stream across cut-offs and horizons, scored against the persistence baseline on both scales, with the CRPS decomposition, coverage and bias columns described above. There is no model column, since only one model is scored.

streamncrpsrel_to_baselinelog_crpslog_rel_to_baselinedispersionoverpredictionunderpredictioncoverage_50coverage_90bias
confirmed cases163275.541.930.2080.62192.6582.620.270.730.870.36
confirmed deaths163251.443.092.1253.8364.454.93182.060.20.93-0.43
isolation beds14392.520.820.1240.5135.649.497.430.360.930.39
onset reports29173.731.040.5170.7872.772.3328.710.410.930.15
recovered63210.772.760.4742.14172.5834.413.780.860.890.13

The same relative skill against the baseline, by horizon, for the frozen cut-offs.

julia
frozen_relative_skill_fig = plot_forecast_relative_skill(
    frozen_score_by_horizon_table);

The same columns as a table, broken out by horizon, and again broken out by frozen cut-off, are behind the two dropdowns below.

Scores by horizon
streamhorizonncrpsrel_to_baselinelog_crpslog_rel_to_baselinedispersionoverpredictionunderpredictioncoverage_50coverage_90bias
confirmed cases780110.561.30.1720.3470.7239.290.560.70.850.35
confirmed cases1437260.643.320.2112.2172.3888.2600.70.920.38
confirmed cases2134510.12.680.2841.52343.42166.6800.760.820.38
confirmed cases2812756.831.290.2240.73640.88115.9500.9210.29
confirmed deaths780167.786.342.7183.6410.383.28154.110.20.85-0.46
confirmed deaths1437254.684.71.7987.92597.37188.310.191-0.39
confirmed deaths2134353.922.291.42.56131.987.89214.050.181-0.38
confirmed deaths2812508.871.571.2383.84250.460258.410.331-0.54
isolation beds74090.212.250.1131.1724.5764.770.880.20.750.62
isolation beds143798.121.920.1241.0336.6961.4300.4110.51
isolation beds213491.840.540.1270.3540.5442.428.880.4110.25
isolation beds283289.650.410.1370.3142.924.0822.670.4410.12
onset reports714124.70.690.8020.7351.2750.5522.880.430.930.09
onset reports1411200.161.260.2640.9280.3283.9735.860.2710.31
onset reports214272.691.840.2151.35126.71116.5129.470.750.75-0.05
recovered72094.643.70.4992.7745.546.652.490.70.80.25
recovered1417178.174.970.4823.56106.9967.383.80.820.820.18
recovered2114218.832.060.4391.64209.33.775.76110
recovered2812441.12.40.4611.28434.453.063.59110.01
Scores by frozen cut-off
made_datestreamncrpsrel_to_baselinelog_crpslog_rel_to_baselinedispersionoverpredictionunderpredictioncoverage_50coverage_90bias
2026-05-23confirmed cases2025.410.310.1580.2623.4601.9511-0.15
2026-05-27confirmed cases20141.150.780.9660.7713.190127.9600.15-0.93
2026-06-08confirmed cases60382.477.430.414.13172.85209.6200.350.830.65
2026-07-16confirmed cases4294.680.460.1680.29283.4610.980.24110.07
2026-07-18confirmed cases4361.260.780.1850.51333.8527.410110.15
2026-07-19confirmed cases4439.621.070.2090.7351.8387.790110.3
2026-07-20confirmed cases4385.620.850.1850.54332.653.020110.22
2026-07-24confirmed cases4422.921.350.1920.96343.779.220110.32
2026-07-25confirmed cases4460.722.050.2091.55356.82103.910110.34
2026-07-26confirmed cases4315.561.480.1641.29284.5930.970110.24
2026-07-27confirmed cases4346.391.710.1761.32306.9739.420110.27
2026-07-28confirmed cases4396.781.370.1560.86361.5235.260110.22
2026-07-31confirmed cases4498.383.850.2132.58369.58128.800.7510.38
2026-08-04confirmed cases4281.573.280.1563.03253.9927.590110.18
2026-08-08confirmed cases4253.513.180.152.33246.676.120.72110.03
2026-08-10confirmed cases3184.423.560.1462.72180.351.652.4211-0.03
2026-08-15confirmed cases3424.167.830.2624.82264.67159.490010.54
2026-08-17confirmed cases2223.283.040.2172.4479.72143.560010.75
2026-08-18confirmed cases2272.282.50.2832.0875.44196.84000.50.92
2026-08-23confirmed cases2293.456.870.2965.8467.34226.11000.50.88
2026-08-24confirmed cases1157.597.350.2446.2241.22116.370000.92
2026-08-25confirmed cases1124.074.570.1963.6435.3888.690010.82
2026-08-30confirmed cases151.441.560.081.3132.8918.550110.48
2026-05-23confirmed deaths209.420.380.4310.313.42060.31-0.56
2026-05-27confirmed deaths2031.180.961.4391.21.88029.300.3-0.94
2026-06-08confirmed deaths6048.520.710.2190.4941.955.790.780.910.03
2026-07-16confirmed deaths43371.11.6183.2483.740253.2501-0.75
2026-07-18confirmed deaths4284.841.311.5164.9111.230173.6101-0.66
2026-07-19confirmed deaths4274.731.341.435.21110.30164.4301-0.64
2026-07-20confirmed deaths4312.411.441.7235.75107.040205.3801-0.64
2026-07-24confirmed deaths4286.292.081.81911.47119.460166.830.51-0.63
2026-07-25confirmed deaths4303.322.521.8812.76135.850167.470.251-0.58
2026-07-26confirmed deaths4358.623.152.39718.65104.710253.9101-0.71
2026-07-27confirmed deaths4363.043.182.38717.17146.080216.960.251-0.58
2026-07-28confirmed deaths4436.852.752.13611.26185.320251.5301-0.67
2026-07-31confirmed deaths4351.183.112.02814.84172.540178.640.51-0.57
2026-08-04confirmed deaths4460.336.213.3438.2874.880385.4500.75-0.78
2026-08-08confirmed deaths4521.227.363.72334.8352.580468.6400.75-0.83
2026-08-10confirmed deaths3484.8212.714.42363.8720.260464.5700.67-0.88
2026-08-15confirmed deaths3272.827.212.36835.52150.30122.520.671-0.47
2026-08-17confirmed deaths235.610.720.0690.6224.2211.390110.38
2026-08-18confirmed deaths251.50.770.1140.7223.6827.830010.64
2026-08-23confirmed deaths268.341.240.1391.2725.1743.170010.69
2026-08-24confirmed deaths120.131.840.0631.7116.613.520110.29
2026-08-25confirmed deaths120.71.110.0650.9413.577.130110.4
2026-08-30confirmed deaths121.081.620.071.4113.27.880110.44
2026-06-08isolation beds8074.140.440.150.3838.172.2933.680.51-0.21
2026-07-16isolation beds473.222.640.0962.5320.0153.220010.6
2026-07-18isolation beds477.782.760.1032.6118.7159.070010.61
2026-07-19isolation beds4100.013.470.1323.3122.3877.630010.64
2026-07-20isolation beds482.233.20.1063.0519.8962.3300.2510.61
2026-07-24isolation beds476.612.770.0962.5910.1466.470010.73
2026-07-25isolation beds4106.432.710.1362.6213.8492.59000.750.73
2026-07-26isolation beds494.882.440.1222.3222.6172.2700.250.750.63
2026-07-27isolation beds485.562.90.1132.8138.9246.6500.510.48
2026-07-28isolation beds4153.082.770.2012.4841.15111.9400.50.750.67
2026-07-31isolation beds4119.063.040.1392.7841.3677.700.50.750.57
2026-08-04isolation beds4110.511.550.1531.4963.2647.2400.750.750.42
2026-08-08isolation beds459.170.890.0780.9151.547.150.480.7510.16
2026-08-10isolation beds346.180.870.0620.8745.960.020.2110
2026-08-15isolation beds3101.761.350.1071.1167.9933.760110.36
2026-08-17isolation beds291.773.550.13.0454.5137.2600.510.54
2026-08-18isolation beds275.081.310.0841.1348.6426.4400.510.44
2026-08-23isolation beds2128.83.940.1383.5447.5481.270010.77
2026-08-24isolation beds172.932.560.082.2531.9640.970010.64
2026-08-25isolation beds144.530.730.0540.7226.99017.5411-0.49
2026-08-30isolation beds131.631.180.0391.1931.5900.0411-0.01
2026-07-26onset reports294.331.630.1931.1965.2229.1100.510.43
2026-07-27onset reports2124.721.10.2050.856.6667.740.320.510.34
2026-07-28onset reports2148.651.240.2360.4661.3986.460.80.510.36
2026-07-31onset reports2100.670.970.1580.8875.7224.960110.33
2026-08-04onset reports3121.10.710.3730.867.940.9612.240.6710.25
2026-08-08onset reports3184.290.950.3080.7371.380112.910.331-0.65
2026-08-10onset reports3169.940.740.30.3966.10103.830.331-0.63
2026-08-15onset reports3390.243.270.3351.52135.96254.2800.330.670.65
2026-08-17onset reports2110.761.150.2190.8571.8234.344.590.510.14
2026-08-18onset reports2170.911.630.3251.3478.0137.2255.6801-0.09
2026-08-23onset reports2199.560.660.8940.8266.89132.680010.71
2026-08-24onset reports1125.580.50.6990.6962.2563.330010.59
2026-08-25onset reports1358.630.725.7071.0235.28323.350000.99
2026-08-30onset reports158.190.310.1660.1135.24022.9511-0.48
2026-07-16recovered4139.341.10.4691.09135.503.8411-0.1
2026-07-18recovered4149.612.290.4552.24146.942.670110.1
2026-07-19recovered4181.472.240.4822.11179.352.120110.07
2026-07-20recovered4171.812.130.4171.93169.632.180110.06
2026-07-24recovered4179.372.260.4091.93175.364.010110.13
2026-07-25recovered4177.72.040.3791.64173.524.180110.12
2026-07-26recovered4139.991.510.3811.57138.990111-0.05
2026-07-27recovered4219.872.190.4831.81218.3801.4911-0.06
2026-07-28recovered4544.073.750.3690.83543.250.50.33110.02
2026-07-31recovered4249.162.390.4021.48242.916.250110.12
2026-08-04recovered4157.982.310.4572.71146.41011.5711-0.19
2026-08-08recovered4125.591.620.4172.13109.72015.8711-0.23
2026-08-10recovered3119.751.970.5242.5885.84033.911-0.32
2026-08-15recovered3241.6413.110.4126.36212.8628.780110.28
2026-08-17recovered2309.2112.680.7736.6557.82251.390000.92
2026-08-18recovered2342.3115.630.8237.8452.64289.670000.94
2026-08-23recovered2303.4925.250.69812.859.73243.76000.50.9
2026-08-24recovered1187.3118.080.69810.5140.92146.390000.92
2026-08-25recovered1182.9616.990.679.8339.71143.250000.93
2026-08-30recovered1170.2517.020.60910.2635.53134.720010.86

The frozen forecasts made at each cut-off against the value observed since, one panel per stream and horizon, the observed value in black. Each panel carries the frozen forecast and the persistence baseline, coloured as in the cross-release overlay above, and the x-axis is the cut-off each forecast was made from. The same per-panel axis crop and overflow marker applies here.

Frozen-fit forecasts-versus-now overlay
julia
frozen_overlay_fig = plot_forecast_overlay(frozen_overlay_df);

The frozen re-fits below freeze the renewal data to an earlier cut-off and re-fit, so that a change driven by newer data can be distinguished from one driven by a change of method. Each uses the full headline settings: 1000 draws across two chains.

Freeze the renewal data to a cut-off and re-fit
julia
# Frozen re-fits and released_df are prepared in the setup block above.

Individual fits against the baseline ​

This section carries the same cross-release forecast scoring as Forecast scoring across releases above, for each stream's own individual fit rather than the joint, against the same persistence baseline. Recovered has no individual fit, so it does not appear here. These are the individual-fit rows of the same scored forecasts, not a separate computation. Only the newest few releases carry an individual-stream forecast, so these tables cover those releases alone rather than the outbreak's history.

Individual-fit rows of the cross-release scores
julia
# The relative skill against a stream's individual fit is only ever
# computed on the joint model's row, so on these rows it is missing by
# construction and the column is dropped rather than shown empty.
individual_score_overview_table = drop_individual_fit_columns(
    select_fit_role(forecast_score_overview_table, "individual"))
individual_score_by_horizon_table = drop_individual_fit_columns(
    select_fit_role(forecast_score_by_horizon_table, "individual"))
# See the comment above `joint_score_by_release_table`'s assignment for why
# this setup chunk's last statement needs a trailing `;`.
individual_score_by_release_table = drop_individual_fit_columns(
    select_fit_role(forecast_score_by_release_table, "individual"));
streamfitncrpsrel_to_baselinelog_crpslog_rel_to_baselinedispersionoverpredictionunderpredictioncoverage_50coverage_90bias
confirmed casesconfirmed53285.741.930.1931.882758.712.04110.02
confirmed deathsconfirmed_deaths53264.862.820.3242.59223.5835.365.920.9110.03
isolation bedstreatment5386.491.610.141.9373.260.1113.120.941-0.21
onset reportsonsets25229.421.220.5340.78109.877.7841.840.40.960.09

The same relative skill against the baseline, by horizon, one panel per stream (dataset), for each stream's own individual fit.

julia
individual_relative_skill_fig = plot_forecast_relative_skill(
    individual_score_by_horizon_table;
    empty_message = "Empty: no release old enough for its targets to " *
                    "have been observed carries an individual-stream " *
                    "forecast. Not a missing forecast.");

The same columns as a table, broken out by horizon, and again broken out by release, are behind the two dropdowns below.

Scores by horizon
streamhorizonfitncrpsrel_to_baselinelog_crpslog_rel_to_baselinedispersionoverpredictionunderpredictioncoverage_50coverage_90bias
confirmed cases7confirmed1797.182.480.1812.692.221.373.5911-0.04
confirmed cases14confirmed14213.222.920.1892.87208.531.413.2811-0.02
confirmed cases21confirmed12369.431.990.1961.67361.417.930.09110.07
confirmed cases28confirmed10607.411.540.2181.13575.0832.330110.14
confirmed deaths7confirmed_deaths1777.3430.3594.1160.791.9214.630.881-0.16
confirmed deaths14confirmed_deaths14169.772.790.2872.85154.3611.773.640.9310.01
confirmed deaths21confirmed_deaths12332.832.870.3052.25293.8937.781.170.9210.11
confirmed deaths28confirmed_deaths10635.22.770.341.61512.86122.3400.910.27
isolation beds7treatment1793.362.520.1553.0380.890.2712.1911-0.16
isolation beds14treatment1484.071.740.1382.0976.370.087.6211-0.15
isolation beds21treatment1283.51.210.1331.4368.08015.420.921-0.25
isolation beds28treatment1081.791.140.1251.3362.13019.660.81-0.37
onset reports7onsets11139.130.70.8160.6779.3240.3619.450.3610.16
onset reports14onsets8301.41.70.3681.26116.28119.3365.780.380.880.15
onset reports21onsets6298.971.660.2361.15157.0390.9850.960.51-0.11
Scores by release
made_datestreamfitncrpsrel_to_baselinelog_crpslog_rel_to_baselinedispersionoverpredictionunderpredictioncoverage_50coverage_90bias
2026-07-23confirmed casesconfirmed4327.810.840.2060.81319.434.843.5511-0.03
2026-07-25confirmed casesconfirmed4346.241.530.2131.57343.882.30.06110.02
2026-07-26confirmed casesconfirmed4299.071.360.2111.59291.390.167.5111-0.12
2026-07-27confirmed casesconfirmed4267.621.320.1711.3262.035.380.21110.04
2026-07-31confirmed casesconfirmed4402.872.880.2192.49383.8219.050110.09
2026-08-01confirmed casesconfirmed4336.312.870.2093322.9813.330110.08
2026-08-02confirmed casesconfirmed4296.982.260.1712.05286.156.724.1111-0.03
2026-08-03confirmed casesconfirmed4291.81.520.1751.2281.4210.080.3110.05
2026-08-04confirmed casesconfirmed4350.843.840.1913.45323.4127.40.04110.11
2026-08-07confirmed casesconfirmed4308.43.450.1922.75294.779.514.1211-0.01
2026-08-11confirmed casesconfirmed3235.245.570.1834.71220.4213.711.11110.06
2026-08-15confirmed casesconfirmed3216.854.130.2023.84213.7303.1211-0.06
2026-08-17confirmed casesconfirmed2144.461.650.2092.06136.6107.8411-0.15
2026-08-22confirmed casesconfirmed2159.722.910.1923.14157.831.890110.07
2026-08-24confirmed casesconfirmed188.923.750.1563.5882.496.430110.14
2026-08-25confirmed casesconfirmed190.863.540.1583.184.136.730110.2
2026-08-30confirmed casesconfirmed188.042.020.1381.7179.798.250110.21
2026-07-23confirmed deathsconfirmed_deaths4701.544.310.4482.27383.29318.2500.2510.56
2026-07-25confirmed deathsconfirmed_deaths4237.781.950.2341.58214.2823.50110.18
2026-07-26confirmed deathsconfirmed_deaths4256.462.20.2551.91231.6524.810110.13
2026-07-27confirmed deathsconfirmed_deaths4310.272.840.2521.88280.8429.430110.21
2026-07-31confirmed deathsconfirmed_deaths4312.362.630.2892.05292.2120.080.07110.08
2026-08-01confirmed deathsconfirmed_deaths4280.922.410.2531.73261.4619.460110.1
2026-08-02confirmed deathsconfirmed_deaths4256.052.830.2762.38238.2316.880.94110.03
2026-08-03confirmed deathsconfirmed_deaths4224.51.640.2451.23212.210.831.46110.02
2026-08-04confirmed deathsconfirmed_deaths4248.773.270.2833.16242.124.292.3511-0.02
2026-08-07confirmed deathsconfirmed_deaths4231.53.060.3252.97223.681.026.7911-0.11
2026-08-11confirmed deathsconfirmed_deaths3166.047.830.3298.55153.95012.0911-0.2
2026-08-15confirmed deathsconfirmed_deaths3169.784.080.3314.53163.2606.5211-0.14
2026-08-17confirmed deathsconfirmed_deaths2108.792.180.2712.4106.1902.611-0.12
2026-08-22confirmed deathsconfirmed_deaths2130.721.790.4833.36111.26019.4611-0.22
2026-08-24confirmed deathsconfirmed_deaths194.678.560.56615.2746.89047.7811-0.49
2026-08-25confirmed deathsconfirmed_deaths1106.136.560.79913.5251.28054.8501-0.55
2026-08-30confirmed deathsconfirmed_deaths1109.7610.960.90424.6545.13064.6301-0.55
2026-07-23isolation bedstreatment455.731.890.0892.347.7907.9411-0.21
2026-07-25isolation bedstreatment451.681.40.0821.6847.140.014.5211-0.14
2026-07-26isolation bedstreatment461.671.640.11.9854.820.236.6311-0.12
2026-07-27isolation bedstreatment459.472.050.0982.4754.6704.811-0.18
2026-07-31isolation bedstreatment482.761.950.1322.4368.17014.60.751-0.28
2026-08-01isolation bedstreatment488.621.380.1531.7873.94014.670.751-0.25
2026-08-02isolation bedstreatment484.371.650.1412.0369.7014.670.751-0.28
2026-08-03isolation bedstreatment479.931.880.1342.3271.9108.0111-0.21
2026-08-04isolation bedstreatment487.061.230.1411.3982.731.223.1111-0.03
2026-08-07isolation bedstreatment4104.981.810.1712.2685.76019.2111-0.29
2026-08-11isolation bedstreatment3103.680.630.1730.7289.95013.7311-0.22
2026-08-15isolation bedstreatment3123.881.620.1982.0395.9027.9811-0.31
2026-08-17isolation bedstreatment298.592.090.1582.5490.907.6911-0.18
2026-08-22isolation bedstreatment2144.283.570.224.54106.64037.6511-0.3
2026-08-24isolation bedstreatment1134.635.570.2167.18100.67033.9611-0.36
2026-08-25isolation bedstreatment1137.862.320.2052.8103.65034.2111-0.26
2026-08-30isolation bedstreatment1117.984.920.1725.8199.19018.7911-0.25
2026-08-02onset reportsonsets3111.070.590.1780.5596.84014.2311-0.2
2026-08-03onset reportsonsets3135.351.50.3121.08100.3235.0300.6710.36
2026-08-04onset reportsonsets3138.850.790.3690.7388.3416.7133.80.3310.01
2026-08-07onset reportsonsets3221.051.880.4672.3782.320138.7301-0.63
2026-08-11onset reportsonsets3237.90.770.4360.4885.910151.990.331-0.59
2026-08-15onset reportsonsets3486.374.020.3881.82205.04281.3300.3310.62
2026-08-17onset reportsonsets2374.543.810.441.82153.61220.9300.50.50.63
2026-08-22onset reportsonsets2271.950.910.5830.85104.12167.830010.69
2026-08-24onset reportsonsets1143.070.620.7640.7782.0960.980010.52
2026-08-25onset reportsonsets1195.350.423.4050.6488.57106.780010.79
2026-08-30onset reportsonsets1112.250.60.6760.4482.46029.7911-0.37

Outbreak size estimated by each data stream ​

Each data stream constrains the latent outbreak size differently. 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);
StreamLower 90%Lower 60%Lower 30%Upper 30%Upper 60%Upper 90%
exports39321555847451347159115755539941175
deaths (DRC)269844499562414106314149314331276
cases (DRC)3013637079419835523067232114167
confirmed (DRC)50324607726984995426117564205410
isolation (DRC)116141342714991190322347642249
onsets (DRC)1904227811359406101183812149746
joint97491120612399145311611119236

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.

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);

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);

Estimate evolution across releases ​

How the outbreak-size estimate has moved as situation reports accrued, three series on one calendar axis. The estimate published at each release is in blue, drawn as a median with nested 30/60/90% interval bars because each release is its own fit rather than one continuous model. The current model frozen at earlier cut-offs is in red, reusing fits already made for the McCabe and Chamla comparisons and the forecast validation. The current model on current data is the green band, drawn day by day so the latest estimate reads against the earlier points. Dotted vertical rules mark the release dates. The published series switches from a closed-form integral model to a renewal model on 7 June, so a step there can reflect the change of method rather than of data.

Released estimates and the current-model frozen re-fits
julia
# Released median and 30/60/90% intervals per release, from
# `data/released_estimates.csv`. Each tuple is
# `(date, median, lo30, hi30, lo60, hi60, lo90, hi90)`.
release_evolution = [(string(r.date), r.median, r.lo30, r.hi30, r.lo60, r.hi60,
                         r.lo90, r.hi90) for r in eachrow(released_df)]

# The current model frozen at earlier cut-offs, each its own discrete
# estimate: the matched-McCabe cut-offs (20, 23, 27 May) already computed
# for the matched-in-time comparison below, the 8 June Chamla
# confirmed-case anchor computed for the Chamla comparison, and the
# one-week-back validation fit (`frozen_lastweek`, at `validation_cutoff`)
# already computed for the forecast validation above. All are reused here so
# the current-model estimate at those earlier cut-offs reads against the
# released overlay, including a recent point one week before the cut-off.
# No extra fits are run. Each tuple carries the median and 30/60/90%
# credible bounds from the frozen draws; `round_fn` rounds to a whole count
# for outbreak size, and is passed through unrounded for a continuous
# quantity such as R0.
function _ci369(xs; round_fn = x -> round(Int, x))
    q(p) = round_fn(quantile(xs, p))
    (q(0.5), q(0.35), q(0.65), q(0.20), q(0.80), q(0.05), q(0.95))
end
frozen_by_cutoff[validation_cutoff] = frozen_lastweek
# The cut-offs every frozen fit above was made at, shared by the
# outbreak-size and R0 by-release overlays below.
_frozen_matched_cutoffs = sort(union(frozen_cutoffs,
    [validation_cutoff, default_chamla_cutoff()]))
frozen_matched = [(c, _ci369(frozen_C(c))...) for c in _frozen_matched_cutoffs]

# The current-data, current-model estimate as the cumulative-infection
# trajectory over the day grid (one calendar date per grid day, day 1 is
# the seeding date), summarised by per-day 30/60/90% credible bounds. This
# is the same latent quantity the cumulative-trajectory figure shows, so
# the current estimate rises over time on the release-date axis instead of
# sitting flat. Drawn against calendar dates, it lines up with the
# release and frozen points.
infection_trajectory = let
    mat = chn_joint[:cumulative_infections]
    trajs = [collect(v) for v in vec(collect(mat))]
    # Only over the comparison window — from the earliest release date to the
    # cut-off — not back to the seeding date.
    start_day = obs.n - value(obs.cutoff - Date(release_evolution[1][1]))
    days = max(start_day, 1):obs.n
    dates = [obs.seeding + Day(d - 1) for d in days]
    q(d, p) = quantile(Float64[t[d] for t in trajs], p)
    (dates,
        [q(d, 0.35) for d in days], [q(d, 0.65) for d in days],
        [q(d, 0.20) for d in days], [q(d, 0.80) for d in days],
        [q(d, 0.05) for d in days], [q(d, 0.95) for d in days])
end

evolution_fig = plot_estimate_evolution(release_evolution;
    renewal = frozen_matched,
    renewal_label = "Current model frozen at earlier cut-offs",
    trajectory = infection_trajectory,
    title = "Outbreak-size estimate as data accrued");

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);

Reproduction number by release ​

The reproduction number estimated at each release, the same kind of release-by-release picture as the outbreak-size evolution above. Each release's cut-off reproduction number is drawn as a discrete estimate, a median with nested 30/60/90% interval bars. The current fit's daily over its established window is drawn as the continuous band, and   is marked.

Reproduction number per release with the current-fit band
julia
rt_release_df = CSV.read(
    joinpath(pkgdir(BVDOutbreakSize), "data", "rt_by_release.csv"), DataFrame)
rt_release = [(string(r.date), r.median, r.lo30, r.hi30, r.lo60, r.hi60,
                  r.lo90, r.hi90) for r in eachrow(rt_release_df)]

# The current fit's daily Rt over its established window, summarised per day
# into a 30/60/90% band, reusing the same walk reconstruction the Rt figure
# uses so the band lines up with the per-release points on the calendar axis.
# The band is drawn only from the first release date onward, so it spans the
# same window as the per-release estimates rather than extending back to the
# renewal start. The first release day is the earliest date in
# `rt_by_release.csv` as a grid day; the walk is still reconstructed from the
# renewal start `_rt_start_plot` (the model knot grid) and the window is
# clamped into the reconstructed range so the quantiles never hit masked days.
rt_release_trajectory = let
    rt_walk_start = clamp(_BREAKPOINT - RT_WALK_LEAD, _rt_start_plot, obs.n)
    mat = reconstruct_rt(chn_joint; n = obs.n, breakpoint = _BREAKPOINT,
        rt_start = _rt_start_plot, rt_walk_start = rt_walk_start,
        ramp = RT_INTERVENTION_RAMP)
    first_release_day = clamp(
        value(minimum(rt_release_df.date) - obs.seeding) + 1,
        _rt_start_plot, obs.n)
    days = first_release_day:obs.n
    dates = [obs.seeding + Day(d - 1) for d in days]
    q(d, p) = quantile(collect(skipmissing(@view mat[:, d])), p)
    (dates,
        [q(d, 0.35) for d in days], [q(d, 0.65) for d in days],
        [q(d, 0.20) for d in days], [q(d, 0.80) for d in days],
        [q(d, 0.05) for d in days], [q(d, 0.95) for d in days])
end

rt_evolution_fig = plot_estimate_evolution(rt_release;
    trajectory = rt_release_trajectory,
    ylabel = "Reproduction number",
    title = "Reproduction number as data accrued",
    released_label = "Released estimate (per project release)",
    trajectory_label = "Current model, current data",
    refline = 1.0);

Reproduction number by release and dataset ​

The same release-by-release reproduction number split into one panel per dataset, so each dataset's history reads against the others and against the joint. Panels share a calendar axis and a y range, and   is marked. Each release's cut-off value is a median with nested 30/60/90% interval bars. A dataset the report also fits on its own carries that fit's current-model band behind its points, built as in the overview above. Confirmed deaths carries no band, so its panel shows release points alone. Only the most recent releases published these per-dataset estimates, so every panel spans a much shorter window than the overview above rather than a different history.

Reproduction number per release by fit
julia
# Schema of the per-release, per-fit estimate tables written by
# scripts/score_releases.jl from each release's stream_estimates.csv.
_by_stream_schema = (; release = String, date = Date, fit = String,
    median = Float64, lo30 = Float64, hi30 = Float64, lo60 = Float64,
    hi60 = Float64, lo90 = Float64, hi90 = Float64)

# Fits in a fixed order, the joint first, so the panels do not reshuffle
# between builds. Labels match the per-stream table above (in "Outbreak
# size estimated by each data stream"). Recovered is absent because it
# has no individual fit.
_fit_order = ["joint", "cases", "deaths", "confirmed", "confirmed_deaths",
    "treatment", "onsets", "exports"]
_fit_labels = Dict("joint" => "joint", "cases" => "cases (DRC)",
    "deaths" => "deaths (DRC)", "confirmed" => "confirmed (DRC)",
    "confirmed_deaths" => "confirmed deaths (DRC)",
    "treatment" => "isolation (DRC)", "onsets" => "onsets (DRC)",
    "exports" => "exports")

# Group a per-fit estimate table into the label => tuples pairs the faceted
# plot takes, keyed on the date so the mixed release tag shapes
# (`results-v1.9.0` and `results-1243`) never reach the axis.
function _fit_groups(df)
    return [get(_fit_labels, f, f) =>
                [(string(r.date), r.median, r.lo30, r.hi30, r.lo60, r.hi60,
                     r.lo90, r.hi90) for r in eachrow(df) if r.fit == f]
            for f in _fit_order]
end

# Per-fit reproduction-number trajectory, reconstructing the walk exactly as
# `plot_rt_streams` does per stream.
function _stream_rt_trajectory(chn, dates; rt_start, rt_walk_start)
    mat = reconstruct_rt(chn; n = obs.n, breakpoint = _BREAKPOINT,
        rt_start = rt_start, rt_walk_start = rt_walk_start,
        ramp = RT_INTERVENTION_RAMP)
    first_date = isempty(dates) ? obs.seeding : minimum(dates)
    first_day = clamp(value(first_date - obs.seeding) + 1, rt_start, obs.n)
    days = first_day:obs.n
    ds = [obs.seeding + Day(d - 1) for d in days]
    q(d, p) = quantile(collect(skipmissing(@view mat[:, d])), p)
    (ds,
        [q(d, 0.35) for d in days], [q(d, 0.65) for d in days],
        [q(d, 0.20) for d in days], [q(d, 0.80) for d in days],
        [q(d, 0.05) for d in days], [q(d, 0.95) for d in days])
end

# The single-stream chains and their renewal-walk starts, keyed on the fit
# id the per-release tables use. Both the joint walk start and the day-1
# per-stream starts are the ones the per-stream implied-Rt figure above
# uses, so the bands here match it. Confirmed deaths has no trajectory
# here: its panel still draws its release points alone.
_stream_chains = (
    "joint" => (; chn = chn_joint, rt_start = _rt_start_plot,
        rt_walk_start = _rt_walk_start_joint),
    "cases" => (; chn = chn_cases, rt_start = 1, rt_walk_start = 1),
    "deaths" => (; chn = chn_deaths, rt_start = 1, rt_walk_start = 1),
    "confirmed" => (; chn = chn_confirmed, rt_start = 1, rt_walk_start = 1),
    "treatment" => (; chn = chn_treatment, rt_start = 1, rt_walk_start = 1),
    "onsets" => (; chn = chn_onsets, rt_start = 1, rt_walk_start = 1),
    "exports" => (; chn = chn_exports, rt_start = 1, rt_walk_start = 1))

# Build a fit label => trajectory dictionary from a per-release table,
# restricted to the fits `_stream_chains` names. A fit with no row in `df`
# gets no trajectory, so its panel still draws its release points alone.
function _rt_trajectories(df)
    trajs = Dict{String, Any}()
    for (fid, cfg) in _stream_chains
        fdates = df.date[df.fit .== fid]
        isempty(fdates) && continue
        trajs[get(_fit_labels, fid, fid)] = _stream_rt_trajectory(
            cfg.chn, fdates; rt_start = cfg.rt_start,
            rt_walk_start = cfg.rt_walk_start)
    end
    return trajs
end

rt_stream_df = _release_data("rt_by_release_by_stream.csv",
    _by_stream_schema)
rt_stream_fig = plot_evolution_by_group(_fit_groups(rt_stream_df);
    trajectories = _rt_trajectories(rt_stream_df),
    ylabel = "Reproduction number",
    title = "Reproduction number as data accrued, by dataset",
    released_label = "Released estimate (per release)",
    refline = 1.0,
    empty_note = "No per-dataset reproduction numbers saved yet.");

Basic reproduction number by release ​

The basic reproduction number estimated at each release, the initial-transmission counterpart of the reproduction number above, before the time-varying decline. Released estimates are blue and the current model frozen at earlier cut-offs is red, each a median with nested 30/60/90% interval bars. The current fit sits behind both as a flat band, and   is marked. Releases only began publishing this quantity recently, so the short blue history reflects that rather than any failed release. The frozen series carries the comparison meanwhile.

Basic reproduction number per release with frozen re-fits and the current-fit band
julia
# Per-release R0 points from r0_by_release.csv, read through the typed
# fallback so a missing or header-only file (until a release carries
# `rt_state.log_R0` in its posterior draws) does not break the build. The
# schema mirrors rt_by_release.csv.
_r0_schema = (; release = String, date = Date, median = Float64,
    lo30 = Float64, hi30 = Float64, lo60 = Float64, hi60 = Float64,
    lo90 = Float64, hi90 = Float64)
r0_release_df = _release_data("r0_by_release.csv", _r0_schema)
r0_release = [(string(r.date), r.median, r.lo30, r.hi30, r.lo60, r.hi60,
                  r.lo90, r.hi90) for r in eachrow(r0_release_df)]

# The current model frozen at earlier cut-offs, one discrete estimate per
# cut-off, reusing the same frozen fits `frozen_matched` above already
# computed. No extra fits are run. Each tuple carries the median and
# 30/60/90% credible bounds of that frozen fit's own R0 draws, unrounded
# since R0 is continuous.
frozen_r0_matched = [(c, _ci369(frozen_R0(c); round_fn = identity)...)
                     for c in _frozen_matched_cutoffs]

# The current fit's R0 posterior is a single distribution rather than a
# daily series, so it summarises into a flat 30/60/90% reference band. The
# window runs from the earliest mark on the axis, the first frozen cut-off
# or release point, to the current cut-off, so the band reads behind both
# series rather than only their recent end.
r0_reference = let
    draws = r0_walk_draws(chn_joint)
    q(p) = quantile(draws, p)
    first_date = min(minimum(Date.(_frozen_matched_cutoffs)),
        isempty(r0_release_df.date) ? obs.cutoff :
        minimum(r0_release_df.date))
    dates = [first_date, obs.cutoff]
    (dates, fill(q(0.35), 2), fill(q(0.65), 2), fill(q(0.20), 2),
        fill(q(0.80), 2), fill(q(0.05), 2), fill(q(0.95), 2))
end

r0_evolution_fig = plot_estimate_evolution(r0_release;
    renewal = frozen_r0_matched,
    renewal_label = "Current model frozen at earlier cut-offs",
    trajectory = r0_reference,
    ylabel = "Basic reproduction number",
    title = "Basic reproduction number as data accrued",
    released_label = "Released estimate (per project release)",
    trajectory_label = "Current model, current data",
    refline = 1.0);

Basic reproduction number by release and dataset ​

The basic reproduction number estimated at each release, one panel per fit, the by-dataset counterpart of the figure above. Panels share a calendar axis and a y range, and   is marked. Each release is a median with nested 30/60/90% interval bars. Every fit the report runs on its own also carries a current-model reference band. Panels fill in from the first release that publishes this quantity per dataset, so a fit with nothing saved yet is left out rather than drawn empty.

Basic reproduction number per release by fit
julia
# Per-fit R0 flat reference band, the by-dataset counterpart of
# `r0_reference` above, a single distribution rather than a daily walk, so
# each fit's band is flat across its own release window. `r0_walk_draws`
# probes for the walk base, so a single-stream model built without its own
# renewal walk drops its band instead of erroring.
function _r0_stream_trajectory(chn, dates)
    draws = r0_walk_draws(chn)
    isnothing(draws) && return nothing
    q(p) = quantile(draws, p)
    first_date = isempty(dates) ? obs.seeding : minimum(dates)
    ds = [first_date, obs.cutoff]
    (ds, fill(q(0.35), 2), fill(q(0.65), 2), fill(q(0.20), 2),
        fill(q(0.80), 2), fill(q(0.05), 2), fill(q(0.95), 2))
end

# Build a fit label => trajectory dictionary from a per-release R0 table,
# restricted to the fits `_stream_chains` names, the same restriction the
# reproduction-number-by-dataset trajectories use. A fit with no row in
# `df`, or whose chain carries no walk base, gets no trajectory, so its
# panel still draws its release points alone.
function _r0_trajectories(df)
    trajs = Dict{String, Any}()
    for (fid, cfg) in _stream_chains
        fdates = df.date[df.fit .== fid]
        isempty(fdates) && continue
        traj = _r0_stream_trajectory(cfg.chn, fdates)
        isnothing(traj) || (trajs[get(_fit_labels, fid, fid)] = traj)
    end
    return trajs
end

r0_stream_df = _release_data("r0_by_release_by_stream.csv",
    _by_stream_schema)
r0_stream_fig = plot_evolution_by_group(_fit_groups(r0_stream_df);
    trajectories = _r0_trajectories(r0_stream_df),
    ylabel = "Basic reproduction number",
    title = "Basic reproduction number as data accrued, by dataset",
    released_label = "Released estimate (per release)",
    refline = 1.0,
    empty_note = "No per-dataset basic reproduction numbers saved yet.");

Comparison with McCabe et al. ​

Our model is a discrete-time renewal model with a time-varying reproduction number and every data stream fitted jointly. McCabe et al. published their estimates as scenarios at fixed situation-report cut-offs, each scenario carrying a 95% confidence interval. We show all three, the 18 May report, the 20 May update and the 27 May Lancet publication, as one panel each, with their intervals kept. Within a panel each method and scenario family is a single line, carrying its sweep over the nuisance assumptions: the case-fatality ratio, the geographic window and the doubling time. The geographic-spread scenarios come from exported cases and travel volume. Their back-calculation-from-deaths scenarios differ between the reports, since the 18 May report used 88 reported deaths and the 20 May update 131. The 20 May update also corrected the case-fatality ratios. McCabe's scenarios estimate cumulative cases at their report dates, though their report is not fully explicit about whether this is symptomatic cases or all infections. We take the like-for-like quantity to be our cumulative symptom onsets on the same dates, not the latent infections (which include the not-yet-symptomatic) or our current cut-off total. We read our value off the joint fit's cumulative-onset trajectory at the grid day for each report date, and show it with its credible interval. Each scenario sits beside our estimate for the date it was made: the 18 May report against our 18 May value, the 20 May update against our 20 May value, and the 27 May Lancet publication against our 27 May value.

McCabe scenarios with uncertainty against our estimates
julia
function _ci90row(xs)
    (round(Int, quantile(xs, 0.5)),
        round(Int, quantile(xs, 0.05)),
        round(Int, quantile(xs, 0.95)))
end

# Our modelled cumulative symptom onsets on a McCabe report date, read off
# the joint fit's per-draw `cumulative_onsets` trajectory. The grid runs to
# the cut-off on day `n`, so the day-index for a date is `n` minus the days
# from that date back to the cut-off (`grid_day("2026-06-07") = n`,
# `"2026-05-20") = n - 18`, `"2026-05-18") = n - 20`).
_onset_trajs = let mat = chn_joint[:cumulative_onsets]
    [collect(v) for v in vec(collect(mat))]
end
# Inverse of `grid_date(day) = obs.cutoff - Day(obs.n - day)`: the day-index
# whose calendar date is `date`, using `value` (imported above) for the
# day count rather than the non-exported `Dates.date2epochdays`.
_grid_day(date) = obs.n - value(obs.cutoff - Date(date))
function _ours_on(date)
    d = _grid_day(date)
    _ci90row(Float64[t[d] for t in _onset_trajs])
end

# Our matched cumulative-onset estimate for each report date, keyed by date so
# it lands beside that vintage's scenarios in its own panel.
mccabe_ours = Dict(
    "2026-05-18" => _ours_on("2026-05-18"),
    "2026-05-20" => _ours_on("2026-05-20"),
    "2026-05-27" => _ours_on("2026-05-27"))

# One panel per report date; within a panel each method-and-family is one row,
# with the case-fatality / window / doubling-time sweep dodged onto that single
# line, so the ~40 scenarios keep their intervals without becoming ~40 rows.
matched_comparison_fig = plot_scenario_comparison(REPORT_SCENARIOS_CI;
    ours = mccabe_ours,
    date_titles = ["2026-05-18" => "18 May report",
        "2026-05-20" => "20 May update",
        "2026-05-27" => "27 May (Lancet)"],
    xlabel = "Cumulative cases");

The McCabe scenarios are outbreak-size estimates, the same quantity our renewal model and the released integral model report. Their 95% confidence intervals come from exact negative-binomial counts for the geographic-spread method and a Poisson likelihood profile for the back-calculation from deaths.

Frozen-fit C_T intervals (kept for the CSV export, not shown)
julia
# The estimate-evolution figure above already shows how the size estimate
# shifts as data accrues, so the side-by-side frozen-fit table is no longer
# rendered in the report; it is kept only to populate the published
# `frozen_matched_cutoffs.csv` export.
frozen_streams_table = streams_table(
    "frozen 20 May" => frozen_C("2026-05-20"),
    "frozen 23 May" => frozen_C("2026-05-23"),
    "frozen 27 May" => frozen_C("2026-05-27"),
    "frozen 8 June" => frozen_C(default_chamla_cutoff()),
    "current data" => posterior_C_joint);

Comparison with Chamla et al. ​

A second group, Chamla et al. (Chamla et al., 2026) at the World Health Organization Regional Office for Africa, published a stochastic compartmental model of the same outbreak on 25 June 2026. Their model is a discrete-time susceptible-exposed-infectious-recovered-dead ensemble, recalibrated by simulation filtering to the laboratory-confirmed case series and anchored on the 598 confirmed cases reported by 8 June. It is then run forward to project the confirmed-case trajectory under a low, central and high transmissibility scenario.

Their published quantity is the cumulative confirmed-case count, with the reporting fraction held at one, so it does not adjust for the cases that are infected but never laboratory-confirmed. This is a different quantity from the cumulative cases this analysis and McCabe et al. estimate, which include the unconfirmed and unascertained. It therefore sits below them: a floor on the true size rather than an estimate of it. The like-for-like comparison is therefore against our own confirmed-case projection, not against our cumulative infection count.

We compare forward projections rather than refitting to their assumptions. We take our fit frozen at 8 June, the exact date of their confirmed-case calibration anchor. We roll its confirmed-case stream forward to the dates Chamla report, using the same machinery as the one-week-ahead forecast. Setting our projection, their projection and the confirmed cases observed since on one timeline shows how each projection has held up against the data.

Project the 8 June fit forward and assemble the Chamla comparison
julia
# The 8 June frozen joint fit matches Chamla's confirmed-case calibration
# anchor exactly and carries the confirmed-case testing history through then,
# so we roll its confirmed-case stream forward with the one-week-ahead forecast
# machinery to the dates Chamla report.
chamla_anchor = frozen_by_cutoff["2026-06-08"]

# Our projected cumulative confirmed cases at a horizon of `h` days past the
# 8 June cut-off: a forward `forecast_reported` run (its reproduction number
# left to keep evolving), summarised as (median, 5%, 95%).
function _our_confirmed_h(h)
    fc = forecast_reported(chamla_anchor.chn;
        horizon = h,
        obs_cases = chamla_anchor.o.reported_cases,
        obs_deaths = chamla_anchor.o.total_deaths,
        obs_confirmed = chamla_anchor.o.confirmed_cases,
        obs_confirmed_deaths = chamla_anchor.o.confirmed_deaths)
    return _ci90row(float.(fc.confirmed_cum))
end

# Our projection at Chamla's forward report dates (10 and 24 June, week 12):
# the anchor day is the fitted confirmed total at 8 June, each later date a
# forward forecast. Reused for the matched-date table and the week-12 figure.
chamla_fan = map(["2026-06-08", "2026-06-10", "2026-06-24"]) do d
    h = value(Date(d) - chamla_anchor.cutoff)
    row = h == 0 ?
          (chamla_anchor.o.confirmed_cases, chamla_anchor.o.confirmed_cases,
        chamla_anchor.o.confirmed_cases) : _our_confirmed_h(h)
    (d, row...)
end
_fan_at(date) =
    let r = first(x for x in chamla_fan if x[1] == date)
        (r[2], r[3], r[4])
    end
ours_10jun = _fan_at("2026-06-10")
ours_24jun = _fan_at("2026-06-24")

# Observed confirmed cases over the comparison window: the daily cumulative
# series read off the chain's grid from 18 May (Chamla's first projected point)
# to the cut-off.
chamla_obs_series = let
    ds = [grid_date(d) for d in obs.confirmed_history.days]
    cs = obs.confirmed_history.counts
    [(string(ds[i]), cs[i]) for i in eachindex(ds) if ds[i] >= Date("2026-05-18")]
end

# Chamla's central confirmed-case projection over the comparison window; their
# later, far-larger horizons are noted in the text rather than plotted so the
# window stays legible.
chamla_central_window = CHAMLA_CONFIRMED_CENTRAL[1:4]

chamla_projection_fig = plot_projection_comparison(;
    external = chamla_central_window,
    ours = chamla_fan,
    observed = chamla_obs_series,
    external_label = "Chamla et al. central (R₀=1.71)",
    ours_label = "Our projection (from 8 June)",
    observed_label = "Observed confirmed",
    title = "Confirmed-case projections versus observed, from mid-May");

By 24 June their central scenario projected just under a thousand confirmed cases, and their low and high scenarios ranged from roughly 870 to 1360. The figure below sets that week-12 scenario spread beside our 8 June projection for the same date and the confirmed count observed by the cut-off, so each reads against their three scenarios at a glance.

Week-12 (24 June) scenario spread against ours and observed
julia
chamla_w12_rows = vcat(
    [(label, m, lo, hi) for (label, m, lo, hi) in CHAMLA_CONFIRMED_W12],
    [("Our projection (from 8 June)", ours_24jun...)],
    [("Observed by 23 June cut-off", obs.confirmed_cases,
        obs.confirmed_cases, obs.confirmed_cases)])
chamla_w12_groups = vcat(fill("Chamla et al. scenarios", 3),
    ["Our projection"], ["Observed"])

chamla_w12_fig = plot_estimate_comparison(chamla_w12_rows;
    xlabel = "Cumulative confirmed cases by 24 June",
    groups = chamla_w12_groups,
    group_colours = ["Chamla et al. scenarios" => :steelblue,
        "Our projection" => :firebrick,
        "Observed" => :black]);

The matched-date numbers behind these figures are in the dropdown below, with the observed column taken to the 23 June cut-off.

Matched-date projection numbers (10 and 24 June)
julia
chamla_comparison_table = let
    fmt(t) = string(t[1], " (", t[2], "–", t[3], ")")
    central(date) =
        let r = first(x for x in CHAMLA_CONFIRMED_CENTRAL
            if x[1] == date)
            fmt((r[2], r[3], r[4]))
        end
    DataFrame(
        "Date" => ["10 June", "24 June"],
        "Chamla central (90% PI)" => [central("2026-06-10"),
            central("2026-06-24")],
        "Our projection (90% CrI)" => [fmt(ours_10jun), fmt(ours_24jun)],
        "Observed confirmed" => [
            string(freeze_observations("2026-06-10").confirmed_cases),
            string(obs.confirmed_cases) * " (23 June)"])
end;
DateChamla central (90% PI)Our projection (90% CrI)Observed confirmed
10 June648 (470–812)679 (619–800)676
24 June990 (709–1293)1805 (1168–3402)6686 (23 June)

Beyond the comparison window their central scenario continues to roughly 8200 confirmed cases by mid-September, with the high scenario far higher. Those longer projections are not set against data here.

Reproduction number behind the projection ​

The forward projection above is carried by the reproduction-number trajectory our 8 June fit estimated, a quantity we report in its own right rather than as a comparison. The figure shows that trajectory, the time-varying reproduction number from the renewal walk with its credible intervals, as the fit saw it at 8 June. It declines over the weeks leading to the cut-off, and that decline is what bends the projected trajectory away from sustained early growth.

Reproduction number as estimated by the 8 June fit
julia
# Reconstruct the reproduction-number trajectory the 8 June fit estimated,
# mirroring the current-data R_t figure but with the frozen vintage's own grid,
# breakpoint and renewal start.
chamla_rt_obs = chamla_anchor.o
chamla_rt_breakpoint = chamla_rt_obs.n - chamla_rt_obs.who_first_sitrep_days
chamla_rt_start = clamp(
    chamla_rt_obs.n - round(Int, chamla_rt_obs.tmrca_days) + RENEWAL_START_LEAD,
    1, chamla_rt_obs.n)
chamla_rt_fig = plot_rt(chamla_anchor.chn;
    n = chamla_rt_obs.n, breakpoint = chamla_rt_breakpoint,
    rt_start = chamla_rt_start,
    rt_walk_start = clamp(chamla_rt_breakpoint - RT_WALK_LEAD,
        chamla_rt_start, chamla_rt_obs.n),
    as_of_date = string(chamla_rt_obs.cutoff),
    seeding = chamla_rt_obs.seeding, ramp = RT_INTERVENTION_RAMP);

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 uses the full headline settings: 1000 draws across two chains.

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._";
StreamLower 90%Lower 60%Lower 30%Upper 30%Upper 60%Upper 90%
baseline (hospital pathway)97491120612399145311611119236
community pathway101421145112468148291637220370
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.

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._";
StreamLower 90%Lower 60%Lower 30%Upper 30%Upper 60%Upper 90%
Skygrid (baseline)97491120612399145311611119236
Exponential growth99681117212214142991566518503
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._";
StreamLower 90%Lower 60%Lower 30%Upper 30%Upper 60%Upper 90%
Skygrid (baseline)163166170176181193
Exponential growth168171173180185196
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._";

Saving sensitivity results ​

The stream-comparison and frozen-fit tables and the per-stream reproduction number figure are written to the shared output directory. The main analysis writes the rest, so the combined release and summary dashboard pick up both pages' outputs.

Write sensitivity 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)
CSV.write(joinpath(output_dir, "frozen_matched_cutoffs.csv"),
    frozen_streams_table)

# The one-week-back validation forecast, in the same archive format as the
# release forecast, so the frozen "last week versus now" forecast is recorded
# as a release asset alongside the forecast it is scored against.
CSV.write(joinpath(output_dir, "forecast_validation.csv"),
    forecast_archive([(7, validation_forecast)];
        made_date = frozen_lastweek.o.cutoff, thin = 5))

# The per-stream reproduction-number figure for the summary dashboard; the
# main analysis writes the other three dashboard figures.
dashboard_dir = joinpath(
    pkgdir(BVDOutbreakSize), "docs", "src", "summary_assets")
mkpath(dashboard_dir)
CairoMakie.save(joinpath(dashboard_dir, "rt_streams.png"), stream_rt_fig)