Estimating the current size of the 2026 DRC Bundibugyo virus outbreak

Authors: Sam Abbott, Kath Sherratt, Samuel Brand and Sebastian Funk.

Last updated: 24 June 2026. This is a live report, re-run as new data arrive, so the estimates change between updates.

Data as of: 21 June 2026. DRC counts come from the situation reports of the Institut National de Santé Publique (INSP); Uganda imports come from WHO. The rendered report fills in the build date and the exact data cut-off automatically.

See: current outbreak size · one-week-ahead forecast · estimate evolution across releases · comparison with McCabe et al. · how the data streams compare · limitations · full joint results.

Abstract. An outbreak of Ebola disease caused by Bundibugyo virus (BVD) is ongoing in the Democratic Republic of the Congo (DRC), with cases also detected across the border in Uganda. This is a real-time joint Bayesian estimate of the current size of that outbreak, refreshed as new data arrive. Most infections are not yet reported, so the current size has to be inferred from the surveillance data that are available. The model is a discrete-time renewal process on a daily grid that fits the surveillance streams jointly in a single posterior: the DRC suspected cases, suspected deaths, laboratory-confirmed cases and confirmed deaths, and the cases and deaths exported to Uganda. It estimates the latent infections, symptom onsets and deaths over time, the reported and confirmed cases, and the time-varying reproduction number with its growth rate and doubling time, alongside the case-fatality ratio, the ascertainment of each surveillance system, and a short-term forecast of each stream over the coming week. The DRC data come from the INSP situation reports and the Uganda exports from the WHO situation reports and Disease Outbreak News, with a genetic bound on the time to the most recent common ancestor and priors taken from the McCabe et al. report.

Scope. This work is motivated by adding an external view of the current situation, based on our understanding of real-time infectious disease dynamics and the infection process that gives rise to observed epidemic surveillance counts. We are actively developing it and encourage feedback, so please get in touch. We fully support reuse and adaptation. Find out more in the contributing guide.

Use of AI. The model code and analysis were drafted by a language model and reviewed and revised under human oversight; the named authors are responsible for that oversight.

This page is generated from docs/examples/analysis.jl; the model code it calls is in src/. See the LLM-driven reimplementation limitation below for the oversight context behind the Use of AI note.

Offline copy. A self-contained single-file HTML version of this report, built from the same run, is attached to each results release: download the latest.

Origins of this work

This work began as a replication of the McCabe et al. (McCabe and others, May 2026) report. It has since evolved into a real-time joint Bayesian estimate of the current outbreak size, a discrete-time renewal process with a time-varying reproduction number fitted to more of the available data streams than the original. The points below summarise how it now differs from the report; the Methods section carries the full treatment, and the later comparison with McCabe et al. sets the current estimates against theirs.

Expand: differences from the report

Latent process and parameters

• Discrete-time renewal model. The whole model runs on a daily grid. Infections follow the discrete renewal equation $I_t=R_t\sum _{s\ge 1}{I}_{t-s}{g}_s$, where $g$ is the discretised generation-interval PMF, and every delay is applied as a discrete convolution. McCabe et al. (McCabe and others, May 2026) use continuous-time closed forms.

• Time-varying reproduction number. $R_t$ is held flat at the established $R_0$ until the first WHO situation report (18 May 2026), then follows a weekly Gaussian random walk on the log scale, interpolated within weeks, with a logistic outbreak-response ramp of about three weeks from that report. McCabe et al. use one constant exponential growth rate.

• Joint posterior rather than scenario estimates. The reproduction number, case-fatality ratio, all delays, traveller volume and surveillance dispersion have priors and are sampled together. McCabe et al. (McCabe and others, May 2026) fix each and report a set of scenarios.

• Two-phase seeding with a wide, genetically-floored outbreak age. A single import grows through an unobserved cryptic exponential phase to a magnitude set by a wide prior on the doubling count, at the growth rate the genetic estimate informs, before the renewal process takes over. The established reproduction number is derived forward from that growth rate. The genetic time to the most recent common ancestor floors the cryptic duration from below. McCabe et al. fix the start from a single seed.

Delays and convolutions

• Delays re-estimated with uncertainty. McCabe et al. (McCabe and others, May 2026) take the onset-to-death delay from the Isiro 2012 point estimate of Rosello et al. (Rosello and others, 2015). We instead use a Bayesian reanalysis of the same line list (Funk and Abbott, 2026) that re-estimates the delay with uncertainty, and we sample every other delay (generation interval, incubation period, onset-to-report, onset-to-confirmation and onset-to-detection abroad) from a prior centred on published Ebola estimates, discretised with double interval censoring (Charniga et al., 2024), so the delay uncertainty propagates.

Likelihoods and data streams

• More streams fitted. McCabe et al. (McCabe and others, May 2026) fit the Uganda export cases and deaths. We add the DRC suspected cases, the laboratory-confirmed cases, the confirmed deaths and the deaths among the Uganda exports.

• Per-vintage time-series fitting. The DRC streams are fitted on the incidence scale, as the between-vintage increments across successive sitreps (the first vintage being the cumulative count to that date), which sharpens $R_t$. McCabe et al. condition on a single cumulative total.

• Ascertainment estimated. We jointly estimate the outbreak size and the fraction of cases each surveillance system reports. McCabe et al. have no ascertainment component.

• Comparison against published scenarios. The model is set beside the McCabe et al. (McCabe and others, May 2026) scenario estimates as an external sense-check, matched in time at the cut-off each scenario was computed, while the cumulative infection count, the running sum of the daily infections, is the headline quantity reported separately.

Extensions

• No-onward-transmission counterfactual and one-week-ahead forecasts. Future expected deaths from infections already seeded, and a posterior-predictive projection of each stream.

Limitations

The limitations are grouped by the data, the model assumptions and design, and the implementation, with the most consequential first in each group.

Expand: limitations

Data and what it can support

• Most quantities rest on weakly-informed priors. Nearly all of the delays, the case-fatality ratio and the laboratory assumptions are set by priors informed at best by a handful of literature sources, often from other outbreaks, and in places by our own prior judgement rather than anything from this outbreak. The data do little to move them, so these posteriors largely track their priors. We fit the between-report increments, so the trajectory informs the change in the reproduction number over the window but is uninformative about the delays, the surveillance dispersion or the reporting fractions on their own.

• Only report-date totals, no epidemiological dating. We have no counts by symptom onset or any other epidemiologically relevant date, only cumulative totals at the report date. The timing of the underlying epidemic is therefore weakly identified, and we recover it only through the assumed delays.

• Fitted to aggregate counts. The DRC data are situation-report totals of suspected cases and deaths, laboratory-confirmed cases and deaths, and specimens received and analysed; the Uganda data are three export cases with one death. We do not have a line list or information on case definitions or reporting completeness. The laboratory testing series gives partial information on testing capacity, but it is incomplete and stops at the cut-off. Every estimate is a model-based extrapolation under strong assumptions, not a measurement.

• Later sitreps revise earlier figures. A later situation report can revise an earlier total up or down as suspects are reclassified and newly-reporting health zones are added, and ascertainment probably rose over the window. We do not model this revision process.

• Streams share one case pool. They are fitted as conditionally independent given latent incidence but observe overlapping people, which can understate uncertainty. Whether the streams imply mutually consistent outbreak sizes is not assessed here.

Model assumptions and design

• Inherits McCabe et al.'s epidemiological assumptions. A single zoonotic seed, an assumed generation interval, no spatial structure beyond the Ituri / Nord Kivu split, and no depletion of susceptibles. The onset-to-death delay is grounded on Isiro 2012 and the genetic seeding bound on an external clock rate, neither propagating cross-outbreak or clock uncertainty.

• Intervention ramp is weakly identified. With only a few sitreps straddling it, the ramp effect and the pre-ramp reproduction number are not well separated.

Implementation

• LLM-driven reimplementation. The model code, priors and analysis were drafted by a language model from the McCabe et al. (McCabe and others, May 2026) report and the companion delay reanalysis, then reviewed and revised. Not independently replicated against the authors' code.

Load packages and seed the RNG

using Turing
using Distributions
using StatsFuns: logistic
using DataFrames: DataFrame, eachrow
import CSV
using Random
using Markdown
using Dates: Date, Day, value
using BVDOutbreakSize
import CairoMakie
# Loading TensorBoardLogger activates the `tensorboard_callback` extension
# so `fit_callback` can stream each fit to TensorBoard as well as a progress
# log. Set `BVD_FIT_LOG=none` to disable all fit logging (CI release builds).
using TensorBoardLogger

# Render figures at higher resolution so they stay crisp in the docs.
CairoMakie.activate!(type = "png", px_per_unit = 3)

Random.seed!(20260518)

Random.TaskLocalRNG()

Methods

Data

The DRC data come from the situation reports of the Institut National de Santé Publique (Institut National de Santé Publique and Centre des Opérations d'Urgence de Santé Publique (COUSP-RDC), May 2026). Each report gives the national cumulative suspected cases and deaths, laboratory-confirmed cases and deaths, and the specimens received and analysed by the laboratory, at the report date. From SitRep 013 (27 May) INSP began reclassifying suspects, so the cumulative suspected count falls. We freeze it at its last stable vintage (26 May) and instead read the daily new-suspect count ("nouveaux cas suspects du jour") that the confirmed-based reports publish from 4 June, fitting it as a daily incidence where the cumulative series stops. The same reports print a daily new suspected-death count alongside it ("cas suspects du jour N (M deces)", from 7 June), fitted the same way where the cumulative suspected-death series stops. The confirmed-based reports also publish a daily "Patients en isolement" count, the number of patients (confirmed plus suspected) in an isolation/treatment bed at the end of the day; we fit it as the suspect inflow carried through a length-of-stay survival into a daily bed count. The fitted series runs from 1 June (SitRep 018), where the column is relabelled to the all-patients "Patients en isolement - hospitalisation"; the narrower suspects-only count in SitReps 016-017 is a different quantity and is left out. The reports also print a cumulative "cumul guéris" total of confirmed cases recorded as recovered, from 6 June; we fit it as survivors among the modelled confirmed cases (a scaled confirmation-to-recovery convolution, the incidence analogue of the isolation prevalence stream). We extracted these figures from the written situation-report PDFs (archived by INRB-UMIE (INRB-UMIE, 2026)) using a language model, with a second pass to re-read them, rather than the published per-zone CSVs. The zone sums in the CSVs are inconsistent with the national headline totals because they drop counts not yet attributed to a zone, so they understate the national totals. The Uganda data are the cases and the one death exported across the border, taken from the WHO situation reports and Disease Outbreak News (World Health Organization, May 2026). The cross-border traveller volume and source population come from McCabe et al. (McCabe and others, May 2026); the source population is fixed and the traveller volume is given a Normal prior around the McCabe et al. figure.

The first table lists each figure at the cut-off, or at the date reporting stopped for that stream. The second table gives the per-date history of each situation-report stream; the model fits the between-report increments of these series, so a single date reduces to the cut-off total.

Loading observations and building the data table

obs = load_observations()
# Grid day-index (day n is the cut-off) back to a calendar date.
grid_date(day) = obs.cutoff - Day(obs.n - day)
# Date a cumulative history last reports (the cut-off for streams that
# run to it, or the freeze date for streams that stop earlier).
hist_last_date(h) = isempty(h.days) ? missing : grid_date(maximum(h.days))
observations_table = DataFrame(
    field = [
        "exported_cases",
        "exports_deaths",
        "suspected_deaths",
        "suspected_cases",
        "confirmed_cases",
        "confirmed_deaths",
        "specimens_analysed",
        "genetic_tmrca_bound",
        "daily_outbound_travellers (prior mean)",
        "daily_outbound_travellers_sd (prior SD)",
        "source_population"
    ],
    date = [
        isempty(obs.export_case_days) ? missing :
        grid_date(maximum(obs.export_case_days)),
        isempty(obs.export_death_days) ? missing :
        grid_date(maximum(obs.export_death_days)),
        hist_last_date(obs.deaths_history),
        hist_last_date(obs.reported_history),
        hist_last_date(obs.confirmed_history),
        hist_last_date(obs.confirmed_deaths_history),
        hist_last_date(obs.lab_history),
        grid_date(obs.n - obs.tmrca_days),
        missing,
        missing,
        missing
    ],
    value = [
        obs.exported_cases,
        obs.exports_deaths,
        obs.total_deaths,
        obs.reported_cases,
        obs.confirmed_cases,
        obs.confirmed_deaths,
        obs.tests_analysed,
        obs.tmrca_days,
        ITURI_DAILY_TRAVEL,
        ITURI_DAILY_TRAVEL_SD,
        ITURI_POPULATION
    ]
);

11×3 DataFrame

Row	field	date	value
	String	Date?	Int64
1	exported_cases	2026-05-23	3
2	exports_deaths	2026-05-14	1
3	suspected_deaths	2026-05-26	246
4	suspected_cases	2026-05-26	1077
5	confirmed_cases	2026-06-21	1048
6	confirmed_deaths	2026-06-21	267
7	specimens_analysed	2026-05-28	755
8	genetic_tmrca_bound	2026-03-25	88
9	daily_outbound_travellers (prior mean)	missing	1871
10	daily_outbound_travellers_sd (prior SD)	missing	200
11	source_population	missing	4392200

The per-date cumulative history of the DRC situation-report streams, the national totals at each report date. The joint model fits the between-report increments of these series, so a single date reduces to the cut-off total. Two columns are the exception. suspected_new_daily is a per-day new-suspect count (not a cumulative total), fitted directly as a daily incidence, and it picks up where the cumulative suspected_cases column freezes on 26 May. patients_isolated is a daily count of patients in an isolation/treatment bed, fitted as the suspect inflow carried through a length-of-stay survival. See data/observations.toml for the per-stream sources.

Building the per-date time-series table

vintage_table = let
    # Each history carries grid day-indices and counts; key the counts
    # by calendar date so every stream lines up in one table.
    bydate(h) = Dict(grid_date(d) => c for (d, c) in zip(h.days, h.counts))
    streams = (
        suspected_cases = bydate(obs.reported_history),
        suspected_new_daily = bydate(obs.suspected_daily_history),
        patients_isolated = bydate(obs.isolation_history),
        suspected_deaths = bydate(obs.deaths_history),
        suspected_new_daily_deaths = bydate(obs.suspected_daily_deaths_history),
        confirmed_cases = bydate(obs.confirmed_history),
        confirmed_deaths = bydate(obs.confirmed_deaths_history),
        recovered_confirmed = bydate(obs.recovered_history),
        specimens_received = bydate(obs.tests_received_history),
        specimens_analysed = bydate(obs.lab_history)
    )
    dates = sort(collect(union((keys(s) for s in streams)...)))
    at(s) = [haskey(s, d) ? s[d] : missing for d in dates]
    DataFrame(
        date = dates,
        suspected_cases = at(streams.suspected_cases),
        suspected_new_daily = at(streams.suspected_new_daily),
        patients_isolated = at(streams.patients_isolated),
        suspected_deaths = at(streams.suspected_deaths),
        confirmed_cases = at(streams.confirmed_cases),
        confirmed_deaths = at(streams.confirmed_deaths),
        recovered_confirmed = at(streams.recovered_confirmed),
        specimens_received = at(streams.specimens_received),
        specimens_analysed = at(streams.specimens_analysed)
    )
end;

Per-date situation-report data table

36×10 DataFrame

Row	date	suspected_cases	suspected_new_daily	patients_isolated	suspected_deaths	confirmed_cases	confirmed_deaths	recovered_confirmed	specimens_received	specimens_analysed
	Date	Int64?	Int64?	Int64?	Int64?	Int64	Int64	Int64?	Int64?	Int64?
1	2026-05-14	missing	missing	missing	missing	8	4	missing	missing	missing
2	2026-05-17	missing	missing	missing	missing	13	4	missing	missing	missing
3	2026-05-18	516	missing	missing	131	33	4	missing	missing	missing
4	2026-05-19	575	missing	missing	148	51	4	missing	missing	missing
5	2026-05-20	672	missing	missing	160	64	6	missing	missing	missing
6	2026-05-21	745	missing	missing	175	83	9	missing	missing	missing
7	2026-05-22	872	missing	missing	204	91	10	missing	missing	missing
8	2026-05-23	904	missing	missing	220	101	10	missing	418	211
9	2026-05-24	906	missing	missing	223	105	10	missing	431	295
10	2026-05-25	998	missing	missing	238	106	12	missing	431	295
11	2026-05-26	1077	missing	missing	246	121	17	missing	662	403
12	2026-05-27	missing	missing	missing	missing	125	17	missing	774	648
13	2026-05-28	missing	missing	missing	missing	210	17	missing	883	755
14	2026-05-29	missing	missing	missing	missing	263	42	missing	missing	missing
15	2026-05-30	missing	missing	missing	missing	282	42	missing	missing	missing
16	2026-05-31	missing	missing	missing	missing	321	48	missing	missing	missing
17	2026-06-01	missing	missing	173	missing	344	60	missing	missing	missing
18	2026-06-02	missing	missing	206	missing	363	62	missing	missing	missing
19	2026-06-03	missing	missing	233	missing	381	64	missing	missing	missing
20	2026-06-04	missing	153	258	missing	452	82	missing	missing	missing
21	2026-06-05	missing	119	267	missing	488	86	missing	missing	missing
22	2026-06-06	missing	117	283	missing	515	91	12	missing	missing
23	2026-06-07	missing	94	309	missing	550	101	19	missing	missing
24	2026-06-08	missing	138	297	missing	598	115	22	missing	missing
25	2026-06-09	missing	119	260	missing	635	127	30	missing	missing
26	2026-06-10	missing	119	262	missing	676	136	32	missing	missing
27	2026-06-11	missing	168	315	missing	689	139	32	missing	missing
28	2026-06-13	missing	136	359	missing	782	181	40	missing	missing
29	2026-06-14	missing	165	363	missing	808	192	48	missing	missing
30	2026-06-15	missing	235	376	missing	837	196	49	missing	missing
31	2026-06-16	missing	192	379	missing	875	202	67	missing	missing
32	2026-06-17	missing	151	383	missing	896	232	78	missing	missing
33	2026-06-18	missing	238	416	missing	933	245	80	missing	missing
34	2026-06-19	missing	162	361	missing	956	247	92	missing	missing
35	2026-06-20	missing	201	365	missing	1003	254	100	missing	missing
36	2026-06-21	missing	202	371	missing	1048	267	112	missing	missing

Model

Model overview

We model a single outbreak seeded by a zoonotic introduction on a daily grid from a seeding date to the cut-off (day $n$). The generating infection process produces daily infection incidence via the discrete renewal equation

$$\begin{array}{}\text{(1)}& I_t=R_t\sum _{s=1}^L{I}_{t-s}{g}_s,\end{array}$$

where $g$ is the discretised probability mass function (PMF) of the generation interval, indexed from lag 1 so an infectee is always infected strictly after its infector, and $R_t$ is the per-day reproduction number. A PMF gives the probability assigned to each whole-day lag. We never observe infections directly. Each data stream observes a thinned, delayed or transformed view of the same latent incidence. This is the class of time-varying renewal model used in EpiNow2 (Abbott et al., 2020), with the streams fitted jointly here rather than in a pipeline.

The model is assembled from modular Turing (Ge et al., 2018) submodels, each holding the maths and priors for one part of the generative process. We describe them in generative order, from the infection process through the epidemiological delays to the observation streams. The implementation uses Mooncake (Tebbutt and Ge, 2024) reverse-mode automatic differentiation, CensoredDistributions for delay discretisation, FlexiChains for chain handling, and PairPlots (Thompson, 2024) with AlgebraOfGraphics (Danisch and Krumbiegel, 2021) for the figures. Each submodel's source is shown in the collapsible block beneath its prose.

The table below shows which parameters inform each observation submodel. The analysed column is the analysed-specimen volume, the single laboratory stream fitted as a count; the confirmed positives are scored as a Binomial of the observed analysed denominator with a positivity linked to the composition of the suspected pool, so the laboratory data help identify the non-BVD background. The conf. deaths column mirrors the laboratory pipeline on the death side, with a death testing fraction and a death-pool composition positivity built from the same assay:

Parameter	Exports	Deaths	Cases	Analysed	Confirmed	Conf. deaths	Export deaths
Reproduction number $R_t$	●	●	●	●	●	●	●
Generation interval	●	●	●	●	●	●	●
Incubation period	●	●	●	●	●	●	●
Seed $I_0$	●	●	●	●	●	●	●
Onset-to-death delay		●				●	●
Case-fatality ratio		●				●	●
Death ascertainment $p_{\text{death}}$		●				●
Background CFR ${\mathrm{cfr}}_{\text{bg}}$		●				●
Onset-to-report delay			●	●	●
Receipt delay				●	●	●
Onset-to-detection delay	●
Assay sensitivity / specificity					●	●
Severity enrichment ${\delta }_0$					●
Death testing fraction ${\tau }_{\text{death}}$						●
Testing fraction ${\tau }_{\text{test}}$				●	●
Background rate ${\lambda }_{\text{bg}}$		●	●	●	●	●
Surveillance dispersion		●	●	●
Ascertainment	●		●	●	●		●
Traveller volume	●						●

Infections

The infection process combines several components. These are the reproduction number, the generation interval that drives the renewal, the seeding that sets the initial infection count, the genetic bound on the outbreak age, the growth rate that fills the unobserved cryptic phase, and the renewal construction that grows the seed forward to the cut-off. Each is described in a subsection below.

Reproduction number

The reproduction number is held flat at the established reproduction number $R_0$ until a month before the first WHO situation report, then follows a non-centred Gaussian random walk on the log scale with weekly knots to the cut-off. The month-long lead lets $R_t$ start moving before the first report, since transmission may already have turned before the outbreak was formally reported; the walk start is floored at the renewal start. The walk starts from $R_0$ at its first knot:

$$\begin{array}{}\text{(2)}& \log R_k=\log R_0+{\sigma }_{\text{rw}}\sum _{j=1}^k{z}_j,z_j\sim \mathrm{Normal}(0,1),{\sigma }_{\text{rw}}\sim {\mathrm{Normal}}^{+}(0,0.1).\end{array}$$

We do not place a prior on $R_0$ directly. We put the prior on the initial growth rate $r$ instead, given in the seeding and growth subsection below, and derive the established reproduction number forward from it through the Euler–Lotka relation under our generation interval $g$:

$$\begin{array}{}\text{(3)}& R_0={\left(\sum _{s\ge 1}{g}_s{e}^{-rs}\right)}^{-1}.\end{array}$$

The step-size prior keeps weekly changes in the reproduction number moderate. We set the half-normal on ${\sigma }_{\text{rw}}$ so that the reproduction number is unlikely to change by more than about 20% from one week to the next: two standard deviations of the weekly log-step is around $0.20$.

Daily $\log R_t$ is the linear interpolation between the weekly knots, so the reproduction number varies piecewise linearly within each week; before the first knot it is held flat at $R_0$ (the interpolation clamps below the first knot rather than extrapolating):

$$\begin{array}{}\text{(4)}& \log R_t=\log R_k+\frac{t-d_k}{d_{k+1}-d_k}(\log R_{k+1}-\log R_k),d_k\le t\le d_{k+1},\end{array}$$

with $d_k$ the day of knot $k$. The outbreak response adds a sampled effect shaped by a logistic ramp at the first WHO situation report on 18 May 2026. We assume the response takes about three weeks (21 days) to take effect, and that it can only reduce transmission, so the effect is constrained to be non-positive:

$$\begin{array}{}\text{(5)}& \log R_t+=\delta \cdot \mathrm{logistic}\left(\frac{t-t_{\text{bp}}}{21}\right),\delta \sim {\mathrm{Normal}}^{-}(0,0.4).\end{array}$$

Submodel: rt_walk_model

@model function rt_walk_model(n::Integer, log_R0_base::Real;
        week::Integer = 7,
        breakpoint::Union{Missing, Real} = missing,
        rt_start::Integer = 1,
        ramp::Real = 21.0,
        sigma_prior = truncated(Normal(0, 0.1); lower = 0),
        effect_prior = truncated(Normal(0, 0.4); upper = 0))
    days = knot_days(n; week, start = rt_start)
    nb = length(days)
    ## The established `R0` at the genetic bound is the base the random walk
    ## grows from. It is DERIVED (forward Euler–Lotka from the sampled growth
    ## rate) and passed in, not sampled here; it is tracked as a deterministic
    ## so the walk base stays available on the chain. The days before the
    ## renewal start (`rt_start`) are filled by the analytic cryptic
    ## exponential in `infection_model`, so the walk values there are unused;
    ## the interpolation clamps to `log_R0` before the first knot, which is
    ## harmless.
    log_R0 := log_R0_base
    sigma_rw ~ sigma_prior
    z ~ product_distribution(fill(Normal(0, 1), max(nb - 1, 1)))
    intervention_effect ~ effect_prior
    steps = sigma_rw .* z[1:(nb - 1)]
    log_R = log_R0 .+ vcat(zero(log_R0), cumsum(steps))
    log_Rt = interpolate_knots(log_R, days, n)
    log_Rt = log_Rt .+ intervention_effect .* sigmoid_ramp(n, breakpoint; ramp)
    Rt = exp.(log_Rt)
    return (; Rt, log_R, days, sigma_rw, log_R0, intervention_effect)
end

Generation interval

We assume the generation interval $g$ is a Gamma distribution with a sampled shape $\alpha$ and scale $\theta$, taken from the Ebola virus disease serial interval used as a generation-time proxy (mean 15.3 d, SD 9.3 d; WHO Ebola Response Team 2014). That distribution maps once to a Gamma shape near $2.71$ and scale near $5.65$, and the priors are centred on those values, with spreads that carry the source's reported uncertainty on the mean rather than a spread we assign ourselves:

$$\begin{array}{}\text{(6)}& \alpha \sim {\mathrm{Normal}}^{+}(2.71,0.70),\theta \sim {\mathrm{Normal}}^{+}(5.65,1.50).\end{array}$$

The Gamma is discretised through the same double-interval-censoring route as every delay, described with the first epidemiological process model below, and the lag-0 bin is dropped and the remainder renormalised so the generation interval starts at one day and an infectee is infected strictly after its infector.

Submodel: generation_interval_model

@model function generation_interval_model(nmax::Integer;
        alpha_prior = truncated(Normal(2.71, 0.7); lower = 0.1),
        theta_prior = truncated(Normal(5.65, 1.5); lower = 0.1))
    α ~ alpha_prior
    θ ~ theta_prior
    dist = Gamma(α, θ)
    pmf = discretise_censored(dist, nmax)
    g = pmf[2:end] ./ sum(pmf[2:end])
    return (; g, gi_mean = α * θ, gi_sd = sqrt(α) * θ,
        gi_alpha = α, gi_theta = θ)
end

Seeding and growth

We assume the outbreak started from a single seed case introduced by a zoonotic spillover. The initial infection count $I_0$ on the last day of the seeding window has a prior centred on a single seed:

$$\begin{array}{}\text{(7)}& I_0\sim {\mathrm{Normal}}^{+}(0.1,0.1).\end{array}$$

From that seed we assume the outbreak grew deterministically through an unobserved cryptic exponential phase, doubling $m$ times before sustained transmission was established. The cryptic phase grows the seed to $2^m$ infections at the renewal start, the day the renewal takes over, over a duration $m\tau$ with $\tau$ the doubling time. The doubling count has a wide prior centred on three cryptic doublings:

$$\begin{array}{}\text{(8)}& m\sim {\mathrm{Normal}}^{+}(3,3),\tau =\frac{\log 2}{r},T_{\text{cryptic}}=m\tau .\end{array}$$

The growth rate $r$ carries the prior the genetic source informs. The genetic reanalysis reports the epidemic doubling time as 15.2 to 24.5 d across substitution-rate assumptions, with a centre near 20 d. We put a log-normal prior on $r$ equivalent to a log-normal prior on the doubling time centred on 20 d, with its log spread read from that range and inflated a little, so the prior is slightly wider in spread than the source but unbiased relative to it:

$$\begin{array}{}\text{(9)}& r\sim \mathrm{LogNormal}(\log \frac{\log 2}{20},0.15).\end{array}$$

This single growth rate fills the cryptic phase and, through the forward Euler–Lotka derivation above, sets the established reproduction number, so the cryptic phase and the renewal share one growth source. The genetic report's own established reproduction number of about $1.31$ to $1.55$ uses its own generation interval; deriving $R_0$ forward from the shared growth rate under our generation interval is the consistent choice.

Submodel: seed_model

@model function seed_model(; i0_prior = truncated(Normal(0.1, 0.1); lower = 0))
    I0 ~ i0_prior
    return (; I0)
end

Submodel: exponential_growth_model

@model function exponential_growth_model(;
        r_prior = LogNormal(log(log(2) / M_PRIOR_DOUBLING_DAYS), 0.15),
        m_prior = truncated(Normal(M_PRIOR_BASE, 3.0); lower = 0))
    r ~ r_prior
    m ~ m_prior
    τ := log(2) / r
    T := m * τ
    C_T := 2.0^m
    return (; τ, r, m, T, C_T)
end

Genetic bound on outbreak age

A BEAST time tree of the first ten sequenced genomes (Amuri-Aziza et al., 2026) places the TMRCA, the age of the oldest internal node of the tree, at a mean of 25 March 2026. The temporal sampling range is too short to estimate the molecular clock, so we fix it to the $1.2\times {10}^{-3}$ substitutions/site/year rate of the 2013-2016 West African Ebola epidemic (Holmes et al., 2016). The TMRCA is a lower bound on the outbreak age. Adding sequences, or more geographically representative ones, can only push it earlier, never later, as the sampled tree is almost entirely from Bunia. Using the genetic TMRCA as a one-sided seeding bound rather than a point estimate follows a suggestion of N. Ferguson (Ferguson, 2026).

We treat the TMRCA day as a right-censored, noisy reading of the total outbreak age $T$ (the cryptic duration plus the observed window, defined in the infection process below):

$$\begin{array}{}\text{(10)}& {\text{tmrca}}_{\text{days}}\sim \mathrm{censored}(\mathrm{Normal}(T,\sigma );\text{upper}={\text{tmrca}}_{\text{days}}),\sigma =15\text{d}.\end{array}$$

The renewal starts on the grid day on which the renewal recursion begins and sustained transmission is treated as established. We place it 14 days after the genetic TMRCA day, past the molecular-clock uncertainty, so the observed window from the renewal start to the cut-off is shorter than the TMRCA age. The bound therefore stays informative on the cryptic duration, pulling the origin to sit at or before the most recent common ancestor and bounding the cryptic phase from below. It is one-sided, leaving the age free above the TMRCA. We fix the clock and do not propagate cross-outbreak or clock uncertainty.

Submodel: genetic_seeding_model

@model function genetic_seeding_model(T::Real,
        tmrca_days::Union{Missing, Real}; tmrca_days_sd::Real = 15.0)
    if !ismissing(tmrca_days)
        tmrca_days ~ censored(Normal(T, tmrca_days_sd); upper = tmrca_days)
    end
    return (; T, tmrca_days_sd)
end

Infection process

The renewal start and observed window from the genetic bound above are

$$\begin{array}{}\text{(11)}& \text{renewal start}=n-{\text{tmrca}}_{\text{days}}+14,{\tau }_{\text{obs}}=n-\text{renewal start}.\end{array}$$

The grid days before the renewal start are filled by the cryptic exponential curve at rate $r$ ending at $2^m$, giving the recursion a full generation interval of history. The renewal then grows the trajectory forward under the time-varying reproduction number. The total outbreak age is the cryptic duration plus the observed window:

$$\begin{array}{}\text{(12)}& T=m\tau +{\tau }_{\text{obs}}.\end{array}$$

Cumulative infections are the running sum of the daily infection series. The cumulative infection count at the cut-off is the headline outbreak size. The current growth rate is the exponential growth implied by the cut-off reproduction number and the generation interval through forward Euler–Lotka, so it is sign-consistent with that number by construction, and the current doubling time is $\log 2$ divided by that rate.

Submodel: infection_model

@model function infection_model(n::Integer;
        breakpoint::Union{Missing, Real} = missing,
        rt_start::Integer = 1,
        rt_walk_start::Integer = rt_start,
        rt = rt_walk_model,
        gi = generation_interval_model,
        growth = exponential_growth_model,
        gi_nmax::Integer = cdf_nmax(Gamma(2.71, 5.65)))
    gi_state ~ to_submodel(gi(gi_nmax))
    g = gi_state.g
    ## ONE growth source: the prior is on the cryptic exponential growth rate
    ## `r` (sampled in `growth`), and the SINGLE established reproduction
    ## number `R0` (= the walk base, the first `R_t`) is derived FORWARD from
    ## that `r` and the generation interval through Euler–Lotka. The cryptic
    ## phase and the established renewal therefore share `r`.
    growth_state ~ to_submodel(growth())
    r_clock = growth_state.r
    R0 = r_to_R0(r_clock, g)
    ## The random walk's first knot sits at `rt_walk_start`, decoupled from
    ## the renewal start `rt_start`: the renewal seeds and grows from the
    ## genetic-TMRCA renewal start, but `R_t` is held flat at `R0` until
    ## `rt_walk_start` (the first situation report). Before any case or death
    ## surveillance the dynamics are unidentified, so a free walk there only
    ## adds unsupported drift. `rt_walk_start` defaults to `rt_start`, the
    ## walk-from-renewal-start case.
    rt_state ~ to_submodel(rt(n, log(R0); breakpoint, rt_start = rt_walk_start))
    Rt = rt_state.Rt
    ## renewal_start = genetic-TMRCA grid day (`rt_start`); the observation
    ## span is τ_obs = n − renewal_start. The renewal-start seed magnitude is
    ## `2^m` DIRECTLY (the cryptic phase grows one import to `2^m` over `m`
    ## doublings, `r`-free). Fill grid days 1…renewal_start with the cryptic
    ## exponential curve at rate `r` ending at `2^m` (a full GI of history),
    ## then run the renewal forward from renewal_start+1.
    renewal_start = clamp(rt_start, 1, n)
    τ_obs = n - renewal_start
    seed0 = seed_at_renewal_start(growth_state.C_T)
    seed_vec = seed_infections(seed0, r_clock, renewal_start)
    infections = renewal_infections(Rt, g, seed_vec)
    cumulative = cumsum(infections)
    ## Total outbreak age: cryptic duration (m·τ) plus the observation span.
    T_total = growth_state.T + τ_obs
    ## Current growth rate at the cut-off, derived from the cut-off
    ## reproduction number `Rt[n]` and the generation interval through forward
    ## Euler–Lotka (the inverse of the `r_to_R0` that derives `R0` from the
    ## clock growth above). This makes the reported current growth rate
    ## consistent with `R_T := Rt[n]` BY CONSTRUCTION: `r < 0` iff `R_T < 1`.
    ## An earlier formulation read `r` off the realised last-two-days slope
    ## `log I[n] − log I[n-1]`, but the intervention ramp depresses the final
    ## renewal step (`I[n] < I[n-1]` while `Rt[n] ≥ 1`), so that realised
    ## slope disagreed in sign with `R_T` at the cut-off — an end-of-
    ## trajectory edge artifact rather than the instantaneous growth.
    r = euler_lotka_r(@inbounds(Rt[n]), g)
    return (; infections, cumulative, Rt, g, seed_at_renewal_start = seed0,
        m = growth_state.m, τ = growth_state.τ, R0, r0 = r_clock, r,
        doubling_time_initial = doubling_time(r_clock),
        T = T_total, C_T = cumulative[n],
        C_T_prior = growth_state.C_T, doubling_time = doubling_time(r),
        seeding_age = seeding_age(cumulative, n))
end

Epidemiological process models

We model each observed stream as a delayed and thinned view of the daily onset incidence. This section gives the delays that map infections to onsets and onsets to each observed endpoint, and the case-fatality ratio that maps onsets to deaths. The incubation period comes first, then the onset-to-report delay (also used for export detection), the onset-to-death delay and the report-to-receipt delay, then the case-fatality ratio.

Incubation period

Infections are convolved with the incubation-period PMF to give daily symptom-onset incidence, computed once and consumed by every downstream observation stream. We use the Bundibugyo virus incubation-period estimate from the 2007 Uganda outbreak (mean 6.3 d, 95% CI 5.2-7.3, $n=24$; (MacNeil et al., 2010)). The mean prior reproduces that 95% CI; the source reports no interval on the spread, so the SD prior is our own choice:

$$\begin{array}{}\text{(13)}& {\mu }_{\text{inc}}\sim {\mathrm{Normal}}^{+}(6.3,0.54),{\sigma }_{\text{inc}}\sim {\mathrm{Normal}}^{+}(3.5,0.8).\end{array}$$

Every delay is discretised to a daily PMF over lags $0,\dots ,n_{max}$ by double interval censoring (Charniga et al., 2024). The delays the companion line-list reanalysis reports, namely the onset-to-admission delay (used for both suspected-case reporting and export detection) and the two onset-to-death components, are carried through on their natural Gamma shape and scale, with the reanalysis's reported uncertainty, like the generation interval above. The incubation period and the laboratory receipt delay are not in the line list, so they keep a mean-and-SD prior moment-matched to a LogNormal. The LogNormal and Gamma CDFs both differentiate cleanly under the reverse-mode automatic differentiation. The maximum lag $n_{max}$ is not hand-set: for each delay it is the 98th percentile of the prior-centre distribution, computed once outside the model, so every delay captures a consistent 98% of its mass before the truncated PMF is renormalised.

Both the primary event (the onset, say) and the secondary event (the report) are observed only to the day, so the discretisation censors both. The primary event is taken uniform over its day and the secondary event is interval-censored to its day, giving the daily PMF

$$\begin{array}{}\text{(14)}& f_s={\int }_0^1[F(s+1-u)-F(s-u)]\mathrm{d}u,F=\text{the delay CDF},\end{array}$$

which is then renormalised over lags $0,\dots ,n_{max}$.

The incubation period also enters the infection-to-detection and infection-to-death delays for the export streams, where the survival clock runs from infection rather than onset.

Submodel: onset_incidence_model

@model function onset_incidence_model(infections::AbstractVector;
        incubation = (nmax) -> censored_delay_model(nmax;
            mean_prior = truncated(Normal(6.3, 0.54); lower = 1),
            sd_prior = truncated(Normal(3.5, 0.8); lower = 1)),
        incubation_nmax::Integer = cdf_nmax(lognormal_meansd(6.3, 3.5)))
    inc_state ~ to_submodel(incubation(incubation_nmax))
    onsets = convolve_delay(infections, inc_state.pmf)
    return (; onsets, incubation_pmf = inc_state.pmf,
        incubation_mean = inc_state.mean, incubation_sd = inc_state.sd)
end

Submodel: censored_delay_model

@model function censored_delay_model(nmax::Integer; mean_prior, sd_prior)
    delay_mean ~ mean_prior
    delay_sd ~ sd_prior
    dist = lognormal_meansd(delay_mean, delay_sd)
    return (; pmf = discretise_censored(dist, nmax), dist,
        mean = delay_mean, sd = delay_sd)
end

Onset-to-report delay

The delay from symptom onset to a suspected case being detected and reported into surveillance. We use a Bayesian reanalysis (Funk and Abbott, 2026) of the 2012 Isiro Bundibugyo virus outbreak line list (Rosello and others, 2015), taking its onset-to-admission delay as a Gamma sampled on its natural shape and scale, with priors centred on the reanalysis posterior (implied mean about 4 d) and carrying its reported uncertainty:

$$\begin{array}{}\text{(15)}& {\alpha }_{\text{rep}}\sim {\mathrm{Normal}}^{+}(1.18,0.28),{\theta }_{\text{rep}}\sim {\mathrm{Normal}}^{+}(3.69,1.20).\end{array}$$

We do not use the reanalysis onset-to-notification delay, a near-exponential Gamma with mean about 20 d. We assume that delay reflects a longer notification pathway, likely including laboratory confirmation and administrative processing, rather than the rapid surveillance report we model. This delay drives the suspected-case, laboratory and confirmed-death streams, and the export model uses the same onset-to-admission delay for detection abroad.

Onset-to-death delay

McCabe et al. take the onset-to-death delay from the same line list as a point estimate (Rosello and others, 2015), fitting a $t$-distributed delay. The reanalysis instead fits it as two atomic Gamma components, onset-to-admission and admission-to-death, and convolves them. We do the same: each component is a Gamma sampled on its natural shape and scale, with priors centred on the reanalysis posteriors:

$$\begin{gathered} \begin{array}{}\text{(16)}& {\alpha }_{\text{oa}}\sim {\mathrm{Normal}}^{+}(1.18,0.28),{\theta }_{\text{oa}}\sim {\mathrm{Normal}}^{+}(3.69,1.20), \\ {\alpha }_{\text{ad}}\sim {\mathrm{Normal}}^{+}(2.15,0.60),{\theta }_{\text{ad}}\sim {\mathrm{Normal}}^{+}(3.91,1.38).\end{array} \end{gathered}$$

and the onset-to-death PMF is the convolution of the two discretised components (implied mean about 13 d). The source is shown with the deaths submodel below, where the delay is injected.

Onset-to-detection delay (exports)

An exported case is detected at a point of entry abroad when it first enters surveillance, the same event as a domestic suspected-case report, so the export model uses the same line-list onset-to-admission delay (Funk and Abbott, 2026) as the onset-to-report delay above, with the same natural shape and scale priors:

$$\begin{array}{}\text{(17)}& {\alpha }_{\text{det}}\sim {\mathrm{Normal}}^{+}(1.18,0.28),{\theta }_{\text{det}}\sim {\mathrm{Normal}}^{+}(3.69,1.20).\end{array}$$

It drives the exports streams; its source is shown with the exports submodel below.

Report-to-analysed delay

The delay from a suspected case being reported to its specimen being analysed by the laboratory, centred on a short turnaround with a heavy right tail allowing for specimen shipment to a confirmatory laboratory and the analysis queue. No per-sample outbreak data grounds this, so the prior is our own choice:

$$\begin{array}{}\text{(18)}& {\mu }_{\text{rec}}\sim {\mathrm{Normal}}^{+}(4.5,1.0),{\sigma }_{\text{rec}}\sim {\mathrm{Normal}}^{+}(4.0,0.75).\end{array}$$

It drives the laboratory analysed-specimen volume; its source is shown with the laboratory submodel below.

Case-fatality ratio

The US Centers for Disease Control and Prevention (CDC) summary for the two previous BVD outbreaks is $55$ deaths in $169$ cases ($\approx 33\mathrm{\%}$; CDC outbreak history), with confidence bands spanning roughly $26$-$40\mathrm{\%}$. The companion Bundibugyo virus (BDBV) reanalysis reports a baseline of $0.47$ ($95\mathrm{\%}$ CrI $0.31$-$0.65$) for non-healthcare-worker (non-HCW) confirmed cases. Based on this we use a prior of

$$\begin{array}{}\text{(19)}& \mathrm{CFR}\sim \mathrm{Beta}(6.6,13.4),\end{array}$$

with mean $0.33$ and $95\mathrm{\%}$ interval roughly $0.15$-$0.54$. The mean matches the CDC $55/169\approx 33\mathrm{\%}$ figure and the corrected central CFR in the 20 May report (McCabe and others, May 2026).

Submodel: cfr_model

@model function cfr_model(; cfr_prior = Beta(6.6, 13.4))
    CFR ~ cfr_prior
    return (; CFR)
end

The prior density, with the CDC $0.33$ figure marked.

Observation models

Each observation submodel takes the shared daily onset incidence, convolves it with a sampled onset-to-event delay, scales by the relevant ascertainment, case-fatality ratio or positivity factor, and reads the modelled count off the daily series at each vintage day. Likelihoods score the between-vintage increments. The surveillance streams come first, then the geographic-spread exports.

Shared observation submodels

Several parameters are assumed shared across the streams: the surveillance dispersion, the ascertainment fractions, the laboratory testing priors and the traveller volume. We assume the passive-surveillance count datasets are overdispersed and share a common dispersion, and that the laboratory testing priors are shared between the suspected-case and laboratory streams. More detail is given in the subsections below.

Surveillance dispersion

Each passive-surveillance count stream has its own negative-binomial dispersion, partially pooled across the streams so the sparse ones borrow strength. Following Stan prior-choice recommendations (Stan Development Team, 2024), the dispersion is sampled on the $1/\sqrt{k}$ scale in non-centred log form:

$$\begin{array}{}\text{(20)}& \log (1/\sqrt{k_s})=\mu +\tau z_s,z_s\sim \mathrm{Normal}(0,1),\mu \sim \mathrm{Normal}(\log 0.6,0.33),\tau \sim {\mathrm{Normal}}^{+}(0,0.3),\end{array}$$

so $k_s=1/\exp (\mu +\tau z_s{)}^2$ per stream, with $\tau$ setting the pooling ($\tau =0$ collapses to one shared dispersion). The population value $k=1/\exp (\mu {)}^2$ is the headline dispersion.

Submodel: pooled_dispersion_model

@model function pooled_dispersion_model(n_streams::Integer;
        mean_prior = Normal(log(0.6), 0.33),
        sd_prior = truncated(Normal(0, 0.3); lower = 0))
    μ_log ~ mean_prior
    τ ~ sd_prior
    z ~ product_distribution(fill(Normal(0, 1), max(n_streams, 1)))
    inv_sqrt_k = exp.(μ_log .+ τ .* z[1:n_streams])
    k = 1.0 ./ (inv_sqrt_k .^ 2 .+ eps(eltype(inv_sqrt_k)))
    k_pop = 1.0 / (exp(μ_log)^2 + eps(typeof(float(μ_log))))
    return (; k, inv_sqrt_k, k_pop, μ_log, τ)
end

Ascertainment

Two surveillance systems detect cases: DRC passive community surveillance (the reported suspected-case count) and Uganda's point-of-entry / hospital surveillance (the exported-case count). Each captures a fraction of the true cases passing through it. The two ascertainment fractions $p_{\text{DRC}}$ and $p_{\text{Uganda}}$ share a logit-scale hyperprior with mean $\mu$ and pooling strength $\tau$, centred on a reporting fraction of $75\mathrm{\%}$, reflecting the active case-finding of a declared Ebola response rather than baseline passive surveillance:

$$\begin{array}{}\text{(21)}& \mu \sim \mathrm{Normal}(\mathrm{logit}(0.75),1),\tau \sim {\mathrm{Normal}}^{+}(0,0.5),\end{array}$$$$\begin{array}{}\text{(22)}& \mathrm{logit}(p_{\text{DRC}})\sim \mathrm{Normal}(\mu ,\tau ),\mathrm{logit}(p_{\text{Uganda}})\sim \mathrm{Normal}(\mu ,\tau ).\end{array}$$

The cases likelihood uses $p_{\text{DRC}}$; the two Uganda-side likelihoods use $p_{\text{Uganda}}$.

Submodel: pooled_ascertainment_model

@model function pooled_ascertainment_model(;
        mu_prior = Normal(logit(0.75), 1.0),
        tau_prior = truncated(Normal(0, 0.5); lower = 1e-4))
    μ_logit ~ mu_prior
    τ_logit ~ tau_prior
    z_drc ~ Normal(0, 1)
    z_uganda ~ Normal(0, 1)
    logit_p_drc = μ_logit + τ_logit * z_drc
    logit_p_uganda = μ_logit + τ_logit * z_uganda
    p_drc := logistic(logit_p_drc)
    p_uganda := logistic(logit_p_uganda)
    return (; μ_logit, τ_logit, p_drc, p_uganda)
end

Laboratory priors

We model the process of confirming cases via laboratory testing. The testing fraction ${\tau }_{\text{test}}$ is the share of suspected cases routed to the laboratory. A truly BVD specimen tests positive with the assay sensitivity $s$, and a non-BVD specimen tests positive with the false-positive rate $1-\mathrm{spec}$ from the assay specificity. We assume that more severe cases, more likely to be Ebola, are preferentially tested, through an enrichment factor ${\delta }_0$ that raises the tested BVD share above the suspect-pool composition early on and relaxes towards it as testing broadens. The confirmed deaths mirror this laboratory pipeline rather than enriching the case composition: a fraction ${\tau }_{\text{death}}$ of suspected deaths reach the laboratory, and they confirm at the assay positivity $p=s{q}_{\text{death}}+(1-\mathrm{spec})(1-q_{\text{death}})$ built from the same sensitivity and specificity but the death-pool BVD share $q_{\text{death}}$. Confirmation runs on the altona RealStar Filovirus Screen RT-PCR rather than the Zaire-specific GeneXpert Ebola assay, which does not reliably detect Bundibugyo virus. Sensitivity for Bundibugyo virus is less well characterised than for Zaire ebolavirus, so we centre the sensitivity prior below the values reported for other strains and give it a wide spread. The specificity is high but imperfect; the severity enrichment is moderate and one-sided (triage upsamples BVD, never down); the death testing-intensity scaling is a tight log-normal centred on one, since no death-testing data grounds it:

$${\tau }_{\text{test}}\sim \mathrm{Beta}(5,2),s\sim \mathrm{Beta}(10,1.76),\mathrm{spec}\sim \mathrm{Beta}(60,2),$$$$\begin{array}{}\text{(23)}& {\delta }_0\sim {\mathrm{Normal}}^{+}(1.5,0.75),\text{scaling}\sim \mathrm{LogNormal}(0,0.25).\end{array}$$

The non-BVD background rate ${\lambda }_{\text{bg}}$ enters the suspected-case stream and is described with it below; the suspected deaths carry a death ascertainment $p_{\text{death}}\sim {\mathrm{logit}}^{-1}\mathrm{Normal}(\mathrm{logit}0.9,0.5)$ and a non-BVD death background tied to the case background by a background CFR ${\mathrm{cfr}}_{\text{bg}}\sim \mathrm{Beta}(2,6)$.

Submodel: test_positivity_model

@model function test_positivity_model(;
        lambda_prior = truncated(Normal(0.0, 1.0); lower = 0),
        fraction_tested_prior = Beta(5.0, 2.0))
    λ_bg ~ lambda_prior
    τ_test ~ fraction_tested_prior
    return (; λ_bg, τ_test)
end

Submodel: confirmed_positivity_model

@model function confirmed_positivity_model(nv::Integer;
        baseline_prior = Normal(logit(0.28), 0.7),
        pooling_prior = truncated(Normal(0.0, 1.0); lower = 0))
    m = max(nv, 1)
    q_mu ~ baseline_prior
    σ_q ~ pooling_prior
    z_q ~ product_distribution(fill(Normal(0, 1), m))
    logit_p = q_mu .+ σ_q .* z_q[1:nv]
    p_pos := logistic.(logit_p)
    return (; p_pos, q_mu, σ_q)
end

Submodel: test_sensitivity_model

@model function test_sensitivity_model(;
        sensitivity_prior = Beta(10.0, 1.76))
    s_test ~ sensitivity_prior
    return (; s_test)
end

Submodel: test_specificity_model

@model function test_specificity_model(; specificity_prior = Beta(60.0, 2.0))
    spec ~ specificity_prior
    return (; spec)
end

Submodel: severity_enrichment_model

@model function severity_enrichment_model(;
        logodds_prior = truncated(Normal(1.5, 0.75); lower = 0),
        decay_prior = truncated(Normal(0.0, 200.0); lower = 0.0))
    δ0 ~ logodds_prior
    decay_scale ~ decay_prior
    return (; δ0, decay_scale)
end

Submodel: death_testing_fraction_model

@model function death_testing_fraction_model(; fraction_prior = Beta(5.0, 2.0))
    τ_death ~ fraction_prior
    return (; τ_death)
end

Submodel: death_ascertainment_model

@model function death_ascertainment_model(;
        ascertainment_prior = Normal(logit(0.9), 0.5))
    logit_p_death ~ ascertainment_prior
    p_death := logistic(logit_p_death)
    return (; p_death, logit_p_death)
end

Submodel: background_cfr_model

@model function background_cfr_model(; cfr_prior = Beta(2.0, 6.0))
    cfr_bg ~ cfr_prior
    return (; cfr_bg)
end

Traveller volume

The number of people crossing from the source area to Uganda each day sets the travel rate in the exports likelihood. We treat it as an estimated quantity rather than a fixed input. McCabe et al. Table 3 records mean weekly passenger counts across seven points of entry; the Ituri-side daily total of $1871$ is a sample mean across roughly $15$-$21$ point-of-entry-weeks. We use a Normal prior centred on $1871$ with SD $200$ ($\approx 10\mathrm{\%}$ CV), truncated at zero, covering point-of-entry variation and the sitrep sampling uncertainty; the source population is kept fixed (census):

$$\begin{array}{}\text{(24)}& N_{\text{travel}}\sim {\mathrm{Normal}}^{+}(1871,200).\end{array}$$

Submodel: traveller_volume_model

@model function traveller_volume_model(;
        mean::Real = ITURI_DAILY_TRAVEL,
        sd::Real = ITURI_DAILY_TRAVEL_SD)
    daily_travellers ~ truncated(Normal(mean, sd); lower = 0)
    return (; daily_travellers)
end

Reported cases

Reported suspected cases are the sum of two parts. The first is a BVD-driven component: the daily onsets convolved with the onset-to-report delay $f_{\text{rep}}$ and scaled by the DRC ascertainment $p_{\text{DRC}}$. The convolution of a daily series $x$ with a delay PMF $f$ is the lagged sum

$$(x\ast f{)}_t=\sum _{s\ge 0}{x}_{t-s}{f}_s,$$

used for every delay below. We write the BVD onset-to-report series at unit ascertainment as

$${\text{bvd}}_t=\sum _{s\ge 0}{\text{onsets}}_{t-s}{f}_{\text{rep},s}.$$

The second part is an additive non-BVD background, so a suspected case need not be a true BVD infection. It is a per-day rate ${\lambda }_{\text{bg},t}$ that follows a lognormal random walk on weekly knots around a baseline ${\lambda }_{\mu }$, linearly interpolated to the daily grid,

$${\lambda }_{\text{bg},t}={\lambda }_{\mu }\exp (w_t),w_t=\mathrm{interp}({\sigma }_{\text{rw}}\sum _{s<k}{z}_s),z_s\sim \mathcal{N}(0,1),$$

gated to zero before the surveillance onset (a report-to-receipt lead before the first suspected-case report — the background does not exist before surveillance began) and shared, with one tight innovation SD ${\sigma }_{\text{rw}}$, between the suspected-case and suspected-death streams. Weekly knots match the reproduction-number walk and keep the background a gentle drift over a small number of innovations. The baseline carries a half-normal ${\mathrm{Normal}}^{+}(0,8)$ prior on the natural scale. A log-scale level would have a heavy right tail the background/outbreak-size degeneracy could exploit, whereas the natural-scale half-normal bounds it. It is wide enough that the laboratory positivity (only $210/755\approx 0.28$ of analysed specimens are positive) identifies the background, which is inferred to be the majority of the suspect pool. The daily expected suspected case count is

$$c_t=p_{\text{DRC}}{\text{bvd}}_t+{\lambda }_{\text{bg},t}.$$

The per-vintage increments are scored with a NegBinomial sharing the dispersion $k$:

$$\begin{array}{}\text{(25)}& Y_{\text{cases},i}-Y_{\text{cases},i-1}\sim \mathrm{NegBinomial}(\sum _{t=d_{i-1}+1}^{d_i}{c}_t,k).\end{array}$$

From SitRep 013 (27 May) INSP reclassifies suspects, so the national cumulative suspected total falls. We freeze it at 26 May and instead fit the daily new-suspect count that the confirmed-based reports publish (the "nouveaux cas suspects du jour" $a_j$ on report day $t_j$, 4-7 June). This is a genuine daily incidence, not a cumulative total, so it is scored against the modelled daily suspected count $c_{t_j}$ on that day directly (a single-day mean, not a between-vintage sum) with a NegBinomial sharing $k$:

$$a_j\sim \mathrm{NegBinomial}(c_{t_j},k).$$

The daily report days fall strictly after the frozen cumulative series ends, so the two suspected likelihoods cover disjoint days and do not double-count. The suspected-death stream is fitted the same way: the cumulative suspected-death total freezes at 26 May and the daily new suspected-death count ("cas suspects du jour N (M deces)", from 7 June) is scored against the modelled daily suspected-death count on each report day with a NegBinomial sharing $k$ (see deaths_model).

Submodel: reported_cases_model

@model function reported_cases_model(
        reported_history,
        reported_cases::Union{Missing, Integer},
        onsets::AbstractVector, k::Real, p_drc::Real;
        suspected_daily_history = (; days = Int[], counts = Int[]),
        positivity = test_positivity_model(),
        background_re = nothing,
        ## Onset to a suspected case being detected/reported, from the
        ## line-list onset→admission delay (d_oa, ~4 d): a case enters
        ## surveillance when first formally seen, so one delay serves both the
        ## suspect-case and export streams. The line-list onset→notification
        ## delay (~20 d) is NOT used — we assume it reflects a longer pathway
        ## (likely confirmation and administrative processing), though what it
        ## captures is uncertain.
        onset_to_report = gamma_delay_model(cdf_nmax(Gamma(1.178, 3.694));
            alpha_prior = truncated(Normal(1.178, 0.285); lower = 0.01),
            theta_prior = truncated(Normal(3.694, 1.198); lower = 0.1)))
    pos_state ~ to_submodel(positivity)
    report_state ~ to_submodel(onset_to_report)
    λ_bg = pos_state.λ_bg
    τ_test = pos_state.τ_test
    report_pmf = report_state.pmf

    ## Unit-ascertainment BVD onset-to-report daily series, reused by the
    ## confirmed stream.
    bvd_reports_daily = convolve_delay(onsets, report_pmf)

    n = length(bvd_reports_daily)
    vobs = vintage_obs(reported_history, reported_cases, n)

    ## Daily non-BVD background. With `background_re === nothing` this is
    ## the constant scalar `λ_bg` over the grid (the renewal default). When
    ## `background_re` is injected it is the smooth daily random-walk
    ## background ([`background_walk_model`](@ref)): a length-`n` daily series
    ## that is zero before the surveillance onset and follows a tight lognormal
    ## random walk after it, so the baseline `λ_bg` from `positivity` is
    ## overridden by the walk's level `λ_mu` and the background varies smoothly
    ## day-to-day rather than in per-vintage steps.
    if background_re === nothing
        λ_bg_base = λ_bg
        bg_sigma = zero(λ_bg)
        bg_daily = fill(λ_bg, n)
    else
        bg_state ~ to_submodel(background_re(n))
        λ_bg_base = bg_state.λ_mu
        bg_sigma = bg_state.σ_bg
        bg_daily = bg_state.λ
    end

    ## Suspected daily cases add the p_drc-scaled BVD signal and the
    ## non-BVD background.
    reports_daily = p_drc .* bvd_reports_daily .+ bg_daily

    modelled_increments = bin_increments(reports_daily, vobs.days)
    reported_increments ~ to_submodel(
        vintage_increments_model(modelled_increments, vobs.obs_increments, k))

    ## Daily new-suspect inflow ("nouveaux cas suspects du jour"): per-day
    ## counts scored against the modelled daily series at each report day. The
    ## mean for day `d` is the single-day `reports_daily[d]` (clamped into the
    ## grid), NOT a between-vintage increment — this is a genuine daily
    ## incidence, so it never differences a falling cumulative. Empty by
    ## default; a `missing` count vector samples (the predictive path).
    sd_days = suspected_daily_history.days
    sd_modelled = [reports_daily[clamp(Int(d), 1, n)] for d in sd_days]
    sd_obs = isempty(suspected_daily_history.counts) ? missing :
             collect(Int.(suspected_daily_history.counts))
    suspected_daily ~ to_submodel(
        vintage_increments_model(sd_modelled, sd_obs, k))

    raw_total = sum(reports_daily)
    expected_reports := safe_rate(raw_total)

    ## Implied per-suspected positivity at the cut-off: BVD share of the
    ## expected suspected total.
    bvd_total = p_drc * sum(bvd_reports_daily)
    positivity := safe_rate(bvd_total) / expected_reports

    ## Cumulative background suspected cases over the grid, exposed for
    ## comparison with the observed suspected total.
    bg_total = sum(bg_daily)

    return (; p_drc, λ_bg = λ_bg_base, τ_test, report_pmf, bvd_reports_daily,
        reports_daily, expected_reports, positivity, bg_daily, bg_sigma,
        bg_total)
end

Isolation occupancy

The "Patients en isolement" figure is the daily count of occupied isolation/treatment beds. Bed occupancy may be supply-driven, with demand for beds able to outstrip supply and occupancy catching up as capacity expands. To allow for that we model a latent bed demand and the supply-limited occupancy it produces rather than occupancy directly.

The latent demand is the suspect inflow carried through a length-of-stay survival $S(\tau )=P(\text{LOS}\ge \tau )$ (the renewal analogue of the convolution secondary-observation model of EpiNow2 (Abbott et al., 2020)). A proportion $p_{\text{iso}}$ of the reported suspects need a bed. The suspects are a BVD/background mixture leaving on different clocks, so the demand is the sum of two survival convolutions. The BVD demand uses the treatment length-of-stay $S_{\text{BVD}}$, the time an admitted BVD case occupies a bed, with the admission-to-death delay from the line-list reanalysis (Funk and Abbott, 2026) as its prior. The non-BVD demand uses the rule-out stay $S_{\text{ruleout}}$, how long a ruled-out suspect occupies a bed before discharge, with the report-to-receipt laboratory turnaround as its prior,

$$D_t=p_{\text{iso}}[\sum _{s\ge 0}{p}_{\text{DRC}}{\text{bvd}}_{t-s}{S}_{\text{BVD}}(s)+\sum _{s\ge 0}{\lambda }_{\text{bg},t-s}{S}_{\text{ruleout}}(s)].$$

The bed capacity $C(t)$ is a non-decreasing random walk on weekly knots, since beds are added over the response and not taken away, pinned by the implied bed count, the reported occupancy (the "Patients en isolement" count) divided by the reported "Taux d'occupation" rate ($\approx 400\to 452$ beds over 9–13 June).

The occupied-bed count is the latent demand right-censored at the bed capacity: while demand is below capacity the count tracks it, and once demand reaches capacity the count is censored there. The censoring bound is fixed at the recorded implied capacity $C_j^{\text{cap}}$ so the latent demand is left uncensored,

$$O_j\sim \mathrm{censored}(\mathrm{NegBinomial}(D_{t_j},k_{\text{iso}});\text{upper}=C_j^{\text{cap}}),C_j^{\text{obs}}\sim \mathrm{NegBinomial}(C_{t_j},k_{\text{iso}}),$$

with a dispersion $k_{\text{iso}}$ of its own (not shared with the other streams). Occupancy below capacity identifies the demand directly; the part of demand above a saturated capacity is only partially identified, since occupancy says demand was at least the beds filled and not how much more, so the bed shortfall above capacity is informed by the demand model and its priors rather than measured by the occupancy. The model exposes the cut-off occupancy, the cut-off bed demand (the need under unconstrained supply), their difference (the bed shortfall) and the utilisation $O_T/C$.

One limitation is that this is a single national model, with one national bed capacity and one national demand, so it cannot represent local saturation. On 13 June Ituri was at 93.9% occupancy while Sud-Kivu was at 21.9%, so beds free in one province cannot serve patients who need them in another, and the national capacity averages over a saturated epicentre and slack elsewhere. The national shortfall therefore understates the local unmet need. The renewal model does not carry per-province inflow, so it cannot be split into the per-province bed model at which the supply constraint actually operates. A second limitation is that capacity is taken as a single (slowly varying) national quantity even though beds are being added.

The exposed BVD share is the true-BVD fraction of demand (BVD-confirmed plus BVD-suspect), not the report's confirmed/suspect split. The fitted occupancy series is the all-patients column from 1 June (SitRep 018) onward.

Submodel: treatment_admission_model

@model function treatment_admission_model(
        isolation_history,
        bvd_reports_daily::AbstractVector,
        bg_daily::AbstractVector,
        p_drc::Real;
        capacity_history = (; days = Int[], counts = Int[]),
        admission = isolation_admission_model(),
        severity = isolation_severity_model(),
        capacity = bed_capacity_walk_model,
        dispersion = surveillance_dispersion_model(),
        ## Dispersion can be injected from the joint composer's pooled set
        ## (`k_external`); standalone it samples its own from `dispersion`.
        k_external::Union{Nothing, Real} = nothing,
        ## BVD treatment stay (admission → outcome): the in-hospital length of
        ## stay an admitted BVD case occupies a bed before leaving by death or
        ## discharge. Carried through on the natural Gamma shape/scale from our
        ## BDBV line-list reanalysis (the `bdbv-linelist-analysis` submodule)
        ## with its posterior uncertainty — the admission→death atomic delay,
        ## the same Gamma the onset→death convolution uses (mean ≈ 8.4 d). The
        ## reanalysis admission→discharge stay is close (≈7.7 d), so this single
        ## distribution stands in for the admission→outcome stay; the
        ## death/discharge mixture is a separate refinement. The truncation
        ## covers the 99th percentile so the long-stay survival tail is not
        ## clipped.
        bvd_los = gamma_delay_model(
            cdf_nmax(Gamma(2.151, 3.906); q = 0.99);
            alpha_prior = truncated(Normal(2.151, 0.604); lower = 0.01),
            theta_prior = truncated(Normal(3.906, 1.381); lower = 0.1)),
        ## Sampled non-BVD rule-out stay: how long a ruled-out suspect occupies
        ## an isolation bed before discharge. Centred on the report-to-receipt
        ## laboratory turnaround but a separate parameter from the
        ## lab-turnaround `receipt_pmf`, so the occupancy identifies it on its
        ## own clock.
        ruleout_los = censored_delay_model(
            cdf_nmax(lognormal_meansd(4.5, 4.0); q = 0.99);
            mean_prior = truncated(Normal(4.5, 2.0); lower = 1),
            sd_prior = truncated(Normal(4.0, 1.5); lower = 1)))
    adm_state ~ to_submodel(admission)
    p_iso = adm_state.p_iso
    sev_state ~ to_submodel(severity)
    ## BVD suspects are admitted at a higher rate than non-BVD rule-outs,
    ## skewed up from the base `p_iso` by the severity log-odds `δ_iso`:
    ## triage admits the sicker patients and BVD presents more severely. The
    ## skew enriches BVD among the admitted without conditioning on the
    ## unobserved BVD status of any individual suspect.
    p_iso_bvd = logistic(logit(p_iso) + sev_state.δ_iso)
    if k_external === nothing
        disp_state ~ to_submodel(dispersion)
        k = disp_state.k
    else
        k = k_external
    end
    n = length(bvd_reports_daily)
    ## Time-varying bed capacity `C(t)` (a random walk), so the ceiling tracks
    ## the beds being added rather than being fixed. The walk starts at the
    ## first day with occupancy or capacity data, so off-window capacity stays
    ## flat at the baseline rather than carrying unidentified innovations.
    cap_obs_days = vcat(Int.(isolation_history.days),
        Int.(capacity_history.days))
    cap_start = isempty(cap_obs_days) ? 1 : minimum(cap_obs_days)
    cap_state ~ to_submodel(capacity(n; start = cap_start))
    C = cap_state.C
    C_T = isempty(C) ? zero(eltype(C)) : C[end]
    bvd_los_state ~ to_submodel(bvd_los)
    ruleout_los_state ~ to_submodel(ruleout_los)

    ## Latent bed demand: BVD suspects are admitted at `p_iso_bvd`, non-BVD
    ## suspects at the base `p_iso`. BVD patients stay for the sampled
    ## treatment length-of-stay; non-BVD suspects leave once they are ruled
    ## out and discharged, after the separately sampled rule-out stay.
    bvd_demand = convolve_survival(p_iso_bvd .* p_drc .* bvd_reports_daily,
        bvd_los_state.pmf)
    bg_demand = convolve_survival(p_iso .* bg_daily, ruleout_los_state.pmf)
    demand = bvd_demand .+ bg_demand
    ## In predict / check-model mode the daily series can widen to
    ## `Vector{Any}`, which the `min.(demand, C)` below cannot broadcast over,
    ## so pin it to the capacity's (always-concrete) element type, leaving the
    ## AD/fit path untouched.
    if eltype(demand) === Any
        demand = convert(Vector{eltype(C)}, demand)
    end

    ## Supply-limited occupancy: the latent demand right-censored at the bed
    ## capacity. The censoring bound is the recorded implied-capacity series
    ## ([`censoring_cap`](@ref)), a fixed data bound rather than a sampled
    ## ceiling, so the censored NegativeBinomial has no moving `-Inf` wall to
    ## diverge against. The latent demand is left UNCENSORED, so demand above
    ## capacity is carried by the renewal / length-of-stay demand model and its
    ## priors; the part of demand above a saturated capacity is only partially
    ## identified from occupancy. `occupancy = min(demand, C)` is the
    ## supply-capped stock the derived quantities report; the capacity walk
    ## `C(t)` carries the implied-capacity likelihood and the forecast cap.
    occupancy = min.(demand, C)

    ## Each report day's occupied-bed count is a NegativeBinomial around the
    ## latent demand, right-censored at the fixed capacity bound. Empty history
    ## is a no-op; a `missing` count vector samples (the predictive path).
    iso_days = isolation_history.days
    iso_obs = isempty(isolation_history.counts) ? missing :
              collect(Int.(isolation_history.counts))
    iso_means = [demand[clamp(Int(d), 1, n)] for d in iso_days]
    iso_ceil = censoring_cap(iso_days, iso_obs, capacity_history)
    isolation ~ to_submodel(
        censored_occupancy_model(iso_means, iso_ceil, iso_obs, k))

    ## Bed capacity: the implied bed count (occupancy / reported occupancy
    ## rate) on the days a rate is published is a noisy observation of `C(t)`
    ## on that day.
    cap_days = capacity_history.days
    cap_modelled = [C[clamp(Int(d), 1, n)] for d in cap_days]
    cap_obs = isempty(capacity_history.counts) ? missing :
              collect(Int.(capacity_history.counts))
    bed_capacity ~ to_submodel(
        vintage_increments_model(cap_modelled, cap_obs, k))

    ## Cut-off occupancy (beds in use), demand (need under unconstrained
    ## supply), the shortfall between them, the utilisation and the true-BVD
    ## share of demand.
    z0 = zero(eltype(C))
    occ_T = isempty(occupancy) ? z0 : occupancy[end]
    dem_T = isempty(demand) ? z0 : demand[end]
    bvd_dem_T = isempty(bvd_demand) ? z0 : bvd_demand[end]
    expected_isolation := safe_rate(occ_T)
    expected_bed_demand := safe_rate(dem_T)
    bed_shortfall := safe_rate(max(dem_T - occ_T, z0))
    bed_utilisation := safe_rate(occ_T) / safe_rate(C_T)
    isolation_bvd_share := safe_rate(bvd_dem_T) / safe_rate(dem_T)
    isolation_severity := sev_state.δ_iso
    isolation_bvd_admission := p_iso_bvd

    return (; p_iso, p_iso_bvd, δ_iso = sev_state.δ_iso,
        capacity = C_T, bvd_los_mean = bvd_los_state.mean,
        ruleout_los_mean = ruleout_los_state.mean,
        k_isolation = k, demand, occupancy, isolation,
        expected_isolation = safe_rate(occ_T),
        expected_bed_demand = safe_rate(dem_T))
end

Suspected deaths

Suspected deaths are the ascertained, CFR-weighted convolution of the daily onsets with the onset-to-death PMF $f_d$, plus a non-BVD background, modelled on the incidence scale. The death history ends at the cut-off, so the cut-off total is the final increment and is not scored separately. A fatal BVD infection enters the suspected-death count only when ascertained, so the BVD deaths carry a death ascertainment $p_{\text{death}}$, the death analogue of the case ascertainment $p_{\text{DRC}}$, with an informative prior centred high (a death is more reliably reported than a living suspect). The non-BVD background suspected deaths are a background CFR ${\mathrm{cfr}}_{\text{bg}}$ applied to the per-day non-BVD suspected-case background ${\lambda }_{\text{bg},t}$, lagged by the same onset-to-death delay so a background death follows its background case; the death background tracks the identified case background rather than a second free, outbreak-size- degenerate rate. The daily death series is

$$m_t=p_{\text{death}}\mathrm{CFR}\sum _{s\ge 0}{\text{onsets}}_{t-s}{f}_{d,s}+{\mathrm{cfr}}_{\text{bg}}\sum _{s\ge 0}{\lambda }_{\text{bg},t-s}{f}_{d,s}.$$

The per-vintage increments are scored with a NegBinomial sharing the dispersion $k$:

$$\begin{array}{}\text{(26)}& Y_{\text{deaths},i}-Y_{\text{deaths},i-1}\sim \mathrm{NegBinomial}(\sum _{t=d_{i-1}+1}^{d_i}{m}_t,k).\end{array}$$

Submodel: deaths_model

@model function deaths_model(
        deaths_history,
        total_deaths::Union{Missing, Integer},
        onsets::AbstractVector, k::Real;
        suspected_daily_deaths_history = (; days = Int[], counts = Int[]),
        cfr = cfr_model(),
        ascertainment = death_ascertainment_model(),
        case_bg_daily = nothing,
        background_cfr = background_cfr_model(),
        death_background = nothing,
        background_re = nothing,
        ## nmax covers 98% of the convolved onset->death sum (the two atomic
        ## Gammas moment-matched to a single Gamma only for the truncation).
        onset_to_death = onset_to_death_model(cdf_nmax(Gamma(3.33, 3.83));
            oa_alpha_prior = truncated(Normal(1.178, 0.285); lower = 0.01),
            oa_theta_prior = truncated(Normal(3.694, 1.198); lower = 0.1),
            ad_alpha_prior = truncated(Normal(2.151, 0.604); lower = 0.01),
            ad_theta_prior = truncated(Normal(3.906, 1.381); lower = 0.1)))
    cfr_state ~ to_submodel(cfr)
    od_state ~ to_submodel(onset_to_death)
    asc_state ~ to_submodel(ascertainment)
    CFR = cfr_state.CFR
    p_death = asc_state.p_death
    bvd_deaths_daily = (p_death * CFR) .* convolve_delay(onsets, od_state.pmf)

    n = length(bvd_deaths_daily)
    vobs = vintage_obs(deaths_history, total_deaths, n)

    ## Daily non-BVD background deaths: the background CFR `cfr_bg` applied to
    ## the non-BVD suspected-case background `case_bg_daily`, lagged by the
    ## onset-to-death delay so a background death follows its background case
    ## the way the BVD deaths follow the onsets. `λ_bg_death` is the mean daily
    ## background death rate. The `background_re`, `death_background` and
    ## pure-BVD branches are sensitivity fallbacks.
    if case_bg_daily !== nothing
        bgcfr_state ~ to_submodel(background_cfr)
        cfr_bg = bgcfr_state.cfr_bg
        bg_death_daily = cfr_bg .* convolve_delay(case_bg_daily, od_state.pmf)
        λ_bg_death = sum(bg_death_daily) / n
        bg_death_sigma = zero(CFR)
    elseif background_re !== nothing
        bg_state ~ to_submodel(background_re(n))
        cfr_bg = zero(CFR)
        λ_bg_death = bg_state.λ_mu
        bg_death_sigma = bg_state.σ_bg
        bg_death_daily = bg_state.λ
    elseif death_background !== nothing
        dbg_state ~ to_submodel(death_background)
        cfr_bg = zero(CFR)
        λ_bg_death = dbg_state.λ_bg_death
        bg_death_sigma = zero(λ_bg_death)
        bg_death_daily = fill(λ_bg_death, n)
    else
        cfr_bg = zero(CFR)
        λ_bg_death = zero(CFR)
        bg_death_sigma = zero(CFR)
        bg_death_daily = fill(zero(CFR), n)
    end

    deaths_daily = bvd_deaths_daily .+ bg_death_daily

    modelled_increments = bin_increments(deaths_daily, vobs.days)
    death_increments ~ to_submodel(
        vintage_increments_model(modelled_increments, vobs.obs_increments, k))

    ## Daily new suspected deaths ("cas suspects du jour N (M deces)"): per-day
    ## counts scored against the modelled daily suspected-death series at each
    ## report day. The mean for day `d` is the single-day `deaths_daily[d]`
    ## (clamped into the grid), NOT a between-vintage increment — this is a
    ## genuine daily count, so it never differences a falling cumulative. Empty
    ## by default; a `missing` count vector samples (the predictive path). The
    ## deaths analogue of the suspected-case daily inflow.
    sdd_days = suspected_daily_deaths_history.days
    sdd_modelled = [deaths_daily[clamp(Int(d), 1, n)] for d in sdd_days]
    sdd_obs = isempty(suspected_daily_deaths_history.counts) ? missing :
              collect(Int.(suspected_daily_deaths_history.counts))
    suspected_daily_deaths ~ to_submodel(
        vintage_increments_model(sdd_modelled, sdd_obs, k))

    raw_total = sum(deaths_daily)
    expected_deaths_T := safe_rate(raw_total)
    bg_death_total = sum(bg_death_daily)

    return (; CFR, p_death, cfr_bg, od_pmf = od_state.pmf, deaths_daily,
        bvd_deaths_daily, expected_deaths_T, λ_bg_death, bg_death_sigma,
        bg_death_daily, bg_death_total)
end

Laboratory pipeline

The laboratory pipeline fits a single analysed-specimen volume. There is no separately-modelled testing capacity: the analysed volume is a deterministic function of the suspected-case incidence. It is the suspected daily pipeline ($p_{\text{DRC}}{\text{bvd}}_t$ plus the non-BVD background ${\lambda }_{\text{bg}}$) carried through the report-to-analysed delay $f_{\text{rec}}$ and thinned by the testing fraction ${\tau }_{\text{test}}$ (the share of suspected cases routed to the laboratory),

$$v_t={\tau }_{\text{test}}\sum _{s\ge 0}(p_{\text{DRC}}{\text{bvd}}_{t-s}+{\lambda }_{\text{bg},t-s})f_{\text{rec},s}.$$

This analysed volume is gated to zero before the testing onset: no specimens are analysed before the laboratory existed, so $v_t$ does not accrue over the pre-surveillance cryptic phase (modelling a pre-testing volume would both invent capacity and roll it into the first laboratory and early-confirmed bins, over-predicting the early confirmed counts). The first confirmed vintage is treated as the baseline and the early confirmed increments are scored from it. The suspected-case count itself is not gated, as those cases did accumulate over the cryptic phase.

This construction, a testing fraction times the suspected pipeline carried to laboratory receipt, gives the modelled case analysed volume that the confirmed deaths reuse: the death volume scales it at the per-day suspected death-to-case ratio (see the confirmed-deaths section below), so the two share the laboratory capacity onset.

The per-vintage increments are scored against the cumulative analysed series with a NegBinomial sharing the dispersion $k$:

$$\begin{array}{}\text{(27)}& Y_{\text{ana},i}-Y_{\text{ana},i-1}\sim \mathrm{NegBinomial}(\sum _{t=d_{i-1}+1}^{d_i}{v}_t,k).\end{array}$$

The confirmed positives in each laboratory window $v$ are scored as a Binomial of the observed specimens-analysed denominator $A_v$ with a per-window tested-positive probability $p_{\text{pos},v}$, and where no analysed count is observed (the early and unanchored windows) the modelled volume $v_t$ is the denominator instead, so the fitted volume and the proxy denominator are the same quantity. We tie that probability to the composition of the tested pool, so the confirmed data help identify the non-BVD background. The suspect-pool composition ${\phi }_v$ is the BVD share among the specimens analysed in the window, carried through the same delay as the volume so composition and volume share one clock:

$${\phi }_v=\frac{(p_{\text{DRC}}\text{bvd}\ast f_{\text{rec}}{)}_v}{(p_{\text{DRC}}\text{bvd}\ast f_{\text{rec}}{)}_v+({\lambda }_{\text{bg}}\ast f_{\text{rec}}{)}_v}.$$

The tested BVD share $q_v$ raises ${\phi }_v$ by the decaying severity enrichment ${\delta }_0$. A truly BVD specimen then tests positive with the sensitivity $s$, and a non-BVD specimen with the false-positive rate $1-\mathrm{spec}$, so the false-positive term carries the non-BVD share and the laboratory data identify the background:

$$q_v=\mathrm{logistic}(\mathrm{logit}({\phi }_v)+{\delta }_0{e}^{-c_v/\text{decay}}),p_{\text{pos},v}=s{q}_v+(1-\mathrm{spec})(1-q_v),$$$$\begin{array}{}\text{(28)}& C_v\sim \mathrm{Binomial}(A_v,p_{\text{pos},v}),\end{array}$$

with $c_v$ the cumulative modelled laboratory volume at window $v$, the clock on which the enrichment decays. The confirmed vintages before the first and after the last laboratory date carry no observed analysed denominator. They are scored as NegBinomial counts against the modelled laboratory volume $V_v$, the daily modelled volume $v_t$ summed over the window, with the same composition-linked positivity, so all the confirmed data are used:

$$\begin{array}{}\text{(29)}& C_v^{\text{no-denom}}\sim \mathrm{NegBinomial}(p_{\text{pos},v}{V}_v,k).\end{array}$$

Submodel: lab_delay_model (receipt delay)

@model function lab_delay_model(nmax::Integer = cdf_nmax(lognormal_meansd(4.5, 4.0));
        mean_prior = truncated(Normal(4.5, 1.0); lower = 1),
        sd_prior = truncated(Normal(4.0, 0.75); lower = 1))
    d ~ to_submodel(censored_delay_model(nmax; mean_prior, sd_prior))
    return (; pmf = d.pmf, dist = d.dist, mean = d.mean, sd = d.sd)
end

Submodel: confirmed_cases_model

@model function confirmed_cases_model(
        confirmed_history,
        confirmed_cases::Union{Missing, Integer},
        onsets::AbstractVector, k::Real, p_drc::Real,
        bg_daily::AbstractVector, τ_test::Real,
        bvd_reports_daily::AbstractVector;
        lab_history = (; days = Int[], counts = Int[]),
        lab_daily_history = (; days = Int[], counts = Int[]),
        tests_analysed::Union{Missing, Integer} = missing,
        receipt = lab_delay_model(),
        positivity = confirmed_positivity_model,
        positivity_link::Symbol = :composition,
        severity_enrichment = severity_enrichment_model(),
        sensitivity = test_sensitivity_model(),
        specificity = test_specificity_model(),
        ## When false, the early/late windows (confirmed vintages with NO
        ## observed analysed denominator) are not scored — only the
        ## observed-denominator Binomial windows contribute, so confirmed
        ## informs positivity without extrapolating a denominator from
        ## incidence. Used to probe the no-test-data extrapolation.
        fit_unanchored::Bool = true)
    n = length(onsets)
    ## `missing` cut-off scalar means generator mode: observed increments are
    ## left missing so `predict` resamples them.
    have_data = !ismissing(confirmed_cases)

    ## Laboratory capacity onset. No specimens are analysed before testing
    ## existed, so the modelled analysed volume is gated to zero before the
    ## first confirmed-case vintage (the earliest evidence of testing; the
    ## first laboratory date is the fallback). Modelling a pre-testing analysed
    ## volume would invent capacity that did not exist AND roll it into the
    ## first laboratory and early-confirmed bins, vastly over-predicting the
    ## early confirmed counts. The suspected-case pipeline feeding the volume is
    ## NOT gated — suspected cases did accumulate over the cryptic phase.
    cap_start = !isempty(confirmed_history.days) ?
                clamp(Int(confirmed_history.days[1]), 1, n) :
                (!isempty(lab_history.days) ?
                 clamp(Int(lab_history.days[1]), 1, n) : 1)

    ## Analysed-specimen volume: the suspected pipeline carried through the
    ## report-to-analysed delay and thinned by the tested fraction, fit to the
    ## analysed series and reused as the denominator in the early and unanchored late
    ## windows below. `bg_daily` is the per-day non-BVD background.
    receipt_state ~ to_submodel(receipt)
    suspected_daily = p_drc .* bvd_reports_daily .+ bg_daily
    analysed_daily = τ_test .* convolve_delay(suspected_daily,
        receipt_state.pmf)
    ## In predict mode (no AD) the daily series can infer as `Vector{Any}`
    ## on some Julia versions, which then makes `reduce_empty` / `zero(Any)`
    ## fail on the empty derived window vectors below. Concretise to the
    ## working scalar type; this runs only when the element type has widened,
    ## so the AD/fit path (concrete eltype) is left untouched.
    if eltype(analysed_daily) === Any
        analysed_daily = convert(Vector{typeof(τ_test)}, analysed_daily)
    end
    ## Gate the capacity to the testing window: zero before `cap_start`.
    analysed_daily = gate_before(analysed_daily, cap_start)
    rvobs = vintage_obs(lab_history, tests_analysed, n)
    analysed_inc = bin_increments(analysed_daily, rvobs.days)
    ## Generator mode leaves the volume increments missing so `predict`
    ## resamples them, like the early/late windows below.
    vol_obs = have_data ? rvobs.obs_increments : missing
    analysed_increments ~ to_submodel(
        vintage_increments_model(analysed_inc, vol_obs, k))

    ## Post-cutoff 24h analysed volume. After the national cumulative analysed
    ## series stops, INSP publishes a 24h analysed count on some days; the
    ## modelled daily analysed volume on each such day is scored against that
    ## count, so the post-cutoff testing throughput is fitted from the same
    ## stream rather than only used as a confirmed denominator. Same
    ## `have_data` gate so `predict` resamples it.
    daily_days = [clamp(Int(d), 1, n) for d in lab_daily_history.days]
    daily_modelled = isempty(daily_days) ? similar(analysed_daily, 0) :
                     [analysed_daily[d] for d in daily_days]
    daily_obs = have_data ? lab_daily_history.counts : missing
    analysed_daily_increments ~ to_submodel(
        vintage_increments_model(daily_modelled, daily_obs, k))

    ## Confirmed positives in three groups sharing one partially-pooled
    ## positivity: early windows (before the first lab date, no observed
    ## analysed) scored as counts against the modelled laboratory volume,
    ## observed windows scored as a Binomial of the observed analysed
    ## denominator, and late windows (after the last lab date, INSP's
    ## confirmed-only format) scored as counts against the modelled volume
    ## like the early windows.
    windows = confirmed_positivity_windows(confirmed_history, lab_history,
        lab_daily_history)
    n_early = length(windows.early_days)
    n_obs = length(windows.obs_analysed)
    n_late = length(windows.late_days)
    nv = n_early + n_obs + n_late

    ## Per-window tested BVD share `p_pos`. Two links:
    ## `:free` — a free partially-pooled per-window random effect
    ## ([`confirmed_positivity_model`](@ref)), decoupled from `λ_bg`.
    ## `:composition` (default) — the tested share is the suspect-pool
    ## `φ_v = (p_drc·BVD)_v / ((p_drc·BVD)_v + λ_bg_v)` over each laboratory
    ## window, upsampled by a decaying severity enrichment δ0 (see
    ## [`severity_enrichment_model`](@ref)), so the lab positivity identifies
    ## the background `λ_bg` rather than absorbing it into a free curve.
    window_days = vcat(windows.early_days, windows.obs_days,
        windows.late_days)
    if positivity_link === :composition
        enrich_state ~ to_submodel(severity_enrichment, false)
        δ0 = enrich_state.δ0
        decay_scale = enrich_state.decay_scale
        ## PCR sensitivity and specificity. The tested-positive probability
        ## is `p = s · q + (1 − spec)(1 − q)` with `q` the tested BVD share:
        ## the false-positive term `(1 − spec)(1 − q)` makes the confirmed
        ## counts respond to the non-BVD share `1 − q`, so the laboratory
        ## data identify the background `λ_bg` through the composition `φ`
        ## rather than the BVD signal alone. Without it the confirmed
        ## positivity tracks only `q`, leaving `λ_bg` weakly identified.
        sens_state ~ to_submodel(sensitivity, false)
        spec_state ~ to_submodel(specificity, false)
        s_test = sens_state.s_test
        spec = spec_state.spec
        ## Suspect-pool composition over each window, carried through the
        ## report-to-analysed delay so it reflects the composition of the
        ## specimens actually analysed in the window, consistent with the
        ## modelled volume `analysed_daily`. The `τ_test` factor cancels in the
        ## ratio φ, so it is omitted here.
        analysed_bvd_daily = convolve_delay(p_drc .* bvd_reports_daily,
            receipt_state.pmf)
        analysed_bg_daily = convolve_delay(bg_daily, receipt_state.pmf)
        if eltype(analysed_bvd_daily) === Any
            analysed_bvd_daily = convert(Vector{typeof(τ_test)},
                analysed_bvd_daily)
            analysed_bg_daily = convert(Vector{typeof(τ_test)},
                analysed_bg_daily)
        end
        ## Gate the tested composition to the testing window too, so the
        ## composition clock and the per-window BVD share start at the testing
        ## onset rather than rolling the cryptic phase.
        analysed_bvd_daily = gate_before(analysed_bvd_daily, cap_start)
        analysed_bg_daily = gate_before(analysed_bg_daily, cap_start)
        bvd_window = bin_increments(analysed_bvd_daily, window_days)
        bg_window = bin_increments(analysed_bg_daily, window_days)
        Tt = eltype(bvd_window)
        ## Testing clock: cumulative modelled analysed volume at each window.
        vol_window = bin_increments(analysed_daily, window_days)
        c_window = cumsum(vol_window)
        lo = convert(Tt, 1e-8)
        hi = one(Tt) - lo
        ## Floor the decay scale so a near-zero `decay_scale` draw cannot make
        ## the clock ratio `0/0` (NaN) and break the downstream Binomial.
        dscale = max(convert(Tt, decay_scale), one(Tt))
        s_t = convert(Tt, s_test)
        sp_t = convert(Tt, spec)
        p_pos = map(eachindex(window_days)) do i
            ## Pool composition φ = (p_drc·BVD) / ((p_drc·BVD) + λ_bg) over the
            ## window, guarded against a zero/negative denominator.
            num = bvd_window[i]
            den = bvd_window[i] + bg_window[i]
            ratio = num / (den + lo)
            φ = clamp(isfinite(ratio) ? ratio : convert(Tt, 0.5), lo, hi)
            δ_i = convert(Tt, δ0) * exp(-c_window[i] / dscale)
            ## Severity-enriched tested BVD share, then the assay
            ## sensitivity/specificity transform to the tested-positive
            ## probability so the false-positive term identifies `λ_bg`.
            q = logistic(logit(φ) + δ_i)
            qf = isfinite(q) ? q : φ
            qe = clamp(qf, lo, hi)
            p = s_t * qe + (one(Tt) - sp_t) * (one(Tt) - qe)
            ## Final guard: clamp into (0,1) and replace any non-finite value
            ## with the composition so the confirmed Binomial always sees a
            ## valid probability even under an AD perturbation.
            clamp(isfinite(p) ? p : φ, lo, hi)
        end
    else
        pos_state ~ to_submodel(positivity(nv))
        p_pos = pos_state.p_pos
    end

    ## Early windows: confirmed increment ~ NegBinomial(positivity ×
    ## modelled analysed volume), the volume binned over each window's OWN day
    ## range pinned at `early_start` (the first confirmed vintage, the testing-
    ## onset baseline), so the first early increment is scored from the data
    ## start rather than rolling the (now-gated) pre-testing volume. Mirrors the
    ## late-window pinning at `late_start`.
    early_p = p_pos[1:n_early]
    early_volume = n_early > 0 ?
                   bin_increments(analysed_daily,
        vcat(windows.early_start, windows.early_days))[2:end] :
                   similar(analysed_daily, 0)
    early_mean = early_p .* early_volume
    early_obs = (have_data && n_early > 0 && fit_unanchored) ?
                windows.early_increments : missing
    early_increments ~ to_submodel(
        vintage_increments_model(early_mean, early_obs, k))

    ## Observed windows: Binomial of the observed analysed denominator.
    obs_p = p_pos[(n_early + 1):(n_early + n_obs)]
    obs_positives = (have_data && n_obs > 0) ? collect(windows.obs_positives) :
                    missing
    confirmed_positives ~ to_submodel(
        confirmed_positives_model(obs_positives, windows.obs_analysed, obs_p))

    ## Late windows: confirmed-only vintages after the last laboratory date.
    ## A day that publishes a 24h analysed count (`late_analysed > 0`) is
    ## scored as a Binomial of that observed denominator — like an observed
    ## window, anchoring its positivity to data — and each remaining unanchored day
    ## as NegBinomial(positivity × modelled volume). The modelled volume is
    ## binned over each late window's own day range, with the running edge
    ## PINNED at the last laboratory day (`late_start`): `bin_increments`
    ## runs its running `prev` from day 0, so prepending `late_start` to the
    ## late day edges and dropping the synthetic first bin starts the
    ## accumulation at `late_start`, avoiding double-counting the
    ## observed-window volume.
    late_p = p_pos[(n_early + n_obs + 1):nv]
    if n_late > 0
        late_edges = vcat(windows.late_start, windows.late_days)
        late_volume = bin_increments(analysed_daily, late_edges)[2:end]
    else
        late_volume = similar(analysed_daily, 0)
    end
    late_mean = late_p .* late_volume
    ## Observed late increments: anchored days (24h denominator) carry the
    ## confirmed increment clamped into the Binomial support and are always
    ## scored; unanchored days are scored only when `fit_unanchored` (the
    ## no-extrapolation probe leaves them latent). A per-entry
    ## `missing`/value vector lets the one submodel observe each accordingly.
    if have_data && n_late > 0
        late_obs = Vector{Union{Missing, Int}}(undef, n_late)
        for i in 1:n_late
            a = windows.late_analysed[i]
            if a > 0
                late_obs[i] = clamp(windows.late_increments[i], 0, a)
            elseif fit_unanchored
                late_obs[i] = windows.late_increments[i]
            else
                late_obs[i] = missing
            end
        end
    else
        late_obs = missing
    end
    late_increments ~ to_submodel(
        late_confirmed_model(late_obs, late_mean, windows.late_analysed,
        late_p, k))

    expected_analysed := safe_rate(sum(analysed_daily))
    ## Expected confirmed at the cut-off and the overall positivity, over the
    ## modelled early volume, the observed cumulative analysed windows and the
    ## late windows (anchored days contribute `p · analysed`, unanchored days the
    ## modelled `p · volume`). The window vectors are empty when a vintage has
    ## no such window, and in predict mode their element type can widen to
    ## `Any`, so each sum is given a concrete `init` to skip `reduce_empty`'s
    ## `zero(Any)`. The init is taken from the scalar `τ_test` (always
    ## concrete), NOT from `eltype(p_pos)`, which can widen to `Any`.
    z = zero(τ_test)
    amask = windows.late_analysed .> 0
    late_den_a = float.(windows.late_analysed)
    late_den = n_late > 0 ? ifelse.(amask, late_den_a, late_volume) :
               similar(late_volume, 0)
    late_expected = n_late > 0 ? ifelse.(amask, late_p .* late_den_a, late_mean) :
                    similar(late_mean, 0)
    denom = sum(early_volume; init = z) + float(sum(windows.obs_analysed)) +
            sum(late_den; init = z)
    expected_positives = sum(early_mean; init = z) +
                         sum(late_expected; init = z) +
                         (n_obs > 0 ? sum(obs_p .* windows.obs_analysed) : z)
    expected_confirmed := safe_rate(expected_positives)
    p_positive := safe_rate(expected_positives) / safe_rate(denom)

    ## Modelled daily confirmed-case incidence: the per-window tested-positive
    ## probability expanded onto the daily grid times the modelled analysed
    ## volume. Exposed so the composer's cumulative-confirmed trajectory and a
    ## survivors-among-confirmed (recovered) stream can reuse one consistent
    ## daily series. In predict / check-model mode `p_pos` can widen to
    ## `Vector{Any}`, so pin it to the analysed volume's (always-concrete)
    ## element type before expanding, matching the other guards here.
    p_pos_daily = p_pos
    if eltype(p_pos_daily) === Any
        p_pos_daily = convert(Vector{eltype(analysed_daily)}, p_pos_daily)
    end
    confirmed_daily = expand_vintage_rate(p_pos_daily, window_days, n) .*
                      analysed_daily

    return (; τ_test, bg_daily, p_pos, windows, analysed_daily,
        confirmed_daily,
        receipt_pmf = receipt_state.pmf,
        expected_analysed, expected_confirmed, p_positive)
end

Confirmed deaths

The confirmed deaths mirror the confirmed-case laboratory pipeline. The confirmed cases fit a modelled analysed-specimen volume and score the positives as that volume times a composition-linked positivity; the death side has no published analysed denominator, so we build the death analogue of that volume and score the confirmed-death increments as NegBinomial counts of it.

Deaths are tested out of the same laboratory as cases, so the death analysed volume tracks the modelled case analysed volume $v_t^{\text{c}}$ at the per-day suspected death-to-case ratio, times a testing-intensity scaling,

$$v_t^{\text{d}}=\text{scaling}{v}_t^{\text{c}}\frac{\sum _{s\ge 0}{m}_{t-s}^{\text{d}}{f}_{\text{rec},s}}{\sum _{s\ge 0}{m}_{t-s}^{\text{c}}{f}_{\text{rec},s}},$$

with $m^{\text{d}}$ and $m^{\text{c}}$ the modelled suspected-death and suspected-case series and $f_{\text{rec}}$ the report-to-receipt delay the confirmed cases use. The death-to-case ratio carries the suspect-pool severity and the suspected-death level, so the scaling is the per-suspect testing-intensity difference between deaths and living suspects alone; with no death-testing data it is a tight log-normal centred on one. Those specimens confirm at the assay positivity built from the death-pool BVD share

$$q_{\text{death},t}=\frac{{\text{bvd}}_t^{\text{d}}}{{\text{bvd}}_t^{\text{d}}+{\text{bg}}_t^{\text{d}}},p_t=s{q}_{\text{death},t}+(1-\mathrm{spec})(1-q_{\text{death},t}),$$

with ${\text{bvd}}^{\text{d}}$ and ${\text{bg}}^{\text{d}}$ the BVD and non-BVD components of the suspected deaths (both at receipt) and $s$, $\mathrm{spec}$ the same assay sensitivity and specificity as the confirmed cases. The false-positive term $(1-\mathrm{spec})(1-q_{\text{death}})$ makes the confirmed deaths respond to the non-BVD death share, the same structural link the confirmed cases use; the death background (the background CFR applied to the case background, lagged by the onset-to-death delay) keeps the composition below one. The daily confirmed deaths are the positivity times the death analysed volume,

$${\text{cd}}_t=p_t{v}_t^{\text{d}},$$

and the per-vintage increments are scored with a NegBinomial sharing the dispersion $k$:

$$\begin{array}{}\text{(30)}& Y_{\text{cd},i}-Y_{\text{cd},i-1}\sim \mathrm{NegBinomial}(\sum _{t=d_{i-1}+1}^{d_i}{\text{cd}}_t,k).\end{array}$$

The death analysed volume inherits the laboratory capacity onset from the case volume $v_t^{\text{c}}$, so ${\text{cd}}_t$ is zero before the first confirmed-case vintage: no deaths are confirmed before the laboratory existed.

Submodel: confirmed_deaths_model

@model function confirmed_deaths_model(
        confirmed_deaths::Union{Missing, Integer},
        total_deaths::Union{Missing, Integer},
        deaths_daily::AbstractVector,
        bvd_deaths_daily::AbstractVector,
        bg_death_daily::AbstractVector, k::Real;
        confirmed_deaths_history = (; days = Int[], counts = Int[]),
        receipt_pmf::AbstractVector = [1.0],
        capacity_start::Integer = 0,
        case_analysed_daily = nothing,
        case_suspected_daily = nothing,
        scaling = death_testing_scaling_model(),
        testing = death_testing_fraction_model(),
        sensitivity = test_sensitivity_model(),
        specificity = test_specificity_model())
    sens_state ~ to_submodel(sensitivity)
    spec_state ~ to_submodel(specificity)
    s = sens_state.s_test
    spec = spec_state.spec
    n = length(deaths_daily)

    ## Suspected deaths carried to laboratory receipt by the same
    ## report-to-receipt delay the confirmed cases use, with the BVD component.
    susp_death = convolve_delay(deaths_daily, receipt_pmf)
    bvd_death = convolve_delay(bvd_deaths_daily, receipt_pmf)
    ## In predict or check-model mode the series can widen to `Vector{Any}`,
    ## which trips `zero(Any)` downstream; pin to the sampled scalar type,
    ## leaving the fit path (concrete dual eltype) untouched.
    if eltype(susp_death) === Any
        susp_death = convert(Vector{typeof(s)}, susp_death)
        bvd_death = convert(Vector{typeof(s)}, bvd_death)
    end

    ## Death-pool BVD composition per day, q = bvd / (bvd + bg), and the assay
    ## tested-positive probability p = s·q + (1 − spec)(1 − q). The false-
    ## positive term lets the confirmed deaths respond to the non-BVD share,
    ## the same link the confirmed cases use.
    lo = eps(typeof(s))
    hi = one(s) - lo
    q_death_daily = map(eachindex(susp_death)) do t
        den = susp_death[t]
        ratio = den > lo ? bvd_death[t] / den : one(s)
        clamp(isfinite(ratio) ? ratio : one(s), lo, hi)
    end
    p_pos_daily = s .* q_death_daily .+ (one(s) - spec) .*
                                        (one(s) .- q_death_daily)

    ## Death analysed volume. Deaths are tested out of the same laboratory as
    ## cases, so the death volume tracks the modelled case analysed volume at
    ## the per-day suspected death-to-case ratio, times a testing-intensity
    ## scaling, v = scaling · analysed_case · susp_death / susp_case. The case
    ## volume already carries the laboratory capacity onset, so the death volume
    ## inherits it. The death-only composer has no case stream and falls back to
    ## a death testing fraction of the suspected deaths, gated at the onset.
    if case_analysed_daily !== nothing
        scale_state ~ to_submodel(scaling)
        sc = scale_state.scaling
        susp_case = convolve_delay(case_suspected_daily, receipt_pmf)
        death_volume = map(eachindex(susp_death)) do t
            den = susp_case[t]
            den > lo ? sc * case_analysed_daily[t] * susp_death[t] / den :
            zero(sc)
        end
        τ_death = susp_death[n] > lo ? death_volume[n] / susp_death[n] : zero(sc)
    else
        test_state ~ to_submodel(testing)
        τ_death = test_state.τ_death
        sc = one(τ_death)
        death_volume = τ_death .* gate_before(susp_death, capacity_start)
    end

    confirmed_death_daily = p_pos_daily .* death_volume
    vobs = vintage_obs(confirmed_deaths_history, confirmed_deaths, n)
    modelled_inc = bin_increments(confirmed_death_daily, vobs.days)
    cdeath_increments ~ to_submodel(
        vintage_increments_model(modelled_inc, vobs.obs_increments, k))

    expected_confirmed_deaths := safe_rate(sum(confirmed_death_daily))
    ## Cut-off death-pool composition and confirmation positivity, surfaced as
    ## `death_composition` and `death_confirmation`.
    q_death := q_death_daily[n]
    p_death_conf := p_pos_daily[n]

    return (; τ_death, scaling = sc, s_test = s, spec, q_death, p_death_conf,
        expected_confirmed_deaths)
end

Recovered among confirmed

Recoveries ("cumul guéris") are the survivors among laboratory-confirmed cases, the incidence analogue of the convolution-and-scaling secondary-observation model of EpiNow2 (Abbott et al., 2020). The modelled daily confirmed incidence ${\text{confirmed}}_t$ (the per-window tested-positive probability on the modelled analysed volume, the same daily series the cumulative-confirmed trajectory uses) is scaled by the recovery proportion $p_{\text{rec}}$ and convolved with a sampled confirmation-to-recovery delay $f_{\text{rec}}$,

$${\text{recovered}}_t=p_{\text{rec}}\sum _{s\ge 0}{\text{confirmed}}_{t-s}{f}_{\text{rec},s}.$$

A recovered case is one that did not die, so the recovery proportion is grounded on the case-fatality ratio rather than estimated independently. It is the complement $1-\mathrm{CFR}$ adjusted on the log-odds scale by a sampled offset ${\delta }_{\text{rec}}\sim \mathrm{Normal}(0,0.5)$, since the confirmed cases are a slightly different population from the one the CFR is defined over,

$$p_{\text{rec}}=\mathrm{logistic}(\mathrm{logit}(1-\mathrm{CFR})+{\delta }_{\text{rec}}).$$

A case is taken to be confirmed before it is recorded as recovered (the report counts recoveries among confirmed cases); a positive result could in principle return after a patient has already recovered, but we assume the reported total reflects confirmed cases recorded as recovered. The cumulative recovered series ends at the cut-off, so its per-vintage increments are fitted, like the confirmed and confirmed-death streams, with a NegBinomial whose dispersion $k_{\text{rec}}$ is its own rather than shared with the other streams:

$$Y_{\text{rec},i}-Y_{\text{rec},i-1}\sim \mathrm{NegBinomial}(\sum _{t=d_{i-1}+1}^{d_i}{\text{recovered}}_t,k_{\text{rec}}).$$

The convolution right-censors recoveries that have not yet resolved by the cut-off, so the small observed totals (12 to 40 over 6-13 June) are consistent with a high eventual survival fraction and a multi-week recovery delay.

Submodel: recovered_model

@model function recovered_model(
        recovered_history,
        recovered_total::Union{Missing, Integer},
        confirmed_daily::AbstractVector, CFR::Real;
        recovery = recovery_probability_model,
        dispersion = surveillance_dispersion_model(),
        ## Confirmation-to-recovery (discharge) delay; an Ebola survivor is
        ## discharged a couple of weeks after confirmation, so the default is
        ## a mean ~14 d stay before recovery is recorded.
        confirmation_to_recovery = censored_delay_model(
            cdf_nmax(lognormal_meansd(14.0, 8.0); q = 0.99);
            mean_prior = truncated(Normal(14.0, 5.0); lower = 1),
            sd_prior = truncated(Normal(8.0, 4.0); lower = 1)),
        ## Dispersion can be injected from the joint composer's pooled set
        ## (`k_external`); standalone it samples its own from `dispersion`.
        k_external::Union{Nothing, Real} = nothing)
    ## Recovery fraction grounded on the CFR complement (see
    ## `recovery_probability_model`), adjusted for the confirmed population.
    rec_state ~ to_submodel(recovery(CFR))
    p_recover = rec_state.p_recover
    if k_external === nothing
        disp_state ~ to_submodel(dispersion)
        k = disp_state.k
    else
        k = k_external
    end
    delay_state ~ to_submodel(confirmation_to_recovery)

    ## Survivors among confirmed cases, lagged by the confirmation-to-recovery
    ## delay: a scaled convolution of the modelled daily confirmed incidence.
    recovered_daily = p_recover .* convolve_delay(confirmed_daily,
        delay_state.pmf)

    n = length(confirmed_daily)
    vobs = vintage_obs(recovered_history, recovered_total, n)
    modelled_inc = bin_increments(recovered_daily, vobs.days)
    recovered_increments ~ to_submodel(
        vintage_increments_model(modelled_inc, vobs.obs_increments, k))

    expected_recovered := safe_rate(sum(recovered_daily))

    return (; p_recover, recovery_delay_mean = delay_state.mean,
        k_recovered = k, recovered_daily, expected_recovered)
end

Exported cases

The exports stream is travel-gated, so the at-risk clock runs from infection. An infected person travels to Uganda at the daily per-capita travel rate $q=N_{\text{travel}}/N_{\text{source}}$ and stays at risk of being exported and detected only until the infection-to-detection delay has elapsed. The daily at-risk export prevalence is the infections still infected and not yet detected, scaled by the Uganda ascertainment and the travel rate. The infection-to-detection delay is the onset-to-detection delay convolved with the incubation period, so the survival clock runs from infection. Write the cumulative infections and the infections that have completed the detection delay as

$$C_t=\sum _{u\le t}{I}_u,{\text{det}}_t=\sum _{s\ge 0}{I}_{t-s}(f_{\text{inc}}\ast f_{\text{det}}{)}_s.$$

Then the daily export intensity and its running sum are

$$\begin{array}{}\text{(31)}& {\lambda }_t=p_{\text{Uganda}}q(C_t-{\text{det}}_t),\mathrm{\Lambda }(t)=\sum _{u\le t}{\lambda }_u.\end{array}$$

We model outbound travel only, not return, so this term would overestimate the infections on its own. Each observed Uganda import is fitted at its reported detection date. An import detected on a given day is scored as a Poisson of the rise in cumulative export intensity between consecutive detection dates, with a term before the earliest detection $d_1$ observed at zero, since no export is expected then. After the last detection date we stop modelling exports rather than scoring further zeros: travellers' reasons for crossing the border change over the outbreak, so the baseline travel rate no longer applies beyond it and the export clock is truncated there:

$$\begin{array}{}\text{(32)}& Y_{\text{exports},i}\sim \mathrm{Poisson}(\mathrm{\Lambda }(d_i)-\mathrm{\Lambda }(d_{i-1})),0\sim \mathrm{Poisson}(\mathrm{\Lambda }(d_1-1)).\end{array}$$

Submodel: exports_model

@model function exports_model(
        exported_cases::Union{Missing, Integer},
        infections::AbstractVector, p_uganda::Real;
        export_case_days::AbstractVector{<:Integer} = Int[],
        pre_detection_exports::Union{Missing, Integer} = 0,
        incubation_pmf::AbstractVector,
        source_population::Real = ITURI_POPULATION,
        traveller = traveller_volume_model(),
        ## Export detection abroad uses the same line-list onset→admission
        ## delay (d_oa) as the suspect-case report: a case is detected at a
        ## point of entry when first formally seen, ~4 days after onset.
        onset_to_detection = gamma_delay_model(cdf_nmax(Gamma(1.178, 3.694));
            alpha_prior = truncated(Normal(1.178, 0.285); lower = 0.01),
            theta_prior = truncated(Normal(3.694, 1.198); lower = 0.1)))
    travel_state ~ to_submodel(traveller)
    daily_travellers = travel_state.daily_travellers
    q = daily_travellers / source_population

    detect_state ~ to_submodel(onset_to_detection)
    ## Infection→detection delay: onset→detection convolved with the shared
    ## incubation PMF, so the survival clock runs from infection.
    f_det = convolve_pmf(incubation_pmf, detect_state.pmf)
    detected_daily = convolve_delay(infections, f_det)
    ## At-risk prevalence (person-days): infected but not yet detected.
    prevalence = cumsum(infections) .- cumsum(detected_daily)
    export_prevalence = p_uganda .* q .* prevalence
    n = length(export_prevalence)

    if isempty(export_case_days)
        ## No dated series: cumulative single-total Poisson at the cut-off.
        raw_exports = sum(export_prevalence)
        expected_exports_T := safe_rate(raw_exports)
        exported_cases ~ Poisson(expected_exports_T)
    else
        ## Dated per-day Poisson. The export clock stops at the last import
        ## `t_last` (the `last_offset` truncation); prevalence past it does
        ## not accrue. `d₁` is the earliest detection day.
        days, counts = dated_event_bins(export_case_days, n)
        d₁ = days[1]
        ## Pre-detection survival weight Λ(d₁−1): the cumulative export
        ## intensity up to the day before the earliest detection.
        pre = d₁ > 1 ? sum(@view export_prevalence[1:(d₁ - 1)]) :
              zero(@inbounds export_prevalence[begin])
        pre_detection_exports ~ Poisson(safe_rate(pre))
        ## Per-day-edge increments between consecutive detection days; the
        ## first is measured from the pre-detection weight `pre`, so the
        ## pre-detection term and the increments partition Λ(t_last).
        raw_inc = bin_increments(export_prevalence, days)
        μ_day = [i == 1 ? raw_inc[1] - pre : raw_inc[i]
                 for i in eachindex(raw_inc)]
        obs = ismissing(exported_cases) ? missing : counts
        export_obs ~ to_submodel(dated_poisson_model(μ_day, obs))
        ## Reported expected count is the cumulative intensity to `t_last`.
        expected_exports_T := safe_rate(pre + sum(μ_day))
    end

    ## Travel-scaled at-risk prevalence WITHOUT the export-case ascertainment
    ## `p_uganda`: a death among an exported case would be reported whether or
    ## not the case itself was ascertained as an import, so the export-death
    ## model accrues over the travelled person-time `q · prevalence`, not the
    ## ascertained `export_prevalence = p_uganda · q · prevalence`.
    travelled_prevalence = q .* prevalence
    return (; p_uganda, daily_travellers, q, prevalence,
        export_prevalence, travelled_prevalence,
        expected_exports = expected_exports_T)
end

Deaths among exports

The expected deaths among exports weight the travelled at-risk prevalence by the infection-to-death delay (the onset-to-death PMF convolved with the incubation period) and scale by the CFR. The travelled prevalence is the export prevalence before the ascertainment factor $p_{\text{Uganda}}$, because a death among an exported case would be reported whether or not the case itself was ascertained as an import. Writing it ${\ell }_t=q(C_t-{\text{det}}_t)$, the daily export-death intensity is

$${\mu }_t=\mathrm{CFR}\sum _{s\ge 0}{\ell }_{t-s}(f_{\text{inc}}\ast f_d{)}_s.$$

Its running sum is the cumulative export-death intensity:

$$\begin{array}{}\text{(33)}& {\mathrm{\Lambda }}_d(t)=\sum _{u\le t}{\mu }_u.\end{array}$$

Each dated Uganda export death is scored at its reported date with a per-day Poisson, the same dated-event likelihood the exports use, with a zero term before the first death day ${\delta }_1$:

$$\begin{array}{}\text{(34)}& Y_{\text{exp-deaths},i}\sim \mathrm{Poisson}({\mathrm{\Lambda }}_d({\delta }_i)-{\mathrm{\Lambda }}_d({\delta }_{i-1})),0\sim \mathrm{Poisson}({\mathrm{\Lambda }}_d({\delta }_1-1)).\end{array}$$

Submodel: exports_deaths_model

@model function exports_deaths_model(
        exports_deaths::Union{Missing, Integer},
        travelled_prevalence::AbstractVector, CFR::Real,
        od_pmf::AbstractVector, incubation_pmf::AbstractVector;
        export_death_days::AbstractVector{<:Integer} = Int[],
        pre_death_exports::Union{Missing, Integer} = 0)
    n = length(travelled_prevalence)
    ## Infection→death PMF by age (age 0 = same day).
    fd_pmf = convolve_pmf(incubation_pmf, od_pmf)
    ## Per-day expected export-death increment: CFR-scaled convolution of
    ## the daily at-risk prevalence with the infection→death PMF. Its
    ## running sum is the cumulative export-death intensity `Λ_d`.
    death_daily = CFR .* convolve_delay(travelled_prevalence, fd_pmf)

    if isempty(export_death_days)
        ## No dated series: cumulative single-total Poisson at the cut-off.
        expected_exports_deaths_T := safe_rate(sum(death_daily))
        exports_deaths ~ Poisson(expected_exports_deaths_T)
    else
        ## Dated per-day Poisson; the clock stops at the last death day.
        days, counts = dated_event_bins(export_death_days, n)
        δ₁ = days[1]
        pre = δ₁ > 1 ? sum(@view death_daily[1:(δ₁ - 1)]) :
              zero(@inbounds death_daily[begin])
        pre_death_exports ~ Poisson(safe_rate(pre))
        raw_inc = bin_increments(death_daily, days)
        μ_day = [i == 1 ? raw_inc[1] - pre : raw_inc[i]
                 for i in eachindex(raw_inc)]
        obs = ismissing(exports_deaths) ? missing : counts
        death_obs ~ to_submodel(dated_poisson_model(μ_day, obs))
        expected_exports_deaths_T := safe_rate(pre + sum(μ_day))
    end

    return (; expected_exports_deaths_T)
end

Joint model

The joint model runs the generating infection process once, stages it to daily onset incidence, and routes the shared onsets into every observation stream. It samples a single dispersion $k$ and the pooled ascertainment fractions, threading $p_{\text{DRC}}$ to the suspected-case, laboratory and confirmed-death likelihoods and $p_{\text{Uganda}}$ to the two Uganda-side likelihoods, and adds the genetic seeding bound on the outbreak age. Each observation stream argument may be dropped, so the same model structure generates prior- and posterior-predictive draws.

Alongside the joint model we write single-stream models for each count-based stream (exported cases, suspected deaths, suspected cases, laboratory-confirmed cases, confirmed deaths and deaths among exports), so each stream's posterior over the outbreak size can be compared with the joint. Other model variants reuse these models with different amounts of data, cutting the data to an earlier date or dropping the counts.

Composer: exports-only fit

@model function exports_only_model(
        n::Integer, exported_cases::Union{Missing, Integer};
        export_case_days::AbstractVector{<:Integer} = Int[],
        breakpoint::Union{Missing, Real} = missing,
        source_population::Real = ITURI_POPULATION,
        infection = infection_model,
        onset_incidence = onset_incidence_model,
        exports = exports_model,
        ascertainment = pooled_ascertainment_model())
    latent ~ to_submodel(
        _latent(n, breakpoint, infection, onset_incidence), false)
    asc_state ~ to_submodel(ascertainment)
    exports_state ~ to_submodel(
        exports(exported_cases, latent.infection_state.infections,
        asc_state.p_uganda; export_case_days,
        incubation_pmf = latent.incubation_pmf,
        source_population))
end

Composer: deaths-only fit

@model function deaths_only_model(
        n::Integer, total_deaths::Union{Missing, Integer};
        deaths_history = (; days = Int[], counts = Int[]),
        suspected_daily_deaths_history = (; days = Int[], counts = Int[]),
        breakpoint::Union{Missing, Real} = missing,
        infection = infection_model,
        onset_incidence = onset_incidence_model,
        deaths = deaths_model,
        dispersion = surveillance_dispersion_model())
    latent ~ to_submodel(
        _latent(n, breakpoint, infection, onset_incidence), false)
    dispersion_state ~ to_submodel(dispersion)
    deaths_state ~ to_submodel(
        deaths(deaths_history, total_deaths, latent.onsets,
        dispersion_state.k; suspected_daily_deaths_history))
end

Composer: cases-only fit

@model function cases_only_model(
        n::Integer, reported_cases::Union{Missing, Integer};
        reported_history = (; days = Int[], counts = Int[]),
        suspected_daily_history = (; days = Int[], counts = Int[]),
        breakpoint::Union{Missing, Real} = missing,
        infection = infection_model,
        onset_incidence = onset_incidence_model,
        cases = reported_cases_model,
        dispersion = surveillance_dispersion_model(),
        ascertainment = pooled_ascertainment_model())
    latent ~ to_submodel(
        _latent(n, breakpoint, infection, onset_incidence), false)
    dispersion_state ~ to_submodel(dispersion)
    asc_state ~ to_submodel(ascertainment)
    cases_state ~ to_submodel(
        cases(reported_history, reported_cases, latent.onsets,
        dispersion_state.k, asc_state.p_drc; suspected_daily_history))
end

Composer: confirmed-only fit

@model function confirmed_only_model(
        n::Integer, confirmed_cases::Union{Missing, Integer};
        confirmed_history = (; days = Int[], counts = Int[]),
        lab_history = (; days = Int[], counts = Int[]),
        lab_daily_history = (; days = Int[], counts = Int[]),
        tests_analysed::Union{Missing, Integer} = missing,
        breakpoint::Union{Missing, Real} = missing,
        infection = infection_model,
        onset_incidence = onset_incidence_model,
        cases = reported_cases_model,
        confirmed = confirmed_cases_model,
        dispersion = surveillance_dispersion_model(),
        ascertainment = pooled_ascertainment_model(),
        confirmed_positivity_link::Symbol = :composition)
    latent ~ to_submodel(
        _latent(n, breakpoint, infection, onset_incidence), false)
    dispersion_state ~ to_submodel(dispersion)
    asc_state ~ to_submodel(ascertainment)
    k = dispersion_state.k
    p_drc = asc_state.p_drc
    cases_state ~ to_submodel(
        cases((; days = Int[], counts = Int[]), missing, latent.onsets,
        k, p_drc))
    confirmed_state ~ to_submodel(
        confirmed(confirmed_history, confirmed_cases, latent.onsets, k,
        p_drc, cases_state.bg_daily, cases_state.τ_test,
        cases_state.bvd_reports_daily;
        lab_history, lab_daily_history,
        tests_analysed,
        positivity_link = confirmed_positivity_link))
end

Composer: exports-deaths-only fit

@model function exports_deaths_only_model(
        n::Integer, exports_deaths::Union{Missing, Integer};
        export_death_days::AbstractVector{<:Integer} = Int[],
        breakpoint::Union{Missing, Real} = missing,
        source_population::Real = ITURI_POPULATION,
        infection = infection_model,
        onset_incidence = onset_incidence_model,
        deaths = deaths_model,
        exports = exports_model,
        dispersion = surveillance_dispersion_model(),
        ascertainment = pooled_ascertainment_model())
    latent ~ to_submodel(
        _latent(n, breakpoint, infection, onset_incidence), false)
    dispersion_state ~ to_submodel(dispersion)
    asc_state ~ to_submodel(ascertainment)
    deaths_state ~ to_submodel(
        deaths((; days = Int[], counts = Int[]), missing, latent.onsets,
        dispersion_state.k))
    exports_state ~ to_submodel(
        exports(missing, latent.infection_state.infections,
        asc_state.p_uganda; incubation_pmf = latent.incubation_pmf,
        source_population))
    exports_deaths_state ~ to_submodel(
        exports_deaths_model(exports_deaths,
        exports_state.travelled_prevalence, deaths_state.CFR,
        deaths_state.od_pmf, latent.incubation_pmf; export_death_days))
end

Composer: joint fit

@model function bvd_joint(
        n::Integer,
        exported_cases::Union{Missing, Integer},
        total_deaths::Union{Missing, Integer},
        reported_cases::Union{Missing, Integer} = missing,
        exports_deaths::Union{Missing, Integer} = missing,
        confirmed_cases::Union{Missing, Integer} = missing,
        tests_analysed::Union{Missing, Integer} = missing;
        confirmed_deaths::Union{Missing, Integer} = missing,
        recovered_cases::Union{Missing, Integer} = missing,
        deaths_history = (; days = Int[], counts = Int[]),
        reported_history = (; days = Int[], counts = Int[]),
        confirmed_history = (; days = Int[], counts = Int[]),
        confirmed_deaths_history = (; days = Int[], counts = Int[]),
        lab_history = (; days = Int[], counts = Int[]),
        lab_daily_history = (; days = Int[], counts = Int[]),
        suspected_daily_history = (; days = Int[], counts = Int[]),
        suspected_daily_deaths_history = (; days = Int[], counts = Int[]),
        isolation_history = (; days = Int[], counts = Int[]),
        bed_capacity_history = (; days = Int[], counts = Int[]),
        recovered_history = (; days = Int[], counts = Int[]),
        export_case_days::AbstractVector{<:Integer} = Int[],
        export_death_days::AbstractVector{<:Integer} = Int[],
        breakpoint::Union{Missing, Real} = missing,
        source_population::Real = ITURI_POPULATION,
        infection = infection_model,
        onset_incidence = onset_incidence_model,
        exports = exports_model,
        deaths = deaths_model,
        cases = reported_cases_model,
        confirmed = confirmed_cases_model,
        confirmed_deaths_stream = confirmed_deaths_model,
        treatment = treatment_admission_model,
        recovered = recovered_model,
        dispersion = pooled_dispersion_model,
        ascertainment = pooled_ascertainment_model(),
        background_re::Bool = false,
        confirmed_positivity_link::Symbol = :composition,
        genetic = nothing,
        tmrca_days::Union{Missing, Real} = missing,
        tmrca_days_sd::Real = 15.0,
        renewal_start_lead::Integer = RENEWAL_START_LEAD,
        rt_walk_lead::Integer = RT_WALK_LEAD)
    ## The renewal start sits `renewal_start_lead` days AFTER the genetic
    ## TMRCA day (`n - tmrca_days + lead`), past the TMRCA's uncertainty where
    ## sustained transmission is confident. The lead keeps the observed span
    ## `τ_obs = n − renewal_start` strictly shorter than `tmrca_days`, so the
    ## genetic bound on the total age `T = m·τ + τ_obs` stays informative (it
    ## bounds the cryptic duration `m·τ` from below). The renewal seeds and
    ## grows from here.
    rt_start = ismissing(tmrca_days) ? 1 :
               clamp(n - round(Int, tmrca_days) + renewal_start_lead, 1, n)
    ## Start the random walk `rt_walk_lead` days (a month by default) BEFORE
    ## the first situation report (`breakpoint`) rather than exactly at it, so
    ## R_t is free to move over the weeks of transmission leading up to the
    ## first report instead of being held flat at R0 right to it (the response
    ## decline can begin before the outbreak is first reported). The start is
    ## floored at the renewal start so the walk never precedes the seeded
    ## trajectory. With no breakpoint the walk falls back to the renewal start.
    rt_walk_start = ismissing(breakpoint) ? rt_start :
                    clamp(round(Int, breakpoint) - rt_walk_lead, rt_start, n)
    latent ~ to_submodel(
        _latent(n, breakpoint, infection, onset_incidence;
            rt_start, rt_walk_start), false)
    infection_state = latent.infection_state
    onsets = latent.onsets

    ## Partially-pooled per-stream dispersions: every count stream draws its
    ## own negative-binomial dispersion from a shared population rather than
    ## sharing one global `k`, so a stream's noise is not pulled around by
    ## whichever stream dominates the likelihood while the sparse streams
    ## still borrow strength. Order: 1 suspected cases, 2 suspected deaths,
    ## 3 confirmed cases, 4 confirmed deaths, 5 isolation occupancy,
    ## 6 recovered. The isolation and recovered dispersions are injected into
    ## their submodels (which sample their own only when run standalone).
    dispersion_state ~ to_submodel(dispersion(6))
    asc_state ~ to_submodel(ascertainment)
    kv = dispersion_state.k
    k_cases = kv[1]
    k_deaths = kv[2]
    k_confirmed = kv[3]
    k_confirmed_deaths = kv[4]
    k_isolation = kv[5]
    k_recovered = kv[6]
    p_drc = asc_state.p_drc
    p_uganda = asc_state.p_uganda

    ## Non-BVD background as a SMOOTH daily lognormal random walk over the
    ## surveillance window ([`background_walk_model`](@ref)), with the tight
    ## innovation SD `σ_rw` driving the suspected-CASE stream. The background is
    ## gated to zero before the surveillance onset (a report-to-receipt lead
    ## before the first suspected-case report) — it does not exist before
    ## surveillance began. The tight innovation SD keeps it fairly constant,
    ## which regularises the background/outbreak-size degeneracy (the prior used
    ## a per-vintage STEP random effect whose multiplicative blow-up opened a
    ## second posterior mode that broke convergence). The suspected-DEATH
    ## background is NOT a separate random effect: it is tied to the case
    ## background by a background CFR (`cfr_bg · case_bg_daily`, see
    ## [`deaths_model`](@ref)), so it inherits the case background's smooth,
    ## gated, ramped level and time-variation rather than competing as a second
    ## free, outbreak-size-degenerate rate. With `background_re = false` (the
    ## renewal default) the case stream keeps its scalar `λ_bg`.
    if background_re
        bg_pool ~ to_submodel(background_pooling_model())
        σ_rw_shared = bg_pool.σ_bg
        ## Onset of the suspected pool's non-BVD background: a report-to-receipt
        ## lead BEFORE the first suspected-case report, not exactly at it. The
        ## suspects in the first report were already in the pipeline, and the
        ## background feeds the laboratory analysed volume through the report-to-
        ## receipt convolution, so it must begin early enough for that
        ## convolution to be fully formed by the first report. The lead is the
        ## MAX lag of the report-to-receipt kernel (its truncation `nmax`, the
        ## default `lab_delay_model` support), not its mean, so no tail
        ## contribution is cut off at the onset.
        bg_lead = cdf_nmax(lognormal_meansd(4.5, 4.0))
        bg_onset = isempty(reported_history.days) ? 1 :
                   clamp(Int(reported_history.days[1]) - bg_lead, 1, n)
        case_bg_re = nn -> background_walk_model(nn, σ_rw_shared;
            onset = bg_onset)
    else
        case_bg_re = nothing
    end

    ## Cases first so the suspected-case background `bg_daily` is available to
    ## the deaths stream (which scales it by `cfr_bg` for the death background)
    ## and to the laboratory pipeline.
    cases_state ~ to_submodel(
        cases(reported_history, reported_cases, onsets, k_cases, p_drc;
        suspected_daily_history, background_re = case_bg_re))
    deaths_state ~ to_submodel(
        deaths(deaths_history, total_deaths, onsets, k_deaths;
        suspected_daily_deaths_history, case_bg_daily = cases_state.bg_daily))
    confirmed_state ~ to_submodel(
        confirmed(confirmed_history, confirmed_cases, onsets, k_confirmed,
        p_drc, cases_state.bg_daily, cases_state.τ_test,
        cases_state.bvd_reports_daily;
        lab_history, lab_daily_history,
        tests_analysed,
        positivity_link = confirmed_positivity_link))
    ## Confirmed deaths mirror the confirmed-case lab pipeline: the death
    ## analysed volume scales the modelled case analysed volume
    ## (`confirmed_state.analysed_daily`) at the per-day suspected
    ## death-to-case ratio, scored through a death-pool composition positivity
    ## from the death series' own BVD and background components. The case
    ## volume carries the laboratory capacity onset, so the death volume
    ## inherits it and no deaths are confirmed before testing began.
    confirmed_deaths_state ~ to_submodel(
        confirmed_deaths_stream(confirmed_deaths, total_deaths,
        deaths_state.deaths_daily, deaths_state.bvd_deaths_daily,
        deaths_state.bg_death_daily, k_confirmed_deaths;
        confirmed_deaths_history, receipt_pmf = confirmed_state.receipt_pmf,
        case_analysed_daily = confirmed_state.analysed_daily,
        case_suspected_daily = cases_state.reports_daily))
    ## Isolation/treatment-bed occupancy: the suspect inflow carried through a
    ## length-of-stay survival into a latent bed demand, soft-capped at the bed
    ## capacity the implied-capacity series pins (see
    ## [`treatment_admission_model`](@ref)). The non-BVD rule-out stay is a
    ## separate parameter from the lab-turnaround `receipt_pmf`.
    treatment_state ~ to_submodel(
        treatment(isolation_history, cases_state.bvd_reports_daily,
        cases_state.bg_daily, p_drc;
        capacity_history = bed_capacity_history, k_external = k_isolation))
    ## Recovered among confirmed ("cumul guéris"): survivors among the modelled
    ## daily confirmed cases, with a recovery fraction grounded on the CFR and
    ## lagged by a confirmation-to-recovery delay (see [`recovered_model`](@ref)).
    recovered_state ~ to_submodel(
        recovered(recovered_history, recovered_cases,
        confirmed_state.confirmed_daily, deaths_state.CFR;
        k_external = k_recovered))
    exports_state ~ to_submodel(
        exports(exported_cases, infection_state.infections, p_uganda;
        export_case_days, incubation_pmf = latent.incubation_pmf,
        source_population))
    exports_deaths_state ~ to_submodel(
        exports_deaths_model(exports_deaths,
        exports_state.travelled_prevalence, deaths_state.CFR,
        deaths_state.od_pmf, latent.incubation_pmf; export_death_days))

    if genetic !== nothing
        genetic_state ~ to_submodel(
            genetic(infection_state.T, tmrca_days; tmrca_days_sd), false)
    end

    ## Daily cumulative trajectories for the headline 3x2 figure: the
    ## modelled expected cumulative infections, symptom onsets and deaths
    ## over the grid. Exposed as vector deterministics so the ribbon panels
    ## reconstruct from the chain without re-running the renewal. All three
    ## are BVD-only latent renewal quantities: deaths uses the BVD death
    ## series (onsets convolved with the onset-to-death delay), NOT the
    ## fitted total, so it stays smooth like infections and onsets. The
    ## additive non-BVD background is a daily random walk and belongs to the
    ## observation side, not this latent trajectory. `cumulative_infections`
    ## and `C_T` are exposed once by the shared `_latent` submodel above.
    cumulative_onsets := cumsum(onsets)
    cumulative_expected_deaths := cumsum(deaths_state.bvd_deaths_daily)
    ## Modelled daily laboratory-confirmed cases (from `confirmed_cases_model`:
    ## the per-window tested-positive probability applied to the modelled,
    ## testing-onset-gated analysed volume), so the cumulative trajectory carries
    ## the confirmed-case timing for the delay-corrected confirmed-CFR
    ## reconstruction. The onset-to-confirmation kernel (onset-to-report ⊕
    ## receipt) and the onset-to-death-confirmation kernel (onset-to-death ⊕
    ## receipt) are exposed alongside so the residual delay between a confirmed
    ## case and its confirmed death can be rebuilt per draw off the chain.
    ## Re-add the testing-onset baseline: the laboratory capacity is gated to
    ## zero before testing began and the first confirmed vintage is treated as
    ## the initial condition (a baseline the early windows do not score), so the
    ## reconstructed cumulative counts only the fitted increments. Adding the
    ## first observed confirmed count back from the testing onset onward makes
    ## the trajectory comparable to the observed confirmed total (and keeps the
    ## delay-corrected confirmed-CFR denominator on the right level).
    _conf_inc_cum = cumsum(confirmed_state.confirmed_daily)
    _conf_base = isempty(confirmed_history.counts) ? 0 :
                 Int(confirmed_history.counts[1])
    _conf_cap = isempty(confirmed_history.days) ? 1 :
                clamp(Int(confirmed_history.days[1]), 1, n)
    _conf_base_vec = [t >= _conf_cap ? _conf_base : 0 for t in 1:n]
    cumulative_confirmed := _conf_inc_cum .+ _conf_base_vec
    onset_to_confirmation_pmf := convolve_pmf(cases_state.report_pmf, confirmed_state.receipt_pmf)
    onset_to_death_confirmation_pmf := convolve_pmf(deaths_state.od_pmf, confirmed_state.receipt_pmf)
    R0 := infection_state.R0
    r := infection_state.r
    r0 := infection_state.r0
    doubling_time := infection_state.doubling_time
    T := infection_state.T
    R_T := infection_state.Rt[n]
    expected_infections_T := infection_state.infections[n]
    CFR := deaths_state.CFR
    ## Population-level dispersion (`k`, the headline scalar) plus the
    ## partially-pooled per-stream dispersions and the pooling SD.
    k := dispersion_state.k_pop
    k_cases := kv[1]
    k_deaths := kv[2]
    k_confirmed := kv[3]
    k_confirmed_deaths := kv[4]
    dispersion_sd := dispersion_state.τ
    p_drc := asc_state.p_drc
    p_uganda := asc_state.p_uganda
    expected_deaths_T := deaths_state.expected_deaths_T
    expected_reports_T := cases_state.expected_reports
    expected_confirmed_T := confirmed_state.expected_confirmed
    expected_analysed_T := confirmed_state.expected_analysed
    _ecd = confirmed_deaths_state.expected_confirmed_deaths
    expected_confirmed_deaths_T := _ecd
    expected_exports_T := exports_state.expected_exports
    expected_exports_deaths_T := exports_deaths_state.expected_exports_deaths_T
    expected_isolation_T := treatment_state.expected_isolation
    expected_bed_demand_T := treatment_state.expected_bed_demand
    bed_shortfall_T := safe_rate(treatment_state.expected_bed_demand -
                                 treatment_state.expected_isolation)
    bed_capacity := treatment_state.capacity
    isolation_admission := treatment_state.p_iso
    isolation_bvd_admission := treatment_state.p_iso_bvd
    isolation_severity := treatment_state.δ_iso
    isolation_bvd_los_mean := treatment_state.bvd_los_mean
    isolation_ruleout_los_mean := treatment_state.ruleout_los_mean
    isolation_dispersion := treatment_state.k_isolation
    expected_recovered_T := recovered_state.expected_recovered
    recovery_probability := recovered_state.p_recover
    recovery_delay_mean := recovered_state.recovery_delay_mean
    recovered_dispersion := recovered_state.k_recovered
    tau_test := cases_state.τ_test
    lambda_bg := cases_state.λ_bg
    bg_sigma := cases_state.bg_sigma
    background_total := cases_state.bg_total
    death_ascertainment := deaths_state.p_death
    background_cfr := deaths_state.cfr_bg
    lambda_bg_death := deaths_state.λ_bg_death
    bg_death_sigma := deaths_state.bg_death_sigma
    background_death_total := deaths_state.bg_death_total
    tau_death := confirmed_deaths_state.τ_death
    death_testing_scaling := confirmed_deaths_state.scaling
    suspected_positivity := cases_state.positivity
    test_positivity := confirmed_state.p_positive
    death_composition := confirmed_deaths_state.q_death
    death_confirmation := confirmed_deaths_state.p_death_conf
end

Model fitting and evaluation

Prior predictive check

Before any observation is taken into account, what does the prior imply about replicated exports, deaths and reported cases? Draws from the prior over the unobserved data should bracket the observed counts.

Sample the joint prior

prior_chn = let
    breakpoint = obs.n - obs.who_first_sitrep_days
    m = bvd_joint(obs.n, missing, missing, missing, missing, missing;
        deaths_history = (; days = Int[], counts = Int[]),
        reported_history = (; days = Int[], counts = Int[]),
        confirmed_history = (; days = Int[], counts = Int[]),
        export_case_days = obs.export_case_days,
        export_death_days = obs.export_death_days,
        breakpoint = breakpoint,
        background_re = true,
        confirmed_positivity_link = :composition,
        genetic = genetic_seeding_model,
        tmrca_days = obs.tmrca_days)
    sample(m, Prior(), 2_000; progress = false)
end;

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	276.0	725.0	1494.0	6152.0	14838.0	80128.0

Pair plot of the prior over the latent quantities.

Prior pair plot

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

Fitting the models

We sample with NUTS (Hoffman and Gelman, 2014) and Mooncake (Tebbutt and Ge, 2024) reverse-mode automatic differentiation, running two chains of 1000 post-warmup draws each after 1000 warmup adaptation steps, at a target acceptance probability of 0.85. Chains initialise from the prior. We fit the joint model and each single-stream model so the per-stream posteriors over the outbreak size can be compared with the joint.

Run the joint and per-stream NUTS fits

const _BREAKPOINT = obs.n - obs.who_first_sitrep_days

# The joint and six single-stream fits are independent, so they run through
# `fit_parallel`: model-level parallelism bounded by the available threads
# (sequential on CI's two threads, several at once on a many-core machine).
# The exports-deaths composer keeps the deaths and exports submodels only
# for their CFR, onset-to-death PMF and export onsets, leaving their own
# counts missing, which leaves two redundant sampled discrete draws, so its
# model check is disabled (see `nuts_sample`).
# Setup for the single master fit pool: relocated cut-offs, frozen-fit and
# sensitivity-variant helpers so every independent fit can run in one pool.
validation_cutoff = string(obs.cutoff - Day(7))

# Published per-release estimates, pulled from the tagged results
# releases by `scripts/refresh_releases.jl` into
# `data/released_estimates.csv`. Columns: tag, date, model (integral or
# renewal), median and the 30/60/90% bounds.
released_df = CSV.read(
    joinpath(pkgdir(BVDOutbreakSize), "data", "released_estimates.csv"),
    DataFrame)

# Integral-era release cut-offs for the frozen renewal overlay: a
# release whose data cut-off already has a renewal release needs no
# re-fit, since that release already is the renewal estimate. 20, 23 and
# 27 May are additionally fit for the matched-cutoff comparison further
# down (27 May matches the Lancet publication's cut-off).
renewal_release_dates = Set(string(r.date)
for r in eachrow(released_df) if r.model == "renewal")
frozen_evolution_cutoffs = sort(unique(string(r.date)
for r in eachrow(released_df)
if r.model == "integral" && string(r.date) ∉ renewal_release_dates))
frozen_cutoffs = sort(union(frozen_evolution_cutoffs,
    ["2026-05-20", "2026-05-23", "2026-05-27"]))

# A joint fit at the full headline settings (1000 draws × 2 chains) to the
# data frozen at `cutoff_date`. The frozen named tuple has the same shape as
# the full `obs`, so the model call mirrors the headline joint fit.
function fit_frozen_joint(cutoff_date; samples = 1000, chains = 2)
    o = freeze_observations(cutoff_date)
    bp = o.n - o.who_first_sitrep_days
    chn = nuts_sample(
        bvd_joint(
            o.n, o.exported_cases, o.total_deaths,
            o.reported_cases, o.exports_deaths, o.confirmed_cases,
            o.tests_analysed;
            confirmed_deaths = o.confirmed_deaths,
            deaths_history = o.deaths_history,
            reported_history = o.reported_history,
            confirmed_history = o.confirmed_history,
            confirmed_deaths_history = o.confirmed_deaths_history,
            lab_history = o.lab_history,
            lab_daily_history = o.lab_daily_history,
            isolation_history = o.isolation_history,
            bed_capacity_history = o.bed_capacity_history,
            export_case_days = o.export_case_days,
            export_death_days = o.export_death_days,
            breakpoint = bp,
            background_re = true,
            confirmed_positivity_link = :composition,
            genetic = genetic_seeding_model,
            tmrca_days = o.tmrca_days);
        samples = samples, chains = chains,
        callback = fit_callback("frozen_$(cutoff_date)"))
    return (; cutoff = o.cutoff, o, chn)
end

# One joint re-fit on the live data at the full headline settings, with hooks
# to override the genetic-seeding bound and the deaths submodel. The deaths
# submodel is passed the same way the genetic-seeding override is, as a
# closure matching the joint's `deaths(history, total, onsets, k;
# background_re)` call, so an alternative onset-to-death delay can be injected
# without touching the package.
function refit_joint_variant(;
        deaths = deaths_model,
        tmrca_days = obs.tmrca_days,
        tmrca_days_sd = 15.0,
        samples = 1000, chains = 2)
    chn = nuts_sample(
        bvd_joint(
            obs.n, obs.exported_cases, obs.total_deaths,
            obs.reported_cases, obs.exports_deaths, obs.confirmed_cases,
            obs.tests_analysed;
            confirmed_deaths = obs.confirmed_deaths,
            recovered_cases = obs.recovered_cases,
            deaths_history = obs.deaths_history,
            reported_history = obs.reported_history,
            confirmed_history = obs.confirmed_history,
            confirmed_deaths_history = obs.confirmed_deaths_history,
            lab_history = obs.lab_history,
            lab_daily_history = obs.lab_daily_history,
            suspected_daily_history = obs.suspected_daily_history,
            suspected_daily_deaths_history =
            obs.suspected_daily_deaths_history,
            isolation_history = obs.isolation_history,
            bed_capacity_history = obs.bed_capacity_history,
            recovered_history = obs.recovered_history,
            export_case_days = obs.export_case_days,
            export_death_days = obs.export_death_days,
            breakpoint = _BREAKPOINT,
            background_re = true,
            confirmed_positivity_link = :composition,
            deaths = deaths,
            genetic = genetic_seeding_model,
            tmrca_days = tmrca_days,
            tmrca_days_sd = tmrca_days_sd);
        samples = samples, chains = chains,
        callback = fit_callback("variant"))
    return chn
end

# Community-pathway onset-to-death delay from the Isiro 2012 line-list
# reanalysis community-death model (a single Gamma, implied mean ≈ 8 d,
# shape ≈ 5.5), a closure that re-injects this delay on its natural Gamma
# shape and scale into the deaths submodel while keeping its other defaults.
deaths_community_delay = (history,
    total,
    onsets,
    k;
    kwargs...) -> deaths_model(history, total, onsets, k;
    onset_to_death = gamma_delay_model(40;
        alpha_prior = truncated(Normal(5.48, 2.0); lower = 0.01),
        theta_prior = truncated(Normal(1.49, 0.5); lower = 0.1)),
    kwargs...)

# The faster early-epidemic clock dates the common ancestor about 17 days
# more recently, so the bound on the outbreak age sits that many days
# closer to the cut-off, with a tighter spread from its narrower interval.
clock_alt_offset = value(Date("2026-04-11") - Date("2026-03-25"))
tmrca_days_alt = obs.tmrca_days - clock_alt_offset

# Sensitivity refits (onset-to-death delay, molecular clock) are slow extra
# joint fits. They are PAUSED due to compute constraints: the serial
# sensitivity re-fits pushed the main docs build toward the 6h CI cap, so they
# are not run on any build for now. The code below is retained unchanged;
# re-enable by restoring the `BVD_RUN_SENSITIVITY` env gate:
#   RUN_SENSITIVITY = lowercase(strip(get(ENV, "BVD_RUN_SENSITIVITY",
#       "false"))) in ("true", "1", "yes", "on")
RUN_SENSITIVITY = false

# Every independent fit runs as one work-stealing pool so the long joint fit
# overlaps the per-stream, frozen and (gated) sensitivity re-fits and keeps
# all cores busy, rather than the short fits idling cores while the joint
# finishes. The fits are data-only independent, so the schedule is free.
_headline_thunks = [
    () -> nuts_sample(
        bvd_joint(
            obs.n, obs.exported_cases, obs.total_deaths,
            obs.reported_cases, obs.exports_deaths, obs.confirmed_cases,
            obs.tests_analysed;
            confirmed_deaths = obs.confirmed_deaths,
            recovered_cases = obs.recovered_cases,
            deaths_history = obs.deaths_history,
            reported_history = obs.reported_history,
            confirmed_history = obs.confirmed_history,
            confirmed_deaths_history = obs.confirmed_deaths_history,
            lab_history = obs.lab_history,
            lab_daily_history = obs.lab_daily_history,
            suspected_daily_history = obs.suspected_daily_history,
            suspected_daily_deaths_history = obs.suspected_daily_deaths_history,
            isolation_history = obs.isolation_history,
            bed_capacity_history = obs.bed_capacity_history,
            recovered_history = obs.recovered_history,
            export_case_days = obs.export_case_days,
            export_death_days = obs.export_death_days,
            breakpoint = _BREAKPOINT,
            background_re = true,
            confirmed_positivity_link = :composition,
            genetic = genetic_seeding_model,
            tmrca_days = obs.tmrca_days);
        callback = fit_callback("joint")),
    # Exports fit cases and deaths jointly as one "exports" stream, sharing
    # the travel-gated at-risk prevalence, rather than as two separate fits.
    () -> nuts_sample(
        exports_joint_only_model(obs.n, obs.exported_cases,
            obs.exports_deaths;
            export_case_days = obs.export_case_days,
            export_death_days = obs.export_death_days,
            breakpoint = _BREAKPOINT);
        check_model = false, callback = fit_callback("exports")),
    () -> nuts_sample(
        deaths_only_model(obs.n, obs.total_deaths;
            deaths_history = obs.deaths_history,
            suspected_daily_deaths_history = obs.suspected_daily_deaths_history,
            breakpoint = _BREAKPOINT);
        callback = fit_callback("deaths")),
    () -> nuts_sample(
        cases_only_model(obs.n, obs.reported_cases;
            reported_history = obs.reported_history,
            suspected_daily_history = obs.suspected_daily_history,
            breakpoint = _BREAKPOINT);
        callback = fit_callback("cases")),
    () -> nuts_sample(
        confirmed_only_model(obs.n, obs.confirmed_cases;
            confirmed_history = obs.confirmed_history,
            lab_history = obs.lab_history,
            lab_daily_history = obs.lab_daily_history,
            breakpoint = _BREAKPOINT);
        callback = fit_callback("confirmed")),
    () -> nuts_sample(
        confirmed_deaths_only_model(obs.n, obs.confirmed_deaths,
            obs.total_deaths;
            deaths_history = obs.deaths_history,
            confirmed_deaths_history = obs.confirmed_deaths_history,
            breakpoint = _BREAKPOINT);
        callback = fit_callback("confirmed_deaths")),
    () -> nuts_sample(
        treatment_only_model(obs.n;
            isolation_history = obs.isolation_history,
            bed_capacity_history = obs.bed_capacity_history,
            breakpoint = _BREAKPOINT);
        callback = fit_callback("treatment"))
]
_frozen_thunks = [() -> fit_frozen_joint(c) for c in frozen_cutoffs]
_sens_thunks = RUN_SENSITIVITY ?
               [() -> refit_joint_variant(deaths = deaths_community_delay),
    () -> refit_joint_variant(
        tmrca_days = tmrca_days_alt, tmrca_days_sd = 9.0)] : []
_all_fits = fit_parallel(vcat(_headline_thunks,
    [() -> fit_frozen_joint(validation_cutoff)], _frozen_thunks, _sens_thunks))

_n_head = length(_headline_thunks)
_n_frozen = length(frozen_cutoffs)
(chn_joint, chn_exports, chn_deaths, chn_cases, chn_confirmed,
    chn_confirmed_deaths, chn_treatment) = _all_fits[1:_n_head]
frozen_lastweek = _all_fits[_n_head + 1]
frozen_results = _all_fits[(_n_head + 2):(_n_head + 1 + _n_frozen)]
frozen_by_cutoff = Dict(zip(frozen_cutoffs, frozen_results))
frozen_C(c) = vec(Array(frozen_by_cutoff[c].chn[:C_T]))
if RUN_SENSITIVITY
    chn_joint_community_delay = _all_fits[_n_head + 2 + _n_frozen]
    chn_joint_fast_clock = _all_fits[_n_head + 3 + _n_frozen]
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_confirmed_deaths = vec(Array(chn_confirmed_deaths[:C_T]));
posterior_C_treatment = vec(Array(chn_treatment[:C_T]));

# Clean display names for the summary tables and pair plots. The submodel
# prefixes (`rt_state.`, `gi_state.`, ...) are kept in the model so the
# nested submodels stay distinct; this map only relabels them for display.
display_names = Dict{Symbol, String}(
    Symbol("rt_state.sigma_rw") => "Rt step size",
    Symbol("rt_state.intervention_effect") => "intervention effect",
    Symbol("gi_state.α") => "generation interval shape",
    Symbol("gi_state.θ") => "generation interval scale",
    Symbol("inc_state.delay_mean") => "incubation period mean",
    Symbol("inc_state.delay_sd") => "incubation period SD",
    Symbol("cases_state.report_state.α") => "onset-to-report shape",
    Symbol("cases_state.report_state.θ") => "onset-to-report scale",
    Symbol("deaths_state.od_state.oa.α") => "onset-to-admission shape",
    Symbol("deaths_state.od_state.oa.θ") => "onset-to-admission scale",
    Symbol("deaths_state.od_state.ad.α") => "admission-to-death shape",
    Symbol("deaths_state.od_state.ad.θ") => "admission-to-death scale",
    Symbol("exports_state.detect_state.α") => "onset-to-detection shape",
    Symbol("exports_state.detect_state.θ") => "onset-to-detection scale",
    Symbol("confirmed_state.receipt_state.d.delay_mean") => "report-to-receipt mean",
    Symbol("confirmed_state.receipt_state.d.delay_sd") => "report-to-receipt SD",
    :isolation_bvd_los_mean => "isolation BVD length-of-stay mean",
    :isolation_ruleout_los_mean => "isolation non-BVD rule-out stay mean",
    :recovery_delay_mean => "confirmation-to-recovery mean",
    Symbol("exports_state.travel_state.daily_travellers") => "daily travellers");

┌ Info: Found initial step size
└   ϵ = 0.05
┌ Info: Found initial step size
└   ϵ = 0.4
┌ Warning: There were 2 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Warning: There were 2 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Info: Found initial step size
└   ϵ = 0.025
┌ Info: Found initial step size
└   ϵ = 0.0125
┌ Info: Found initial step size
└   ϵ = 0.05
┌ Info: Found initial step size
└   ϵ = 0.0125
┌ Warning: There were 10 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Warning: There were 10 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Info: Found initial step size
└   ϵ = 0.025
┌ Info: Found initial step size
└   ϵ = 0.00625
┌ Warning: There were 29 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Warning: There were 2 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Warning: There were 31 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Info: Found initial step size
└   ϵ = 0.0125
┌ Info: Found initial step size
└   ϵ = 0.00625
┌ Warning: There were 28 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Warning: There were 20 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Warning: There were 48 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Info: Found initial step size
└   ϵ = 0.05
┌ Info: Found initial step size
└   ϵ = 0.0125
┌ Warning: There were 1 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Warning: There were 1 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Info: Found initial step size
└   ϵ = 0.025
┌ Info: Found initial step size
└   ϵ = 0.00625
┌ Warning: There were 174 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Warning: There were 78 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Warning: There were 252 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Info: Found initial step size
└   ϵ = 0.025
┌ Info: Found initial step size
└   ϵ = 0.0125
┌ Warning: There were 19 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Warning: There were 7 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Warning: There were 26 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Info: Found initial step size
└   ϵ = 0.05
┌ Info: Found initial step size
└   ϵ = 0.2
┌ Info: Found initial step size
└   ϵ = 0.05
┌ Info: Found initial step size
└   ϵ = 0.2
┌ Warning: There were 2 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Warning: There were 10 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Warning: There were 12 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Info: Found initial step size
└   ϵ = 0.2
┌ Info: Found initial step size
└   ϵ = 0.05
┌ Warning: There were 16 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Warning: There were 7 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Warning: There were 23 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Info: Found initial step size
└   ϵ = 0.05
┌ Info: Found initial step size
└   ϵ = 0.025
┌ Warning: There were 12 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Warning: There were 11 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483
┌ Warning: There were 23 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4hMHm/src/mcmc/hmc.jl:483

Fit diagnostics

Fit-quality diagnostics for the joint and per-stream fits: the worst R-hat, the smallest bulk effective sample size, and the number of divergent transitions.

Fit diagnostics

8×4 DataFrame

Row	fit	max_rhat	min_ess_bulk	divergences
	String	Float64	Float64	Int64
1	joint	1.011	503.0	23
2	exports	1.004	889.0	2
3	deaths (DRC)	1.021	120.0	10
4	cases (DRC)	1.006	717.0	31
5	confirmed (DRC)	1.009	433.0	48
6	confirmed deaths (DRC)	1.01	991.0	1
7	isolation (DRC)	1.018	425.0	252
8	frozen (1wk back)	1.006	845.0	26

No-onward-transmission counterfactual

To bound the deaths already committed at the cut-off, we project the deaths that would still occur if all transmission stopped on the report date. Every infection present by the cut-off still dies with probability CFR, so the committed future deaths are the CFR-weighted cumulative infection count net of the deaths already expected, $\mathrm{\Delta }D=\mathrm{CFR}\cdot I_T-\mathbb{E}[D_T]$, where $I_T$ is the cumulative infection count to the cut-off. The figure is shown in the counterfactual results below.

Delay-corrected confirmed case-fatality ratio

The case-fatality ratio above is the onset-level CFR, the share of symptomatic infections that die. It is hard to read directly off the data because the case and death streams are ascertained differently, so a reader who wants a figure anchored in the observed counts is left with the naive confirmed ratio, the cumulative confirmed deaths over the cumulative confirmed cases. That naive ratio is biased low in real time. A case confirmed close to the cut-off has not yet had time to die, so it enters the denominator before it can enter the numerator.

We report a delay-corrected confirmed CFR that debiases the real-time ratio following (Nishiura et al., 2009). The denominator is shrunk from all confirmed cases to those expected to have had their death confirmed by the cut-off, weighting each day of confirmed-case incidence by the probability that a case confirmed that day, if it is going to die, has had its death confirmed by the cut-off:

$$\begin{array}{}\text{(35)}& {\mathrm{cCFR}}_{\text{corr}}(T)=\frac{D_{\text{conf}}(T)}{\sum _t{c}_{\text{conf}}(t)Pr(X_d-X_c\le T-t)},\end{array}$$

with $D_{\text{conf}}(T)$ the cumulative confirmed deaths, $c_{\text{conf}}(t)$ the modelled daily confirmed-case incidence, and $X_d-X_c$ the residual delay between a confirmed case and its confirmed death. $X_d$ is the onset-to-death-confirmation lag (onset-to-death convolved with the report-to-receipt laboratory delay) and $X_c$ the onset-to-confirmation lag (onset-to-report convolved with the same laboratory delay), so the common receipt delay cancels in the mean and the residual centres on onset-to-death minus onset-to-report. Both lags and the confirmed trajectories are taken per posterior draw from the joint fit, so the corrected ratio carries the joint uncertainty. As the outbreak matures and recent incidence resolves the correction shrinks and the corrected ratio approaches the eventual confirmed CFR. It is the confirmed-case counterpart of the structural CFR, anchored in the confirmed counts rather than the latent infections, and the gap between the two reflects the difference in case and death ascertainment the structural CFR has to absorb. The result is shown in the confirmed case-fatality ratio results below.

One-week-ahead forecast

We project each DRC stream seven days beyond the cut-off, letting the reproduction number keep evolving over the horizon by continuing the recent trend of its trajectory rather than holding it fixed, with no further interventions and no saturation imposed. The projection carries both parameter and observation uncertainty. We forecast the two confirmed DRC streams (laboratory-confirmed cases and confirmed deaths) as the forecast targets, and also the isolation/treatment beds and the cumulative recovered total. For the beds we project both the bed demand (the need a week ahead, under unconstrained supply, the cut-off demand grown by the horizon factor like the case inflow) and the supply-limited occupancy that demand produces against the bed capacity. The gap between them is the projected bed shortfall, the quantity of interest if bed occupancy is supply-constrained. The suspected case and death streams are no longer published, so they are not shown as targets. Exports are not forecast either, since cross-border travel is unlikely to continue at its baseline rate, so the forward travel rate the export model relies on no longer holds. The figure is shown in the one-week-ahead forecast results below.

Forecast-versus-frozen evaluation

We assess the forecast against data observed since by freezing the data to roughly one week before the current cut-off, re-fitting, and projecting one week ahead with the same forecast machinery, then comparing that projection against the counts observed by the current cut-off. The frozen re-fit cuts the data to an earlier cut-off and re-fits the joint model, so that a change driven by newer data can be distinguished from one driven by a change of method. Each frozen re-fit uses the full headline settings (1000 draws across two chains). The same frozen re-fit is reused to compare against McCabe et al. at the cut-offs they used, and to trace how the estimate moved as the situation reports accrued, re-fitting the renewal model frozen at each release date. The helper below performs one frozen joint re-fit and is reused by the forecast validation, estimate-evolution and matched-in-time results.

Frozen-fit helper (reused by the forecast validation, evolution and matched-in-time sections)

# fit_frozen_joint and the frozen re-fits are defined and run in the setup
# block above.

Results

Summary

The numbers below are our estimate of the underlying infections to date, reported and unreported, from the joint posterior. Each is given as equal-tailed 30%, 60% and 90% credible intervals.

Compute the headline ranges

summary_ranges = let
    med(x) = quantile(x, 0.5)
    iqr(x) = quantile(x, 0.75) - quantile(x, 0.25)
    # Posterior-minus-prior shift in units of the parameter's prior IQR,
    # reusing the prior draws so nothing is respecified here.
    shift(post, prior) = round((med(post) - med(prior)) / iqr(prior);
        digits = 2)

    C = posterior_C_joint
    Td = vec(Array(chn_joint[:T]))
    r0d = vec(Array(chn_joint[:r0]))
    rd = vec(Array(chn_joint[:r]))
    dt0 = log(2) ./ r0d
    dt = vec(Array(chn_joint[:doubling_time]))
    R0d = vec(Array(chn_joint[:R0]))
    RTd = vec(Array(chn_joint[:R_T]))
    cfrd = vec(Array(chn_joint[:CFR]))
    sC = posterior_summary(C)
    sT = posterior_summary(Td)
    sr0 = posterior_summary(r0d)
    sr = posterior_summary(rd)
    sdt0 = posterior_summary(dt0)
    sdt = posterior_summary(dt)
    sR0 = posterior_summary(R0d)
    sRT = posterior_summary(RTd)
    scfr = posterior_summary(cfrd)

    ints_i(s) = string(
        "30% ", round(Int, s.lo30), "–", round(Int, s.hi30),
        ", 60% ", round(Int, s.lo60), "–", round(Int, s.hi60),
        ", 90% ", round(Int, s.lo90), "–", round(Int, s.hi90))
    ints_f(s,
        d) = string(
        "30% ", round(s.lo30; digits = d), "–", round(s.hi30; digits = d),
        ", 60% ", round(s.lo60; digits = d), "–", round(s.hi60; digits = d),
        ", 90% ", round(s.lo90; digits = d), "–", round(s.hi90; digits = d))
    start_from(t) = obs.cutoff - Day(round(Int, t))
    ints_d(s) = string(
        "30% ", start_from(s.hi30), "–", start_from(s.lo30),
        ", 60% ", start_from(s.hi60), "–", start_from(s.lo60),
        ", 90% ", start_from(s.hi90), "–", start_from(s.lo90))
    f_lo = round(sC.lo90 / obs.confirmed_cases; digits = 1)
    f_hi = round(sC.hi90 / obs.confirmed_cases; digits = 1)

    # How far the data has moved each estimate from its prior, in prior
    # interquartile ranges, reusing the prior draws.
    moves = [
        "cumulative infection count" => shift(C, vec(Array(prior_chn[:C_T]))),
        "outbreak age" => shift(Td, vec(Array(prior_chn[:T]))),
        "doubling time" => shift(dt, vec(Array(prior_chn[:doubling_time])))]
    biggest = argmax(p -> abs(p.second), moves)

    Markdown.parse("""
    - **Cumulative infections:** the outbreak is estimated to have caused
      $(ints_i(sC)) infections to date, reported and unreported.
    - Against the $(obs.confirmed_cases) laboratory-confirmed cases by the
      cut-off that is roughly $(f_lo)–$(f_hi)× as many infections, so
      confirmed cases are estimated to capture only a small share of the
      outbreak.
    - **Outbreak start and age:** the outbreak is estimated to have begun on
      a start date of $(ints_d(sT)), an elapsed age to the cut-off of
      $(ints_i(sT)) days.
    - **Growth rate and doubling time:** the initial growth rate is
      estimated to have been $(ints_f(sr0, 3)) per day, an initial doubling
      time of $(ints_f(sdt0, 1)) days.
      The latest growth rate is estimated to be $(ints_f(sr, 3)) per day, a
      latest doubling time of $(ints_f(sdt, 1)) days.
    - **Reproduction number:** the initial reproduction number is estimated
      to have been $(ints_f(sR0, 2)) and the latest to be $(ints_f(sRT, 2)).
    - **Case-fatality ratio:** the case-fatality ratio is estimated to be
      $(ints_f(scfr, 2)).
    - **Shift from priors:** how far the data has moved each estimate from
      its prior, in prior interquartile ranges, where a value of one means
      the posterior median sits one prior interquartile range from the prior
      median, zero means unchanged, and the sign gives the direction.
      The fit moves the cumulative infection count by $(moves[1].second),
      the outbreak age by $(moves[2].second) and the doubling time by
      $(moves[3].second); the largest move is in the $(biggest.first).
    """)
end;

• Cumulative infections: the outbreak is estimated to have caused 30% 2104–2623, 60% 1896–3024, 90% 1632–3838 infections to date, reported and unreported.

• Against the 1048 laboratory-confirmed cases by the cut-off that is roughly 1.6–3.7× as many infections, so confirmed cases are estimated to capture only a small share of the outbreak.

• Outbreak start and age: the outbreak is estimated to have begun on a start date of 30% 2026-02-24–2026-03-06, 60% 2026-02-17–2026-03-12, 90% 2026-02-05–2026-03-22, an elapsed age to the cut-off of 30% 107–117, 60% 101–124, 90% 91–136 days.

• Growth rate and doubling time: the initial growth rate is estimated to have been 30% 0.033–0.037, 60% 0.031–0.04, 90% 0.028–0.045 per day, an initial doubling time of 30% 18.8–21.0, 60% 17.5–22.4, 90% 15.4–25.2 days. The latest growth rate is estimated to be 30% -0.088–-0.042, 60% -0.129–-0.02, 90% -0.267–0.02 per day, a latest doubling time of 30% -12.2–-5.9, 60% -19.2–-3.4, 90% -64.3–41.7 days.

• Reproduction number: the initial reproduction number is estimated to have been 30% 1.15–1.25, 60% 1.12–1.33, 90% 1.07–1.51 and the latest to be 30% 0.6–0.78, 60% 0.51–0.89, 90% 0.36–1.12.

• Case-fatality ratio: the case-fatality ratio is estimated to be 30% 0.38–0.45, 60% 0.34–0.49, 90% 0.27–0.57.

• Shift from priors: how far the data has moved each estimate from its prior, in prior interquartile ranges, where a value of one means the posterior median sits one prior interquartile range from the prior median, zero means unchanged, and the sign gives the direction. The fit moves the cumulative infection count by -0.07, the outbreak age by -0.5 and the doubling time by -0.56; the largest move is in the doubling time.

Joint model estimates

This section reports the joint posterior over the cumulative infection count to date, fitting every data stream together.

Cumulative infection count summary table

cumulative_cases_summary = summary_table(
    chn_joint, [:C_T]; digits = 0);

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	1632.0	1896.0	2104.0	2623.0	3024.0	3838.0

The figure below shows three modelled cumulative quantities over time, one per row, all latent. The top row is infections, the middle row symptom onsets and the bottom row deaths. The left column is the modelled cumulative trajectory with its 50% and 90% intervals; the right column is the posterior density of the current cut-off cumulative. The infection density on the top right is the headline outbreak size, a count of infections rather than reported cases. No observed counts are overlaid, since each quantity sits upstream of the ascertainment, confirmation and reporting that produce the observations.

Cumulative infections, onsets and deaths figure

cumulative_traj_fig = plot_cumulative_trajectories(chn_joint;
    n = obs.n, seeding = obs.seeding);

The cumulative infection count is set by the reproduction number trajectory and the outbreak age, the elapsed time from the import that started the outbreak to the cut-off. Read as a calendar date, the age places the outbreak start at the cut-off date minus the age. The left panel below shows the posterior for that start date; the right panel shows the joint posterior of the outbreak age and the early doubling time.

Outbreak start date and seeding-time posterior

start_date_fig = plot_start_date_pair(chn_joint;
    as_of_date = string(obs.cutoff));

The summary table reports the credible intervals on the infection-process parameters: the growth rate and doubling time, the reproduction number, the outbreak age, the case-fatality ratio and the cumulative infection count. The pair plot beside it shows their joint distribution, with the prior overlaid so the data's contribution to each marginal is visible.

Infection-parameter summary table

infection_summary = summary_table(chn_joint,
    [:r, :doubling_time, :T, :R_T, :CFR, :C_T]; digits = 2);

6×7 DataFrame

Row	Quantity	Lower 90%	Lower 60%	Lower 30%	Upper 30%	Upper 60%	Upper 90%
	String	Float64	Float64	Float64	Float64	Float64	Float64
1	r	-0.27	-0.13	-0.09	-0.04	-0.02	0.02
2	doubling_time	-64.31	-19.22	-12.21	-5.94	-3.39	41.73
3	T	90.79	101.1	106.89	116.78	123.75	135.89
4	R_T	0.36	0.51	0.6	0.78	0.89	1.12
5	CFR	0.27	0.34	0.38	0.45	0.49	0.57
6	C_T	1631.93	1896.35	2104.22	2622.67	3023.56	3837.88

Infection-parameter pair plot (prior overlaid)

infection_pair_fig = plot_pair(chn_joint,
    [:R_T, :r, :T, :CFR,
        Symbol("rt_state.sigma_rw"), Symbol("rt_state.intervention_effect")];
    prior = prior_chn, labels = display_names);

The infection model carries two delays: the generation interval, the time between an infector's and an infectee's onset that drives the renewal recursion, and the incubation period, the time from infection to symptom onset that turns infections into onsets. The table reports their posterior means and standard deviations; the pair plot beside it shows their joint posterior with the prior overlaid.

Infection-delay summary table

infection_delay_summary = summary_table(chn_joint,
    [Symbol("gi_state.α"), Symbol("gi_state.θ"),
        Symbol("inc_state.delay_mean"), Symbol("inc_state.delay_sd")];
    digits = 2, labels = display_names);

Show infection-delay summary table

4×7 DataFrame

Row	Quantity	Lower 90%	Lower 60%	Lower 30%	Upper 30%	Upper 60%	Upper 90%
	String	Float64	Float64	Float64	Float64	Float64	Float64
1	generation interval shape	0.51	1.05	1.39	2.02	2.4	2.99
2	generation interval scale	0.72	1.98	2.92	4.52	5.31	6.74
3	incubation period mean	5.48	5.93	6.19	6.62	6.88	7.29
4	incubation period SD	2.06	2.7	3.06	3.63	3.99	4.63

Infection-delay pair plot (prior overlaid)

infection_delay_pair_fig = plot_pair(chn_joint,
    [Symbol("gi_state.α"), Symbol("gi_state.θ"),
        Symbol("inc_state.delay_mean"), Symbol("inc_state.delay_sd")];
    prior = prior_chn, labels = display_names);

Reproduction number over time

The daily reproduction number over the period we estimate it for, the established outbreak from the genetic bound to the cut-off. The 30%, 60% and 90% credible ribbons are shown with about a hundred sampled trajectories, and the no-growth threshold at one as a grey dashed line. The first situation report on 18 May 2026 marks the start of the response scale-up (red dashed) and the end of the three-week scale-up is the red dotted line; the data cut-off is grey dashed.

Reproduction-number trajectory

# `rt_start` is the renewal/established-window start the plot shows from;
# `rt_walk_start` is where the random walk's knots begin — `RT_WALK_LEAD`
# days (a month) before the first situation report, matching `bvd_joint`'s
# `rt_walk_lead` — so the chain reconstruction uses the same knot grid the
# model did, floored at the renewal start. R_t is flat at R0 between the two.
_rt_start_plot = clamp(
    obs.n - round(Int, obs.tmrca_days) + RENEWAL_START_LEAD, 1, obs.n);
rt_fig = plot_rt(chn_joint;
    n = obs.n, breakpoint = _BREAKPOINT,
    rt_start = _rt_start_plot,
    rt_walk_start = clamp(_BREAKPOINT - RT_WALK_LEAD, _rt_start_plot, obs.n),
    as_of_date = string(obs.cutoff), seeding = obs.seeding,
    ramp = 21.0);

The table reports the posterior of the response effect on the reproduction number as a multiplier, where a value below one is the factor by which the response lowers the reproduction number once the scale-up completes.

Intervention-effect summary table

intervention_effect = vec(Array(
    chn_joint[Symbol("rt_state.intervention_effect")]));
intervention_table = streams_table(
    "Rt multiplier exp(effect)" => exp.(intervention_effect);
    digits = 2);

1×7 DataFrame

Row	Stream	Lower 90%	Lower 60%	Lower 30%	Upper 30%	Upper 60%	Upper 90%
	String	Float64	Float64	Float64	Float64	Float64	Float64
1	Rt multiplier exp(effect)	0.49	0.62	0.7	0.84	0.9	0.97

Observation delays

The delays from symptom onset to each observed event: onset to a suspected case being reported, onset to death, onset to detection abroad (the export model uses the same onset-to-admission delay as the report), and the report-to-laboratory receipt delay. The onset-to-report and onset-to-detection delays are the same line-list onset-to-admission delay, sampled on its natural Gamma shape and scale, and onset-to-death is the convolution of two atomic Gamma delays, onset to admission and admission to death, each with its own shape and scale. The report-to-receipt delay is sampled by its mean and standard deviation. The length-of-stay delays are also shown: the isolation-bed BVD treatment length-of-stay — for how long an admitted BVD patient occupies a bed, with the line-list admission-to-death delay as its prior; the non-BVD rule-out stay — how long a ruled-out suspect occupies a bed before discharge, with the report-to-receipt turnaround as its prior; and the confirmation-to-recovery delay (how long after confirmation a case is recorded as recovered). The table reports their posteriors; the pair plot beside it shows their joint posterior with the prior overlaid, so the data's contribution to each marginal is visible.

Observation-delay summary table

obs_delay_summary = summary_table(chn_joint,
    [Symbol("cases_state.report_state.α"),
        Symbol("cases_state.report_state.θ"),
        Symbol("deaths_state.od_state.oa.α"),
        Symbol("deaths_state.od_state.oa.θ"),
        Symbol("deaths_state.od_state.ad.α"),
        Symbol("deaths_state.od_state.ad.θ"),
        Symbol("exports_state.detect_state.α"),
        Symbol("exports_state.detect_state.θ"),
        Symbol("confirmed_state.receipt_state.d.delay_mean"),
        Symbol("confirmed_state.receipt_state.d.delay_sd"),
        :isolation_bvd_los_mean,
        :isolation_ruleout_los_mean,
        :recovery_delay_mean];
    digits = 2, labels = display_names);

Show observation-delay summary table

13×7 DataFrame

Row	Quantity	Lower 90%	Lower 60%	Lower 30%	Upper 30%	Upper 60%	Upper 90%
	String	Float64	Float64	Float64	Float64	Float64	Float64
1	onset-to-report shape	0.7	0.93	1.06	1.29	1.42	1.65
2	onset-to-report scale	1.5	2.59	3.11	4.01	4.59	5.67
3	onset-to-admission shape	0.73	0.93	1.05	1.26	1.4	1.62
4	onset-to-admission scale	1.6	2.65	3.18	4.11	4.67	5.56
5	admission-to-death shape	1.16	1.58	1.83	2.24	2.51	2.93
6	admission-to-death scale	1.85	2.85	3.38	4.33	4.95	5.95
7	onset-to-detection shape	0.74	0.98	1.1	1.31	1.44	1.67
8	onset-to-detection scale	1.98	2.93	3.44	4.27	4.77	5.6
9	report-to-receipt mean	2.78	4.01	4.53	5.26	5.64	6.24
10	report-to-receipt SD	2.0	2.75	3.13	3.82	4.24	4.9
11	isolation BVD length-of-stay mean	4.95	7.64	9.08	11.88	13.67	17.31
12	isolation non-BVD rule-out stay mean	1.4	2.33	3.0	4.41	5.3	6.78
13	confirmation-to-recovery mean	14.11	17.59	19.27	21.76	23.12	25.8

Observation-delay pair plot (prior overlaid)

obs_delay_pair_fig = plot_pair(chn_joint,
    [Symbol("cases_state.report_state.α"),
        Symbol("deaths_state.od_state.oa.α"),
        Symbol("exports_state.detect_state.α"),
        Symbol("confirmed_state.receipt_state.d.delay_mean"),
        :isolation_bvd_los_mean,
        :isolation_ruleout_los_mean,
        :recovery_delay_mean];
    prior = prior_chn, labels = display_names);

Surveillance parameters

The surveillance-data parameters: the reporting fractions for the DRC and Uganda, the surveillance dispersions, and the laboratory pipeline (the testing fraction and receipt delay, the per-suspected and per-test positivity, the non-BVD background rate, and the death-confirmation probability). The six passive-surveillance count streams (suspected cases, suspected deaths, confirmed cases, confirmed deaths, isolation occupancy and recovered) each have their own negative-binomial dispersion partially pooled from a shared population: $k$ is the population-level dispersion, $k_{\text{cases}}$, $k_{\text{deaths}}$, $k_{\text{confirmed}}$ and $k_{\text{confirmed deaths}}$ the per-stream values for the four DRC count streams, and a pooling spread. The isolation and recovered streams add the proportion of suspects admitted to a bed and the recovery probability among confirmed cases, with their dispersions ($k_{\text{iso}}$, $k_{\text{rec}}$) drawn from the same pooled population (see the length-of-stay delays in the observation-delay table above). The table reports their credible intervals; the pair plot beside it shows their joint posterior with the prior overlaid.

Surveillance-parameter summary table

surveillance_summary = summary_table(chn_joint,
    [:p_drc, :p_uganda, :k, :k_cases, :k_deaths, :k_confirmed,
        :k_confirmed_deaths, :dispersion_sd, :tau_test, :lambda_bg,
        :suspected_positivity, :test_positivity, :expected_confirmed_T,
        :expected_analysed_T, :death_ascertainment, :background_cfr,
        :tau_death, :death_composition,
        :death_confirmation, :expected_confirmed_deaths_T,
        :isolation_admission, :isolation_dispersion, :expected_isolation_T,
        :expected_bed_demand_T, :bed_capacity, :bed_shortfall_T,
        :recovery_probability, :recovered_dispersion, :expected_recovered_T];
    digits = 3);

Show surveillance-parameter summary table

29×7 DataFrame

Row	Quantity	Lower 90%	Lower 60%	Lower 30%	Upper 30%	Upper 60%	Upper 90%
	String	Float64	Float64	Float64	Float64	Float64	Float64
1	p_drc	0.422	0.555	0.644	0.779	0.85	0.933
2	p_uganda	0.436	0.582	0.667	0.802	0.871	0.941
3	k	2.009	3.15	3.884	5.54	6.747	9.36
4	k_cases	2.683	3.499	3.975	5.2	6.118	7.796
5	k_deaths	4.594	6.47	7.974	11.777	14.748	21.685
6	k_confirmed	1.329	1.674	1.901	2.34	2.608	3.165
7	k_confirmed_deaths	0.683	0.873	0.997	1.272	1.479	1.909
8	dispersion_sd	0.531	0.625	0.685	0.788	0.859	0.992
9	tau_test	0.792	0.853	0.884	0.934	0.955	0.98
10	lambda_bg	11.453	14.876	17.399	21.624	24.265	28.831
11	suspected_positivity	0.224	0.241	0.253	0.277	0.292	0.331
12	test_positivity	0.236	0.242	0.246	0.253	0.256	0.264
13	expected_confirmed_T	949.663	974.997	991.171	1018.12	1033.77	1063.76
14	expected_analysed_T	3248.24	3490.13	3607.85	3851.3	4008.82	4276.02
15	death_ascertainment	0.805	0.858	0.886	0.918	0.933	0.954
16	background_cfr	0.153	0.206	0.239	0.297	0.339	0.403
17	tau_death	0.537	0.643	0.709	0.829	0.922	1.098
18	death_composition	0.312	0.396	0.445	0.544	0.602	0.714
19	death_confirmation	0.269	0.336	0.374	0.454	0.505	0.59
20	expected_confirmed_deaths_T	224.245	256.873	275.518	315.269	338.808	396.613
21	isolation_admission	0.162	0.21	0.247	0.32	0.368	0.473
22	isolation_dispersion	125.786	201.359	254.354	384.611	516.729	1011.58
23	expected_isolation_T	346.712	360.798	369.392	383.41	392.805	411.696
24	expected_bed_demand_T	346.712	360.798	369.392	383.41	392.805	411.696
25	bed_capacity	417.635	425.808	430.762	439.804	445.822	456.523
26	bed_shortfall_T	0.0	0.0	0.0	0.0	0.0	0.0
27	recovery_probability	0.305	0.388	0.446	0.558	0.624	0.727
28	recovered_dispersion	1.083	1.668	2.103	3.245	4.293	7.252
29	expected_recovered_T	99.715	113.932	123.879	144.439	160.121	194.376

Surveillance-parameter pair plot (prior overlaid)

surveillance_pair_fig = plot_pair(chn_joint,
    [:p_drc, :p_uganda, :k, :tau_test, :lambda_bg, :test_positivity,
        :death_confirmation];
    prior = prior_chn);

Posterior predictive checks

A posterior predictive check draws replicated observations from the fitted joint model and compares them to the observed counts. The checks read in three groups, following the generative order of the model. The first is the infection process, the latent infections, symptom onsets and deaths that drive every stream; these are not observed directly, so they carry no genuine replicate and are shown as the estimated cumulative trajectories in the joint model estimates figure rather than checked against data here. The second is the surveillance data, the dated DRC streams that are real per-vintage observations: cumulative suspected cases, the daily new-suspect inflow, the daily isolation-bed occupancy, confirmed cases, suspected deaths, confirmed deaths, recovered-among-confirmed and specimens analysed. The third is the exports, the cross-border imported cases and deaths detected in Uganda.

The surveillance group is checked first. Each panel is shown over its own reporting dates with the observed series overlaid: the cumulative streams as replicated cumulative trajectories, and the daily new-suspect inflow and the daily isolation-bed occupancy on a daily scale (each day's replicated count against the observed count). The cumulative suspected case and death streams stop at their last stable vintage on 26 May; the daily new-suspect inflow then runs 4-11 June, where the cumulative suspected series freezes, and the isolation occupancy runs 1-11 June; the laboratory-confirmed streams keep reporting to the cut-off.

Joint posterior predictive plot

# Drop the increment counts but keep each stream's vintage day grid, so
# `predict` resamples the per-vintage increments rather than holding them
# at the observed values. The confirmed-case windows and the per-window
# positivity random effect are defined by the confirmed and laboratory
# histories, so those are passed with their counts intact (only the
# cut-off scalars are set to `missing`) to keep the generator's latent
# dimensions identical to the fitted chain.
_days_only(h) = (; days = h.days, counts = Int[]);

pp_joint = predict(
    bvd_joint(
        obs.n, missing, missing, missing, missing, missing, missing;
        confirmed_deaths = missing,
        recovered_cases = missing,
        deaths_history = _days_only(obs.deaths_history),
        reported_history = _days_only(obs.reported_history),
        suspected_daily_history = _days_only(obs.suspected_daily_history),
        suspected_daily_deaths_history =
        _days_only(obs.suspected_daily_deaths_history),
        isolation_history = _days_only(obs.isolation_history),
        bed_capacity_history = _days_only(obs.bed_capacity_history),
        recovered_history = _days_only(obs.recovered_history),
        confirmed_history = obs.confirmed_history,
        confirmed_deaths_history = _days_only(obs.confirmed_deaths_history),
        lab_history = obs.lab_history,
        lab_daily_history = obs.lab_daily_history,
        export_case_days = obs.export_case_days,
        export_death_days = obs.export_death_days,
        breakpoint = _BREAKPOINT,
        background_re = true,
        confirmed_positivity_link = :composition,
        genetic = genetic_seeding_model,
        tmrca_days = obs.tmrca_days),
    chn_joint);

# `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 = (;
    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 (4-7 June) pick up where the cumulative suspected panel freezes on
# 26 May.
suspected_daily_panel = (;
    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 daily count, 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. 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.
isolation_panel = (;
    title = "Patients in isolation",
    dates = _vintage_dates(obs.isolation_history.days),
    replicates = _vintage_replicates(
        pp_joint, @varname(isolation.obs)),
    observed = obs.isolation_history.counts,
    colour = :darkorange, cumulative = false);
deaths_panel = (;
    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
# (7-14 June) pick up where the cumulative suspected-death panel freezes on
# 26 May, the deaths analogue of the new-suspects-per-day panel.
suspected_daily_deaths_panel = (;
    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 = (;
    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 = (;
    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 = (;
    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 = (;
    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 = (;
    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);

# Each panel runs to its own last vintage: the suspected case and death
# streams freeze at 26 May (their last stable vintage) while the
# laboratory-confirmed streams keep reporting to the cut-off, so the
# confirmed panels show the full series the model is fitting, not just the
# window the suspected streams cover.
joint_vintage_ppc_fig = plot_vintage_conditional_ppc(
    [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]);

The same check as per-vintage incidence: the count reported between consecutive situation reports rather than the running cumulative. Plotting the increment makes the trend in each stream read directly off the height of each step, so a rise or a slowdown is visible where the near-straight cumulative line hides it. The replicates are the modelled per-vintage increments, shown as 30/60/90% credible ribbons with the observed increment overlaid.

Per-vintage incidence posterior predictive plot

joint_vintage_incidence_fig = plot_vintage_incidence_ppc(
    [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]);

The exports group is checked next. The Uganda export and export-death streams are dated per-day series, each import or death scored as a Poisson at its detection day. The scalar posterior predictive sums each replicate's per-day count vector across the dated days, giving the cumulative export and death total to compare with the observed count.

Scalar posterior predictive plot

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

Counterfactual: lower bound under no further transmission

The committed future deaths $\mathrm{\Delta }D$ if transmission stopped at the report date, defined in the methods counterfactual.

Project no-onward deaths and summarise

no_onward = predict_no_onward_deaths(
    chn_joint; obs_deaths = obs.total_deaths);

no_onward_table = streams_table(
    "no-onward total" => no_onward.total_projected;
    digits = 0);

1×7 DataFrame

Row	Stream	Lower 90%	Lower 60%	Lower 30%	Upper 30%	Upper 60%	Upper 90%
	String	Float64	Float64	Float64	Float64	Float64	Float64
1	no-onward total	246.0	246.0	246.0	246.0	246.0	495.0

Two panels: the still expected deaths $\mathrm{\Delta }D$ (future deaths in cases already infected by $T$, net of those already observed) on the left, and the projected total $D(T)+\mathrm{\Delta }D$ on the right with the observed death count marked as a dashed black rule.

No-onward projected-deaths plot

no_onward_fig = plot_no_onward_deaths(
    no_onward; obs_deaths = obs.total_deaths);

Confirmed case-fatality ratio

The delay-corrected confirmed CFR defined in the methods delay-corrected confirmed CFR, set against the structural (infection-based) CFR and the naive confirmed ratio. The corrected ratio debiases the naive confirmed ratio for the real-time delay between a case being confirmed and a death being confirmed; the structural CFR is the onset-level estimate the joint model fits. Reading the three together separates the real-time delay bias (naive versus corrected) from the case/death ascertainment difference (corrected versus structural).

Compute the confirmed-CFR comparison

confirmed_cfr = delay_corrected_confirmed_cfr(chn_joint;
    obs_confirmed = obs.confirmed_cases,
    obs_confirmed_deaths = obs.confirmed_deaths);

confirmed_cfr_summary = confirmed_cfr_table(confirmed_cfr);

# Summary line carrying the data-anchored corrected estimate alongside the
# infection-based structural CFR, so both can be quoted together.
confirmed_cfr_line = let r = confirmed_cfr
    pct(x) = round(100 * x; digits = 1)
    corr = filter(isfinite, r.corrected)
    struc = filter(isfinite, r.structural)
    cs = posterior_summary(corr)
    ss = posterior_summary(struc)
    Markdown.parse(string(
        "**Delay-corrected confirmed CFR:** ",
        pct(quantile(corr, 0.5)), "% (90% CrI ",
        pct(cs.lo90), "–", pct(cs.hi90), "%), versus a naive confirmed ratio ",
        "of ", pct(r.naive_observed), "% and a structural (infection-based) ",
        "CFR of ", pct(quantile(struc, 0.5)), "% (90% CrI ",
        pct(ss.lo90), "–", pct(ss.hi90), "%)."))
end;

Delay-corrected confirmed CFR: 41.9% (90% CrI 30.7–60.2%), versus a naive confirmed ratio of 25.5% and a structural (infection-based) CFR of 41.5% (90% CrI 27.5–57.3%).

The four quantities side by side: the delay-corrected confirmed CFR, the structural CFR, the uncorrected modelled confirmed ratio and the naive observed confirmed ratio.

4×3 DataFrame

Row	Quantity	Central estimate	Narrowest interval
	String	String	String
1	Delay-corrected confirmed CFR	41.9%	30.7–60.2%
2	Structural (infection-based) CFR	41.5%	27.5–57.3%
3	Uncorrected modelled confirmed ratio	29.4%	22.2–39.5%
4	Naive observed confirmed ratio	25.5%	—

The posterior densities of the delay-corrected confirmed CFR and the structural CFR, with the naive observed confirmed ratio drawn as a solid vertical rule and the median uncorrected modelled confirmed ratio as a dashed rule. The gap from the naive rule to the corrected density is the real-time delay debiasing; the gap to the structural density is the residual case/death ascertainment difference.

Confirmed-CFR density plot

confirmed_cfr_fig = plot_confirmed_cfr(confirmed_cfr);

One-week-ahead forecast results

The cumulative and new expected counts by $T+7$ for the two confirmed DRC streams (laboratory-confirmed cases and confirmed deaths), from the no-change projection defined in the methods one-week-ahead forecast. The suspected reported cases and deaths are no longer reported, so they are not shown as forecast targets. The forecast also projects the isolation/treatment beds — both the bed demand a week ahead (the need under unconstrained supply) and the supply-limited occupancy it produces, whose gap is the projected bed shortfall — and the cumulative recovered total.

Generate the one-week-ahead forecast

forecast = forecast_reported(chn_joint;
    horizon = 7,
    obs_cases = obs.reported_cases,
    obs_deaths = obs.total_deaths,
    obs_confirmed = obs.confirmed_cases,
    obs_confirmed_deaths = obs.confirmed_deaths,
    obs_recovered = obs.recovered_cases);
forecast_summary = forecast_table(forecast);

8×8 DataFrame

Row	Stream	Quantity	Lower 90%	Lower 60%	Lower 30%	Upper 30%	Upper 60%	Upper 90%
	String	String	Float64	Float64	Float64	Float64	Float64	Float64
1	DRC confirmed cases	cumulative by T+7	66.0	256.0	398.0	748.0	1028.0	1837.0
2	DRC confirmed cases	new this week	0.0	0.0	0.0	0.0	0.0	789.0
3	DRC confirmed deaths	cumulative by T+7	16.0	71.0	117.0	220.0	305.0	553.0
4	DRC confirmed deaths	new this week	0.0	0.0	0.0	0.0	38.0	286.0
5	DRC isolation beds	demand at T+7	33.0	133.0	191.0	286.0	360.0	579.0
6	DRC isolation beds	occupancy at T+7	33.0	133.0	191.0	286.0	360.0	437.0
7	DRC recovered	cumulative by T+7	6.0	26.0	45.0	98.0	147.0	298.0
8	DRC recovered	new this week	0.0	0.0	0.0	0.0	35.0	186.0

The coming week at a glance, split into the latent quantities and the observations. The latent figure shows the new infections, symptom onsets and deaths over the horizon, with the reproduction number left to keep evolving across it.

One-week-ahead latent forecast plot

forecast_latent_fig = plot_forecast_latent(forecast);

The observation figure shows the new confirmed cases and confirmed deaths over the horizon.

One-week-ahead observed forecast plot

forecast_fig = plot_forecast(forecast);

The bed figure shows the projected isolation/treatment-bed demand (the need a week ahead, under unconstrained supply) against the supply-limited occupancy the beds can actually meet; the gap between the two is the projected bed shortfall, shown in the right panel. The reported "Patients en isolement" count is the occupied-bed count (the report computes the "Taux d'occupation" as that count over the bed capacity), so isolation is bed usage, gated by supply; the demand is its unobserved counterpart, the number who need a bed. Because the model carries a single national bed capacity it cannot represent local saturation, so the national shortfall understates local unmet need. On 13 June Ituri was at 93.9% occupancy while Sud-Kivu was at 21.9%, and beds free in one province cannot serve patients in another.

One-week-ahead isolation-bed forecast plot

forecast_beds_fig = plot_forecast_beds(forecast);

Forecast validation (last week versus now)

How last week's forecast held up against the data since observed, using the frozen re-fit and one-week projection defined in the methods forecast-versus-frozen evaluation. 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, widening the projected bed interval.

Fit one week back and validate the one-week-ahead forecast

# frozen_lastweek is computed in the setup block above.
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);

# 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.
_obs_beds = isempty(obs.isolation_history.counts) ? missing :
            obs.isolation_history.counts[end]
validation_table = forecast_vs_truth(validation_forecast;
    confirmed = obs.confirmed_cases,
    confirmed_deaths = obs.confirmed_deaths,
    isolation = _obs_beds);

3×9 DataFrame

Row	Stream	Observed	Lower 90%	Lower 60%	Lower 30%	Upper 30%	Upper 60%	Upper 90%	Within 90% PI
	String	Float64	Float64	Float64	Float64	Float64	Float64	Float64	String
1	DRC confirmed cases	1048.0	304.0	573.0	797.0	1258.0	1634.0	2547.0	yes
2	DRC confirmed deaths	267.0	57.0	127.0	173.0	293.0	385.0	615.0	yes
3	DRC isolation beds	371.0	333.0	406.0	419.0	436.0	447.0	465.0	yes

The observation panels histogram the one-week-ahead cumulative forecast made from the frozen fit, with the 90% predictive interval shaded and the count actually observed by the current cut-off drawn as a dashed black rule.

Forecast-versus-observed plot

validation_fig = plot_forecast_vs_truth(validation_forecast;
    confirmed = obs.confirmed_cases,
    confirmed_deaths = obs.confirmed_deaths,
    baseline_confirmed = frozen_lastweek.o.confirmed_cases,
    baseline_confirmed_deaths = frozen_lastweek.o.confirmed_deaths);

The bed panel scores last week's projected isolation-bed occupancy against the beds actually occupied now (the dashed rule). At a one-week-back freeze the capacity has no implied-capacity anchor (the occupancy rate starts only on 9 June), so the projection rides the capacity random walk back to the freeze date and its interval is wide.

Bed forecast-versus-observed plot

validation_beds_fig = plot_forecast_beds_vs_truth(validation_forecast;
    isolation = _obs_beds);

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

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

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 on its own and what the joint adds.

Per-stream infection-count table

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,
    "joint" => posterior_C_joint);

6×7 DataFrame

Row	Stream	Lower 90%	Lower 60%	Lower 30%	Upper 30%	Upper 60%	Upper 90%
	String	Float64	Float64	Float64	Float64	Float64	Float64
1	exports	1068.0	2084.0	3096.0	6673.0	11801.0	45024.0
2	deaths (DRC)	5730.0	7430.0	8667.0	11533.0	13926.0	20008.0
3	cases (DRC)	7594.0	8691.0	9592.0	11843.0	14250.0	22409.0
4	confirmed (DRC)	10113.0	12312.0	14254.0	18861.0	23120.0	36074.0
5	isolation (DRC)	2237.0	3191.0	4021.0	6247.0	8828.0	17324.0
6	joint	1632.0	1896.0	2104.0	2623.0	3024.0	3838.0

Each single-stream fit projects its outbreak size out to the cut-off, even for streams whose data stops earlier. The first figure shows each fit's cumulative-infection trajectory over the grid as 50% and 90% credible ribbons, with a dotted vertical rule in each stream's colour where that stream's data stops reporting. The suspected case and death streams freeze at 26 May while the exports and confirmed streams run on, so the part of each ribbon beyond its rule is the model projecting forward from the last data it saw.

Per-stream projected-trajectory plot

# 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)];
    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 confirmed-cases-only stream is ill-defined on its own, so its posterior runs far wider than the rest. The horizontal axis is scaled to a multiple of the joint-fit 90% upper bound, the estimate that constrains every stream together, so the bulk of the joint and the other streams stays visible rather than being flattened by the confirmed-only tail.

Cut-off infection-count density plot

# 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,
    "joint" => posterior_C_joint;
    scenarios = [], xmax = density_xmax);

The third figure is the reproduction number each stream implies on its own, one panel per stream with the joint fit overlaid in grey as the reference. Each panel draws 30/60/90% credible ribbons with no median line, matching the band style used elsewhere; the window, the response markers and the cut-off match the joint Rt figure above.

Per-stream implied-Rt plot

# 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)];
    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 = 21.0);

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), reusing the frozen-fit helper defined above.

Freeze the renewal data to a cut-off and re-fit

# Frozen re-fits and released_df are prepared in the setup block above.

Estimate evolution across releases

How the outbreak-size estimate moved as the situation reports accrued.

The project publishes a tagged results release at each data cut-off (https://github.com/epiforecasts/BVDOutbreakSize/releases), bundling the posterior draws and input data. The released series, in blue, is the project's published estimate at each release: the closed-form integral model up to v1.3.0, then the renewal model from v1.4.0 on. Each release is its own fit, so it is drawn as a discrete estimate, a median with nested 30/60/90% interval bars, rather than a ribbon. The renewal series, in red, is the renewal model re-fit frozen at each integral-era release cut-off. The renewal-era releases already are renewal fits, so they carry no frozen re-fit. The current-data, current-model estimate is drawn in green as the cumulative-infection trajectory over time, a single fit shown across the period so the latest estimate reads against the earlier ones. Each release date is marked with a dotted vertical rule.

Released estimates and frozen renewal re-fits

# 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 renewal model re-fit frozen at each integral-era release
# cut-off, each its own discrete estimate. Each tuple carries the median
# and 30/60/90% credible bounds.
function _ci369(xs)
    q(p) = round(Int, 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
# Reuse the one-week-back frozen fit (already run for forecast validation)
# as an additional recent renewal point, so the evolution plot shows the
# current-vintage frozen estimate without an extra re-fit.
frozen_by_cutoff[validation_cutoff] = frozen_lastweek
renewal_frozen = [(c, _ci369(frozen_C(c))...)
                  for c in sort(union(frozen_evolution_cutoffs,
    [validation_cutoff]))]

# 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-renewal 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 = renewal_frozen,
    trajectory = infection_trajectory,
    title = "Outbreak-size estimate as data accrued");

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 both their reports, the 18 May report and the 20 May update, with those intervals. The two reports share the same geographic-spread scenarios from exported cases and travel volume, so those appear once. 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, with a corrected set of 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, the symptomatic cases, on the same dates, rather than 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, 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, so each scenario sits beside our estimate for the date it was made.

McCabe scenarios with uncertainty against our estimates

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

# McCabe scenarios with their reported 95% confidence intervals, grouped by
# report date, then our modelled cumulative onsets on the matching
# report date with its credible interval.
mccabe_rows = [(label, mean, lo, hi)
               for (_, label, mean, lo, hi) in REPORT_SCENARIOS_CI]
mccabe_groups = [date == "2026-05-18" ? "McCabe et al. (18 May)" :
                 date == "2026-05-20" ? "McCabe et al. (20 May update)" :
                 "McCabe et al. (27 May, Lancet)"
                 for (date, _, _, _, _) in REPORT_SCENARIOS_CI]

ours_rows = [("Renewal onsets on 18 May", _ours_on("2026-05-18")...),
    ("Renewal onsets on 20 May", _ours_on("2026-05-20")...),
    ("Renewal onsets on 27 May", _ours_on("2026-05-27")...)]
ours_groups = fill("Our renewal estimate", 3)

matched_rows = vcat(mccabe_rows, ours_rows)
matched_groups = vcat(mccabe_groups, ours_groups)

matched_comparison_fig = plot_estimate_comparison(matched_rows;
    xlabel = "Cumulative cases / infections",
    groups = matched_groups,
    group_colours = ["McCabe et al. (18 May)" => :grey,
        "McCabe et al. (20 May update)" => :black,
        "McCabe et al. (27 May, Lancet)" => :steelblue,
        "Our renewal estimate" => :firebrick]);

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.

Side-by-side outbreak-size intervals for the two frozen fits and the current-data fit, so the shift with the data cut-off reads off directly.

Frozen-fit C_T table

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"),
    "current data" => posterior_C_joint);

4×7 DataFrame

Row	Stream	Lower 90%	Lower 60%	Lower 30%	Upper 30%	Upper 60%	Upper 90%
	String	Float64	Float64	Float64	Float64	Float64	Float64
1	frozen 20 May	797.0	1105.0	1338.0	1885.0	2272.0	3372.0
2	frozen 23 May	1099.0	1580.0	1868.0	2420.0	2855.0	4023.0
3	frozen 27 May	280.0	379.0	450.0	602.0	720.0	1028.0
4	current data	1632.0	1896.0	2104.0	2623.0	3024.0	3838.0

Delay sensitivity

Paused due to compute constraints. This sensitivity re-fit is not run in the current build to keep the docs build within compute limits. The method and code are unchanged and it can be re-enabled (see the sensitivity gate in the setup block); the description below documents what it does.

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, which 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

# refit_joint_variant, deaths_community_delay and the sensitivity re-fits are
# defined and run (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

delay_sensitivity_table = RUN_SENSITIVITY ?
                          streams_table("baseline (hospital pathway)" => posterior_C_joint,
    "community pathway" => posterior_C_community_delay) :
                          Markdown.md"_Delay sensitivity is paused due to compute constraints._"

Delay sensitivity is paused due to compute constraints.

Delay sensitivity is paused due to compute constraints.

Delay-sensitivity infection-count density plot

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 is paused due to compute constraints._"

Delay sensitivity is paused due to compute constraints.

Delay sensitivity is paused due to compute constraints.

Clock-rate sensitivity

Paused due to compute constraints. This sensitivity re-fit is not run in the current build to keep the docs build within compute limits. The method and code are unchanged and it can be re-enabled (see the sensitivity gate in the setup block); the description below documents what it does.

The whole outbreak-age estimate rests on the genetic bound, the oldest date the common ancestor of the sequenced cases can sit, which is set by the assumed molecular clock rate. The baseline uses the slower clock rate of $1.2\times {10}^{-3}$ substitutions per site per year, the rate of the 2013-2016 West African Ebola epidemic, which dates the common ancestor to 25 March 2026. The sequencing source also reports a faster early-epidemic rate of $1.9\times {10}^{-3}$ substitutions per site per year, which dates the common ancestor about two and a half weeks more recently, to 11 April 2026, without favouring either (Amuri-Aziza et al., 2026). We re-fit the joint model under the faster clock and compare the infection count to date and the outbreak age. The re-fit uses the full headline settings (1000 draws across two chains).

Re-fit the joint under the faster clock rate

# clock_alt_offset, tmrca_days_alt and the faster-clock re-fit are defined and
# run (when enabled) in the setup block above.
posterior_C_fast_clock = RUN_SENSITIVITY ?
                         vec(Array(chn_joint_fast_clock[:C_T])) : nothing
T_baseline_clock = vec(Array(chn_joint[:T]))
T_fast_clock = RUN_SENSITIVITY ? vec(Array(chn_joint_fast_clock[:T])) : nothing

The infection count to date under the two clock rates, side by side.

Clock-rate infection-count table

clock_sensitivity_C_table = RUN_SENSITIVITY ?
                            streams_table("baseline clock" => posterior_C_joint,
    "faster clock" => posterior_C_fast_clock) :
                            Markdown.md"_Clock-rate sensitivity is paused due to compute constraints._"

Clock-rate sensitivity is paused due to compute constraints.

Clock-rate sensitivity is paused due to compute constraints.

Clock-rate infection-count density plot

clock_sensitivity_C_fig = RUN_SENSITIVITY ?
                          plot_cumulative_cases("baseline clock" => posterior_C_joint,
    "faster clock" => posterior_C_fast_clock; scenarios = []) :
                          Markdown.md"_Clock-rate sensitivity is paused due to compute constraints._"

Clock-rate sensitivity is paused due to compute constraints.

Clock-rate sensitivity is paused due to compute constraints.

The outbreak age, the number of days from seeding to the cut-off, under the two clock rates. A more recent common ancestor permits a younger outbreak.

Clock-rate outbreak-age table

clock_sensitivity_T_table = RUN_SENSITIVITY ?
                            streams_table("baseline clock" => T_baseline_clock,
    "faster clock" => T_fast_clock; digits = 0) :
                            Markdown.md"_Clock-rate sensitivity is paused due to compute constraints._"

Clock-rate sensitivity is paused due to compute constraints.

Clock-rate sensitivity is paused due to compute constraints.

Clock-rate outbreak-age density plot

clock_sensitivity_T_fig = RUN_SENSITIVITY ?
                          plot_density_overlay("baseline clock" => T_baseline_clock,
    "faster clock" => T_fast_clock;
    xlabel = "Outbreak age (days before cut-off)",
    title = "Posterior outbreak age by clock rate") :
                          Markdown.md"_Clock-rate sensitivity is paused due to compute constraints._"

Clock-rate sensitivity is paused due to compute constraints.

Clock-rate sensitivity is paused due to compute constraints.

Saving results

The tables above are written to an output/ directory at the repo root so they can be archived and shared. On every push to main a GitHub Actions workflow regenerates these files and publishes them as a GitHub Release, downloadable from the repository's releases page (https://github.com/epiforecasts/BVDOutbreakSize/releases). The release bundles the summary tables, a thinned set of posterior draws, the latent symptom-onset ("symptomatic cases") trajectory over time, and a copy of the input observations.toml so the exact data that produced each result is recorded alongside it.

Write outputs to output/

# Outputs default to `output/` in the package directory (where the
# docs build and Release workflow expect them). Set `BVD_OUTPUT_DIR`
# to redirect them, e.g. when running from a read-only package
# install.
output_dir = get(ENV, "BVD_OUTPUT_DIR",
    joinpath(pkgdir(BVDOutbreakSize), "output"))
mkpath(output_dir)

# Full parameter summary for the published CSV (infection, surveillance and
# export parameters together).
joint_summary = summary_table(chn_joint,
    [:r, :r0, :doubling_time, :T, :R_T, :CFR, :C_T,
        :p_drc, :p_uganda, :k, :tau_test, :lambda_bg,
        Symbol("exports_state.travel_state.daily_travellers")]; digits = 2)
CSV.write(joinpath(output_dir, "posterior_summary.csv"), joint_summary)
CSV.write(joinpath(output_dir, "confirmed_cfr_summary.csv"),
    confirmed_cfr_summary)
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)

# Copy the input data so the release records what produced these
# results.
cp(joinpath(pkgdir(BVDOutbreakSize), "data", "observations.toml"),
    joinpath(output_dir, "observations.toml"); force = true)

# Thinned posterior draws of the key joint parameters (every 10th
# draw) so downstream users can recompute their own summaries.
# `cumulative_onsets_T` is the cumulative symptom onsets by the cut-off,
# the latent "symptomatic cases" outcome (the onset analogue of `C_T`),
# read off the last day of each draw's `cumulative_onsets` trajectory.
_cum_onset_draws = vec(collect(chn_joint[:cumulative_onsets]))
cumulative_onsets_T = Float64[v[end] for v in _cum_onset_draws]
posterior_draws = DataFrame(
    r = vec(Array(chn_joint[:r])),
    r0 = vec(Array(chn_joint[:r0])),
    doubling_time = vec(Array(chn_joint[:doubling_time])),
    T = vec(Array(chn_joint[:T])),
    R_T = vec(Array(chn_joint[:R_T])),
    CFR = vec(Array(chn_joint[:CFR])),
    p_drc = vec(Array(chn_joint[:p_drc])),
    p_uganda = vec(Array(chn_joint[:p_uganda])),
    C_T = vec(Array(chn_joint[:C_T])),
    cumulative_onsets_T = cumulative_onsets_T,
    confirmed_cfr_corrected = confirmed_cfr.corrected
)[1:10:end, :]
CSV.write(joinpath(output_dir, "posterior_draws.csv"), posterior_draws);

# Latent symptom-onset trajectory over time, the "symptomatic cases" curve,
# showing outbreak growth: one row per grid day with the 30/60/90%
# credible intervals of both the daily new and cumulative onsets.
onsets_over_time_table = onsets_over_time(chn_joint;
    n = obs.n, seeding = obs.seeding)
CSV.write(joinpath(output_dir, "onsets_over_time.csv"),
    onsets_over_time_table);

Summary-page assets

The one-page Summary dashboard reuses the results computed above rather than re-fitting. Here we save its headline text, headline tables and the four figures it shows (reproduction number, the reproduction number each data stream implies on its own, infections over time, and modelled versus observed reported cases) into docs/src/summary_assets/, so the static dashboard page can embed them after this build step has run.

Write the dashboard assets

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

# Figures: estimated R(t), the R(t) each data stream implies on its own,
# latent infections over time, and the modelled versus observed reported
# cases. All are produced in the Results sections above; here we just write
# them out at the dashboard size.
CairoMakie.save(joinpath(dashboard_dir, "rt.png"), rt_fig)
CairoMakie.save(joinpath(dashboard_dir, "rt_streams.png"), stream_rt_fig)
CairoMakie.save(joinpath(dashboard_dir, "infections.png"),
    cumulative_traj_fig)
CairoMakie.save(joinpath(dashboard_dir, "reported_cases.png"),
    joint_vintage_ppc_fig)

# Headline prose: the same bullet summary shown at the top of the Results
# section, serialised to markdown so the dashboard renders it verbatim.
open(joinpath(dashboard_dir, "headline.md"), "w") do io
    print(io, sprint(Markdown.plain, summary_ranges))
end

# Headline tables: outbreak size and timing as whole numbers, and the
# growth and severity parameters to two decimals, each with reader-friendly
# quantity names.
dashboard_counts = summary_table(chn_joint, [:C_T, :T]; digits = 0,
    labels = Dict(:C_T => "Cumulative infections",
        :T => "Outbreak age (days)"))
dashboard_rates = summary_table(chn_joint,
    [:R0, :R_T, :r, :doubling_time, :CFR]; digits = 2,
    labels = Dict(:R0 => "Initial reproduction number",
        :R_T => "Latest reproduction number",
        :r => "Latest growth rate (per day)",
        :doubling_time => "Latest doubling time (days)",
        :CFR => "Case-fatality ratio"))
open(joinpath(dashboard_dir, "headline_counts.md"), "w") do io
    print(io, markdown_table(dashboard_counts))
end
open(joinpath(dashboard_dir, "headline_rates.md"), "w") do io
    print(io, markdown_table(dashboard_rates))
end

# The data cut-off the dashboard reports as of, written as a plain date.
open(joinpath(dashboard_dir, "cutoff.md"), "w") do io
    print(io, string(obs.cutoff))
end

---

The full analysis code, data and model definitions are in the epiforecasts/BVDOutbreakSize repository. Issues, corrections and suggestions are welcome there. Maintained by Sam Abbott, Kath Sherratt, Samuel Brand and Sebastian Funk.