Estimating the current size of the 2026 DRC Bundibugyo virus outbreak

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

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

Data as of: 4 July 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. · comparison with Chamla 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. A forward projection from a frozen fit is also set against the Chamla et al. (Chamla et al., 2026) confirmed-case projection and the data observed since.

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.

• Single national bed capacity. The treatment-centre model carries 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% — and the national bed shortfall understates local unmet need.

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, data and fitted chains

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

45

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). From 13 June the reports add a Tableau 6 patient-movement table for the treatment centres, and we read its daily admissions, in-care deaths, rule-outs and absconded flows as four count streams feeding the same treatment-centre model. 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

observations_table = DataFrame(
    field = [
        "exported_cases",
        "exports_deaths",
        "suspected_deaths",
        "suspected_cases",
        "confirmed_cases",
        "confirmed_deaths",
        "specimens_analysed",
        "treatment_admissions",
        "treatment_deaths",
        "treatment_ruleouts",
        "treatment_absconded",
        "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),
        hist_last_date(obs.treatment_admissions_history),
        hist_last_date(obs.treatment_deaths_history),
        hist_last_date(obs.treatment_ruleout_history),
        hist_last_date(obs.treatment_absconded_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,
        isempty(obs.treatment_admissions_history.counts) ? missing :
        obs.treatment_admissions_history.counts[end],
        isempty(obs.treatment_deaths_history.counts) ? missing :
        obs.treatment_deaths_history.counts[end],
        isempty(obs.treatment_ruleout_history.counts) ? missing :
        obs.treatment_ruleout_history.counts[end],
        isempty(obs.treatment_absconded_history.counts) ? missing :
        obs.treatment_absconded_history.counts[end],
        obs.tmrca_days,
        ITURI_DAILY_TRAVEL,
        ITURI_DAILY_TRAVEL_SD,
        ITURI_POPULATION
    ]
);

15×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-07-04	1561
6	confirmed_deaths	2026-07-04	506
7	specimens_analysed	2026-05-28	755
8	treatment_admissions	2026-07-02	72
9	treatment_deaths	2026-07-02	16
10	treatment_ruleouts	2026-07-02	55
11	treatment_absconded	2026-07-02	1
12	genetic_tmrca_bound	2026-03-25	101
13	daily_outbound_travellers (prior mean)	missing	1871
14	daily_outbound_travellers_sd (prior SD)	missing	200
15	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

47×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
37	2026-06-22	missing	131	387	missing	1094	277	115	missing	missing
38	2026-06-23	missing	138	408	missing	1118	291	122	missing	missing
39	2026-06-24	missing	154	385	missing	1155	304	138	missing	missing
40	2026-06-25	missing	265	419	missing	1203	321	148	missing	missing
41	2026-06-27	missing	239	502	missing	1274	360	178	missing	missing
42	2026-06-29	missing	309	609	missing	1333	399	189	missing	missing
43	2026-06-30	missing	301	missing	missing	1406	438	208	missing	missing
44	2026-07-01	missing	150	641	missing	1460	452	213	missing	missing
45	2026-07-02	missing	213	628	missing	1502	473	229	missing	missing
46	2026-07-03	missing	185	missing	missing	1528	492	239	missing	missing
47	2026-07-04	missing	354	missing	missing	1561	506	missing	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.6),\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.6); lower = 0),
        centred::Bool = true)
    μ_log ~ mean_prior
    τ ~ sd_prior
    m = max(n_streams, 1)
    if centred
        ## Draw each stream's `log(1/sqrt(k))` directly from the population.
        ## `eps` floors the SD so a `τ ≈ 0` draw stays a proper distribution.
        log_isk ~ product_distribution(
            fill(Normal(μ_log, τ + eps(typeof(τ))), m))
        inv_sqrt_k = exp.(log_isk[1:n_streams])
    else
        z ~ product_distribution(fill(Normal(0, 1), m))
        inv_sqrt_k = exp.(μ_log .+ τ .* z[1:n_streams])
    end
    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,
        ## Observed daily contact-tracing follow-up rate ("taux de suivi des
        ## contacts", a fraction) as `(; days, values)`. When present it drives
        ## a shared LATENT follow-up-rate process ([`contact_tracing_model`](@ref))
        ## defined over the whole grid, so the case-finding anchor is inferred
        ## before the reported window and projected into forecasts rather than
        ## stopping with the data. Empty (the default) leaves the anchor off.
        contact_followup_history = (; days = Int[], values = Float64[]),
        contact_tracing = contact_tracing_model,
        ## Non-BVD background elasticity to the contact rate (detection scales
        ## the non-BVD background: `exp(β_contact · (q − q̄))`).
        contact_effect = contact_background_model(),
        ## Optional suspected-case reporting-effort (BVD ascertainment) submodel
        ## constructor ([`reporting_effort_walk_model`](@ref)), anchored to the
        ## same latent contact rate plus a walk and applied to the BVD suspected
        ## component only. `nothing` (the default) keeps a constant effort.
        reporting_effort = 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

    ## Shared LATENT contact-tracing follow-up rate. When the follow-up-rate
    ## series is supplied it drives a logistic random-walk process
    ## ([`contact_tracing_model`](@ref)) defined over the whole grid `q_t ∈
    ## (0, 1)`, fitted to the observed rate at its report days. Both the non-BVD
    ## background and the BVD reporting effort anchor to this ONE rate, so
    ## case-finding intensity is carried consistently and continues under the
    ## walk where the data stop. Empty history → no latent process (anchor off).
    has_contact = !isempty(contact_followup_history.days)
    if has_contact
        ct_days = [clamp(Int(d), 1, n) for d in contact_followup_history.days]
        ct_logit = [logit(clamp(Float64(v), 1.0e-3, 1.0 - 1.0e-3))
                    for v in contact_followup_history.values]
        contact_state ~ to_submodel(
            contact_tracing(n; obs_days = ct_days, logit_obs = ct_logit))
        q_contact = contact_state.q
        q_mean = sum(q_contact) / n
    else
        q_contact = nothing
    end

    ## Non-BVD background scaled by contact-tracing intensity. Contact tracing
    ## is a detection process, so it SCALES the non-BVD background rather than
    ## adding to it: `bg_t · exp(β_contact · (q_t − q̄))`, mean-centred so
    ## `β_contact` is an elasticity around mean intensity and the level stays
    ## with λ_bg / the walk. It scales only the non-BVD background (not
    ## `p_drc`), so it does not reopen the ascertainment / outbreak-size
    ## degeneracy a latent-scaled background would (issue #374).
    if q_contact === nothing
        β_contact = zero(λ_bg_base)
    else
        contact_bg_state ~ to_submodel(contact_effect)
        β_contact = contact_bg_state.β_contact
        bg_daily = bg_daily .* exp.(β_contact .* (q_contact .- q_mean))
    end

    ## Suspected daily cases add the p_drc-scaled BVD signal and the non-BVD
    ## background. With a reporting-effort submodel injected, the BVD component
    ## carries a case-finding-effort multiplier ([`reporting_effort_walk_model`]
    ## (@ref)) anchored to the SAME latent contact rate plus a walk, applied
    ## HERE only (not to the returned `bvd_reports_daily` the laboratory /
    ## treatment / death streams reuse), so it is a suspected-case-specific
    ## ascertainment term. `nothing` leaves the constant-ascertainment form.
    if reporting_effort === nothing
        effort_series = nothing
        reports_daily = p_drc .* bvd_reports_daily .+ bg_daily
    else
        effort_state ~ to_submodel(
            reporting_effort(n; contact_covariate = q_contact))
        effort_series = effort_state.effort
        reports_daily = (p_drc .* effort_series) .* bvd_reports_daily .+ bg_daily
    end

    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, β_contact, reporting_effort = effort_series,
        contact_rate = q_contact)
end

Treatment-centre flow

The treatment-centre stream models the daily patient flow through the isolation/treatment centres (CTE/CT/CI): the occupied-bed count ("Patients en isolement"), the daily admissions, and the daily discharges split by outcome — in-care deaths ("décédés"), rule-outs ("non-cas") and absconded ("évadés") — read from the situation-report Tableau 6 patient-movement table. Recoveries ("cumul guéris") are the confirmed-and-discharged subset and are modelled as a separate confirmed-recovery stream (below).

Beds can be supply-driven, with demand outstripping supply and occupancy catching up as capacity expands, so the occupied beds are the suspect admissions carried through a length-of-stay survival $S(\tau )=P(\text{LOS}\ge \tau )$ into a supply-limited occupancy (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, splitting into BVD and non-BVD admissions that leave on different clocks, so the latent bed demand sums two survival convolutions,

$$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)],$$

where the BVD stay $S_{\text{BVD}}$ is itself an outcome mixture: an admitted BVD patient leaves by death (weight ${\text{CFR}}_{\text{iso}}$, the admission-to-death stay) or by recovery (weight $1-{\text{CFR}}_{\text{iso}}$, a longer admission-to-recovery stay). The death-stay prior is the admission-to-death delay from the line-list reanalysis (Funk and Abbott, 2026), and the non-BVD rule-out stay $S_{\text{ruleout}}$ takes the report-to-receipt laboratory turnaround.

The in-care fatality is a sampled log-odds modifier ${\beta }_{\text{iso}}$ on the infection CFR,

$${\text{CFR}}_{\text{iso}}={\mathrm{logit}}^{-1}(\mathrm{logit}\text{CFR}+{\beta }_{\text{iso}}),$$

identified by the in-care death flow. It is an in-care fatality, sitting below the infection CFR where treatment reduces mortality, and is reported with ${\beta }_{\text{iso}}$ and the overall length-of-stay (the mixture mean).

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 beds are scored as the latent demand right-censored at the recorded implied capacity $C_j^{\text{cap}}$ (so the demand above a saturated capacity is left uncensored), and the daily admissions and discharges (in-care deaths, rule-outs, absconded) are scored as additional count streams $F_j$ sharing the dispersion $k_{\text{iso}}$,

$$O_j\sim \mathrm{censored}(\mathrm{NegBinomial}(D_{t_j},k_{\text{iso}});\text{upper}=C_j^{\text{cap}}),F_j\sim \mathrm{NegBinomial}({\mu }_{t_j}^F,k_{\text{iso}}),$$

with each ${\mu }_t^F$ the matching branch of the demand (the BVD and non-BVD inflow, ${\text{CFR}}_{\text{iso}}$ of BVD admissions through the death stay, the non-BVD admissions through the rule-out stay, and a small fraction of occupancy), and the implied capacity carried by a NegBinomial of its own.

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$.

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.

A small set of curated days carry an overnight reclassification of who is counted in the occupancy. On such a day the reported start-of-day in-bed count ("Patients au lit (J-1)") falls below the previous report's end-of-day occupancy ("Fin J") by more than that day's admissions and discharges can bridge — an operational de-registration of ruled-out suspects rather than a fall in transmission — which the smooth survival demand cannot reproduce. We read these days off the situation reports rather than flagging them automatically. Two are carried: 9 June, where the start-of-day stock 268 sits 29 below the previous day's 297 and the drop is concentrated in the suspected pool (dont suspects $184\to 143$ against dont confirmés $113\to 117$), and 19 June, where the start-of-day stock 342 sits 74 below the previous 416. A level step $b_j$ is fitted at each reclassification day $d_j$, and the cumulative offset added to the modelled occupancy on day $t$ is

$$\begin{array}{}\text{(26)}& {\mathrm{\Delta }}_t=\sum _{j:d_j\le t}{b}_j,b_j={\sigma }_{\text{brk}}{z}_j,z_j\sim \mathrm{Normal}(0,1),\end{array}$$

so the censored occupancy likelihood above scores $D_{t_j}+{\mathrm{\Delta }}_{t_j}$ in place of $D_{t_j}$. Each step is centred on zero with a weakly informative scale ${\sigma }_{\text{brk}}=25$ beds, so the fit partitions each reclassification into reporting artefact and genuine change in demand, and the persistent offset, carried forward to every later day, absorbs the overnight gap without bending the reproduction number to chase it.

Submodel: treatment_flow_model

@model function treatment_flow_model(
        isolation_history,
        bvd_reports_daily::AbstractVector,
        bg_daily::AbstractVector,
        p_drc::Real,
        CFR::Real;
        capacity_history = (; days = Int[], counts = Int[]),
        admissions_history = (; days = Int[], counts = Int[]),
        deaths_history = (; days = Int[], counts = Int[]),
        ruleout_history = (; days = Int[], counts = Int[]),
        absconded_history = (; days = Int[], counts = Int[]),
        ## Manually specified occupancy reclassification-break days (grid day-
        ## indices, opt-in). A level step is fitted into the modelled occupancy
        ## mean at each, absorbing a between-report measurement-basis
        ## discontinuity in the observed isolation series. Empty → no break
        ## days, a no-op (the default). See `cumulative_occupancy_offset`.
        occupancy_break_days::AbstractVector{<:Integer} = Int[],
        ## Prior sd of each occupancy break step (beds). Weakly informative and
        ## centred on zero, so the fit decides how much of the step is a
        ## reporting artifact vs real demand change.
        occupancy_break_sd::Real = 25.0,
        admission = isolation_admission_model(),
        severity = isolation_severity_model(),
        capacity = bed_capacity_walk_model,
        dispersion = surveillance_dispersion_model(),
        ## Occupancy / flow dispersion can be injected from the joint composer's
        ## pooled set (`k_external`); standalone it samples its own.
        k_external::Union{Nothing, Real} = nothing,
        ## In-care fatality modifier prior: β_iso on the infection CFR.
        cfr_modifier_prior = Normal(0.0, 0.5),
        ## Small abscond / loss-to-follow-up fraction of occupancy per day.
        abscond_prior = truncated(Normal(0.01, 0.01); lower = 0),
        ## Short suspected→admission delay (report → reaching a bed: triage,
        ## transport, bed-wait), distinct from the report→lab receipt delay.
        admission_delay = censored_delay_model(
            cdf_nmax(lognormal_meansd(2.0, 1.5); q = 0.99);
            mean_prior = truncated(Normal(2.0, 1.0); lower = 0.1),
            sd_prior = truncated(Normal(1.5, 1.0); lower = 0.3)),
        ## Outcome-mixture BVD bed stay: admission→death (the admission→death
        ## atomic delay the onset→death convolution also uses, mean ≈ 8.4 d) and
        ## the longer admission→recovery stay (mean ≈ 14 d). Built to a common
        ## nmax so the two PMFs align for the elementwise mixture.
        death_los = gamma_delay_model(
            cdf_nmax(lognormal_meansd(14.0, 8.0); 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)),
        recovery_los = 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)),
        ## Non-BVD rule-out stay (report→receipt turnaround plus sign-off).
        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 `p_iso` by the severity log-odds `δ_iso`.
    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)
    ## In-care fatality CFR_iso = logistic(logit(CFR) + β_iso): a log-odds
    ## modifier on the infection CFR, identified by the in-care death flow
    ## (Tableau 6 décédés) relative to admissions and occupancy. β_iso < 0 means
    ## treatment lowers the in-care fatality below the infection CFR. This is a
    ## conditional-on-admission (in-care) fatality, not a causal treatment
    ## effect. The recovered-among-confirmed ("cumul guéris") stream is modelled
    ## separately off the confirmed cases ([`recovered_model`](@ref)) — guéris is
    ## the confirmed-and-discharged subset, not all in-care recoveries.
    β_iso ~ cfr_modifier_prior
    CFR_iso = logistic(logit(CFR) + β_iso)
    ## Time-varying bed capacity `C(t)` (a random walk), started at the first
    ## day with occupancy or capacity data.
    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]
    adm_delay_state ~ to_submodel(admission_delay)
    death_los_state ~ to_submodel(death_los)
    recovery_los_state ~ to_submodel(recovery_los)
    ruleout_los_state ~ to_submodel(ruleout_los)
    abscond_frac ~ abscond_prior

    ## Admission inflow carried through the short suspected→admission delay,
    ## split into BVD (admitted at `p_iso_bvd`) and non-BVD (`p_iso`) demand.
    bvd_adm = convolve_delay(p_iso_bvd .* p_drc .* bvd_reports_daily,
        adm_delay_state.pmf)
    bg_adm = convolve_delay(p_iso .* bg_daily, adm_delay_state.pmf)

    ## Outcome-mixture (not competing risks) BVD bed stay: a fraction `CFR_iso`
    ## of BVD admissions leave on the admission→death stay, the rest on the
    ## admission→recovery stay. The two PMFs share an nmax, so the mixture is
    ## elementwise.
    dpmf = death_los_state.pmf
    rpmf = recovery_los_state.pmf
    f_bvd = CFR_iso .* dpmf .+ (one(CFR_iso) - CFR_iso) .* rpmf

    ## Latent bed demand: admissions carried through the survival of their
    ## length-of-stay (occupancy = admissions still in a bed).
    bvd_demand = convolve_survival(bvd_adm, f_bvd)
    bg_demand = convolve_survival(bg_adm, ruleout_los_state.pmf)
    demand = bvd_demand .+ bg_demand
    ## Predict / check-model mode can widen the series to `Vector{Any}`, which
    ## `min.(demand, C)` cannot broadcast over; pin to the capacity's concrete
    ## element type, leaving the AD/fit path untouched.
    if eltype(demand) === Any
        demand = convert(Vector{eltype(C)}, demand)
    end
    occupancy = min.(demand, C)

    ## Occupancy likelihood: each day's bed count is a NegativeBinomial around
    ## the latent demand, right-censored at the fixed implied-capacity bound.
    ## Day-to-day reporting noise is absorbed by the NegativeBinomial
    ## dispersion. A manually specified, opt-in cumulative reclassification
    ## offset Δ is added to the modelled mean only on the supplied
    ## `occupancy_break_days`, so the fit tracks a between-report measurement-
    ## basis discontinuity in the observed series without bending Rt to chase
    ## it. Each step is sampled non-centred and centred on zero, so the fit
    ## partitions it into reporting artifact vs real demand. Empty (the
    ## default) → no sampled step, Δ = 0, a no-op.
    iso_days = isolation_history.days
    iso_obs = isempty(isolation_history.counts) ? missing :
              collect(Int.(isolation_history.counts))
    iso_demand = [demand[clamp(Int(d), 1, n)] for d in iso_days]
    ## Only break days that fall on or before an observed occupancy day can move
    ## the likelihood; drop any later ones so no inert step is sampled.
    iso_last = isempty(iso_days) ? 0 : maximum(iso_days)
    brk_days = isempty(iso_days) ? Int[] :
               [Int(d) for d in occupancy_break_days if Int(d) <= iso_last]
    occ_break_val = zero(eltype(C))
    if isempty(brk_days)
        iso_means = iso_demand
    else
        occupancy_step ~ product_distribution(
            fill(Normal(0, 1), length(brk_days)))
        b = occupancy_break_sd .* occupancy_step
        occ_offset = cumulative_occupancy_offset(iso_days, brk_days, b)
        iso_means = iso_demand .+ occ_offset
        occ_break_val = isempty(occ_offset) ? occ_break_val : last(occ_offset)
    end
    iso_ceil = censoring_cap(iso_days, iso_obs, capacity_history)
    isolation ~ to_submodel(
        censored_occupancy_model(iso_means, iso_ceil, iso_obs, k))

    ## Capacity likelihood: the implied bed count is a noisy observation of C(t).
    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))

    ## Modelled daily in-care outcome flows. The death and recovery flows sum
    ## to the total BVD discharge (consistent with the mixture survival).
    deaths_daily = convolve_delay(CFR_iso .* bvd_adm, dpmf)
    recover_daily = convolve_delay((one(CFR_iso) - CFR_iso) .* bvd_adm, rpmf)
    ruleout_daily = convolve_delay(bg_adm, ruleout_los_state.pmf)
    admit_daily = bvd_adm .+ bg_adm
    abscond_daily = abscond_frac .* occupancy

    ## Optional daily Tableau 6 flow likelihoods. Each is a no-op when its
    ## history is empty (no days → no scored terms).
    dth_days = deaths_history.days
    dth_obs = isempty(deaths_history.counts) ? missing :
              collect(Int.(deaths_history.counts))
    incare_deaths ~ to_submodel(vintage_increments_model(
        [deaths_daily[clamp(Int(d), 1, n)] for d in dth_days], dth_obs, k))
    ro_days = ruleout_history.days
    ro_obs = isempty(ruleout_history.counts) ? missing :
             collect(Int.(ruleout_history.counts))
    ruleouts ~ to_submodel(vintage_increments_model(
        [ruleout_daily[clamp(Int(d), 1, n)] for d in ro_days], ro_obs, k))
    adm_h_days = admissions_history.days
    adm_h_obs = isempty(admissions_history.counts) ? missing :
                collect(Int.(admissions_history.counts))
    admissions ~ to_submodel(vintage_increments_model(
        [admit_daily[clamp(Int(d), 1, n)] for d in adm_h_days], adm_h_obs, k))
    ab_days = absconded_history.days
    ab_obs = isempty(absconded_history.counts) ? missing :
             collect(Int.(absconded_history.counts))
    absconded ~ to_submodel(vintage_increments_model(
        [abscond_daily[clamp(Int(d), 1, n)] for d in ab_days], ab_obs, k))

    ## Cut-off reported quantities.
    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]
    overall_los = CFR_iso * death_los_state.mean +
                  (one(CFR_iso) - CFR_iso) * recovery_los_state.mean
    expected_isolation := safe_rate(occ_T)
    expected_bed_demand := safe_rate(dem_T)
    ## Cut-off daily flows: the end-of-grid value of each modelled daily
    ## series, the one-week-ahead forecast base for admissions, in-care
    ## deaths and rule-outs.
    expected_admissions := safe_rate(isempty(admit_daily) ? z0 : admit_daily[end])
    expected_incare_deaths := safe_rate(isempty(deaths_daily) ? z0 :
                                        deaths_daily[end])
    expected_ruleouts := safe_rate(isempty(ruleout_daily) ? z0 :
                                   ruleout_daily[end])
    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
    incare_cfr := CFR_iso
    incare_cfr_modifier := β_iso
    treatment_overall_los := overall_los
    ## Cumulative occupancy reclassification offset at the cut-off (fitted; can
    ## be negative, so reported raw). Zero unless `occupancy_break_days` is set.
    occupancy_break := occ_break_val

    return (; p_iso, p_iso_bvd, δ_iso = sev_state.δ_iso,
        CFR_iso, β_iso, capacity = C_T,
        death_los_mean = death_los_state.mean,
        recovery_los_mean = recovery_los_state.mean,
        ruleout_los_mean = ruleout_los_state.mean,
        admission_delay_mean = adm_delay_state.mean,
        overall_los, abscond_frac, k_isolation = k,
        demand, occupancy, isolation, occupancy_break = occ_break_val,
        deaths_daily, recover_daily, ruleout_daily, admit_daily,
        expected_isolation = safe_rate(occ_T),
        expected_bed_demand = safe_rate(dem_T),
        expected_admissions = safe_rate(isempty(admit_daily) ? z0 :
                                        admit_daily[end]),
        expected_incare_deaths = safe_rate(isempty(deaths_daily) ? z0 :
                                           deaths_daily[end]),
        expected_ruleouts = safe_rate(isempty(ruleout_daily) ? z0 :
                                      ruleout_daily[end]))
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{(27)}& 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{(28)}& 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{(29)}& 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{(30)}& 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(),
        overdispersion = confirmed_overdispersion_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)

    ## Intra-window overdispersion for the confirmed positives. The observed
    ## and anchored-late windows are scored as an overdispersed BetaBinomial
    ## of the observed analysed denominator (`safe_betabinomial`), so the
    ## confirmed intervals are not the far-too-tight plain Binomial. Sampled
    ## once and shared across all confirmed windows regardless of the
    ## positivity link.
    od_state ~ to_submodel(overdispersion, false)
    ρ_conf = od_state.ρ

    ## 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 BetaBinomial 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: overdispersed BetaBinomial of the observed analysed
    ## denominator (`ρ_conf` the intra-window overdispersion).
    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,
        ρ_conf))

    ## Late windows: confirmed-only vintages after the last laboratory date.
    ## A day that publishes a 24h analysed count (`late_analysed > 0`) is
    ## scored as an overdispersed BetaBinomial 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, ρ_conf))

    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{(31)}& 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 of an independent dispersion $k_{\text{rec}}$:

$$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{(32)}& {\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{(33)}& 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{(34)}& {\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{(35)}& 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[]),
        treatment_admissions_history = (; days = Int[], counts = Int[]),
        treatment_deaths_history = (; days = Int[], counts = Int[]),
        treatment_ruleout_history = (; days = Int[], counts = Int[]),
        treatment_absconded_history = (; days = Int[], counts = Int[]),
        occupancy_break_days::AbstractVector{<:Integer} = 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_flow_model,
        recovered = recovered_model,
        dispersion = pooled_dispersion_model,
        ascertainment = pooled_ascertainment_model(),
        background_re::Bool = false,
        ## Observed daily contact-tracing follow-up rate `(; days, values)`.
        ## Passed to the suspected-case stream, where it drives a shared latent
        ## follow-up-rate process ([`contact_tracing_model`](@ref)) that anchors
        ## both the non-BVD background scaling and the BVD reporting effort.
        ## Empty (the default) leaves the anchor off.
        contact_followup_history = (; days = Int[], values = Float64[]),
        ## Give the suspected-case stream a case-finding reporting-effort
        ## multiplier ([`reporting_effort_walk_model`](@ref)) on the BVD
        ## ascertainment, anchored to the latent contact rate plus a walk.
        ## Default `false` keeps the constant-ascertainment suspected likelihood.
        suspected_reporting_effort::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

    ## Optional suspected-case reporting-effort walk, gated to the surveillance
    ## window (from the first suspected-case report). Injected into the cases
    ## submodel so a suspected-specific reporting slowdown is absorbed as
    ## changing case-finding effort rather than by the shared reproduction-number
    ## walk. `false` (default) passes `nothing`, leaving the suspected likelihood
    ## unchanged.
    ## The effort closure accepts the latent contact rate at call time (built
    ## inside the cases submodel), anchoring the BVD ascertainment to it plus a
    ## walk gated to the surveillance window.
    if suspected_reporting_effort
        eff_onset = isempty(reported_history.days) ? 1 :
                    clamp(Int(reported_history.days[1]), 1, n)
        case_effort_re = (nn; contact_covariate = nothing) -> reporting_effort_walk_model(
            nn; onset = eff_onset, contact_covariate)
    else
        case_effort_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,
        contact_followup_history, reporting_effort = case_effort_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_flow_model`](@ref)). The non-BVD rule-out stay is a
    ## separate parameter from the lab-turnaround `receipt_pmf`.
    ## Treatment-centre patient flow: occupancy plus the in-care outcome flows
    ## (admissions, in-care deaths, rule-outs, absconded), with the in-care
    ## fatality CFR_iso (a modifier on the infection CFR) identified by the
    ## in-care death flow. The Tableau 6 flow histories are optional refinements
    ## (empty → no-op).
    treatment_state ~ to_submodel(
        treatment(isolation_history, cases_state.bvd_reports_daily,
        cases_state.bg_daily, p_drc, deaths_state.CFR;
        capacity_history = bed_capacity_history,
        admissions_history = treatment_admissions_history,
        deaths_history = treatment_deaths_history,
        ruleout_history = treatment_ruleout_history,
        absconded_history = treatment_absconded_history,
        occupancy_break_days = occupancy_break_days,
        k_external = k_isolation))
    ## Recovered among confirmed ("cumul guéris"): survivors among the modelled
    ## daily confirmed cases (the confirmed-and-discharged subset, not all
    ## in-care recoveries), 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)
    ## Cut-off daily treatment flows surfaced for the one-week-ahead forecast.
    expected_admissions_T := treatment_state.expected_admissions
    expected_incare_deaths_T := treatment_state.expected_incare_deaths
    expected_ruleouts_T := treatment_state.expected_ruleouts
    bed_capacity := treatment_state.capacity
    isolation_admission := treatment_state.p_iso
    isolation_bvd_admission := treatment_state.p_iso_bvd
    isolation_severity := treatment_state.δ_iso
    ## BVD bed stay is now the outcome mixture; `isolation_bvd_los_mean`
    ## reports the mixture mean (overall length-of-stay), with the death and
    ## recovery branch means surfaced separately.
    isolation_bvd_los_mean := treatment_state.overall_los
    isolation_death_los_mean := treatment_state.death_los_mean
    isolation_recovery_los_mean := treatment_state.recovery_los_mean
    isolation_ruleout_los_mean := treatment_state.ruleout_los_mean
    isolation_admission_delay_mean := treatment_state.admission_delay_mean
    isolation_dispersion := treatment_state.k_isolation
    ## In-care fatality CFR_iso (a modifier on the infection CFR) and the
    ## abscond fraction.
    incare_cfr := treatment_state.CFR_iso
    incare_cfr_modifier := treatment_state.β_iso
    abscond_fraction := treatment_state.abscond_frac
    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
    ## Contact-tracing background elasticity (0 when no follow-up data).
    background_contact_coef := cases_state.β_contact
    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	297.0	901.0	1848.0	8631.0	23423.0	143768.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.

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	2.126	5.0	1024
2	exports	1.01	683.0	6
3	deaths (DRC)	1.005	763.0	4
4	cases (DRC)	1.005	443.0	31
5	confirmed (DRC)	1.015	136.0	114
6	confirmed deaths (DRC)	1.01	359.0	3
7	isolation (DRC)	1.008	401.0	12
8	frozen (1wk back)	1.011	436.0	58

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{(36)}& {\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. The helper below performs one frozen joint re-fit and is reused by the forecast validation and matched-in-time results.

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

# The frozen re-fits are defined in the fit registry (`docs/fits/registry.jl`) and
# loaded through the cache 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% 4961–5571, 60% 4755–5896, 90% 3922–7091 infections to date, reported and unreported.

• Against the 1561 laboratory-confirmed cases by the cut-off that is roughly 2.5–4.5× 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-03-05–2026-03-16, 60% 2026-02-24–2026-03-16, 90% 2026-02-11–2026-03-16, an elapsed age to the cut-off of 30% 110–121, 60% 110–130, 90% 110–143 days.

• Growth rate and doubling time: the initial growth rate is estimated to have been 30% 0.038–0.039, 60% 0.036–0.039, 90% 0.031–0.046 per day, an initial doubling time of 30% 17.9–18.1, 60% 17.6–19.2, 90% 15.2–22.5 days. The latest growth rate is estimated to be 30% 0.033–0.073, 60% 0.018–0.077, 90% -0.007–0.084 per day, a latest doubling time of 30% 9.3–15.4, 60% 8.7–26.2, 90% -49.3–77.6 days.

• Reproduction number: the initial reproduction number is estimated to have been 30% 1.15–1.42, 60% 1.14–1.64, 90% 1.14–2.0 and the latest to be 30% 1.26–1.34, 60% 1.22–1.41, 90% 0.93–1.87.

• Case-fatality ratio: the case-fatality ratio is estimated to be 30% 0.38–0.38, 60% 0.37–0.44, 90% 0.33–0.51.

• 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.1, the outbreak age by -0.67 and the doubling time by -0.18; the largest move is in the outbreak age.

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	3922.0	4755.0	4961.0	5571.0	5896.0	7091.0

The figure below shows the cumulative trajectories and current-cut-off densities for three latent quantities: infections, symptom onsets and deaths. The infection density is the headline outbreak size, a count of infections rather than reported cases.

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. 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.01	0.02	0.03	0.07	0.08	0.08
2	doubling_time	-49.32	8.71	9.29	15.42	26.17	77.58
3	T	109.58	109.76	110.14	121.37	129.56	142.58
4	R_T	0.93	1.22	1.26	1.34	1.41	1.87
5	CFR	0.33	0.37	0.38	0.38	0.44	0.51
6	C_T	3922.04	4755.21	4961.42	5571.13	5895.95	7090.51

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	1.66	2.33	2.36	2.51	2.74	3.51
2	generation interval scale	1.42	1.44	1.47	4.69	5.81	7.3
3	incubation period mean	5.66	6.18	6.62	6.72	6.79	7.05
4	incubation period SD	2.33	2.37	2.4	3.05	3.64	4.48

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_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 reproduction number is identified alongside the suspected-case case-finding, not in place of it. Both the non-BVD suspected-case background and the suspected-case ascertainment are grounded in the observed contact-tracing follow-up rate (the "taux de suivi des contacts"), a shared latent case-finding-intensity process that scales detection of non-BVD and BVD suspects alike. A change in case-finding reporting is therefore carried by that channel rather than by the reproduction number, so the reproduction number tracks transmission rather than reporting. The reporting-effort multiplier on the BVD ascertainment carries a random-walk deviation on top of the contact anchor, and the anchor is defined over the whole grid so it continues under the walk when the follow-up data stop (see the contact-tracing case-finding sensitivity on the sensitivity page).

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.61	0.63	0.64	0.76	0.88	0.97

Observation delays

The delays carry latent infections through to each observed event: reporting, death, detection abroad and laboratory receipt. 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.81	1.1	1.34	1.38	1.41	1.56
2	onset-to-report scale	2.27	2.58	2.69	3.14	3.98	5.18
3	onset-to-admission shape	0.79	1.08	1.2	1.21	1.23	1.5
4	onset-to-admission scale	2.08	3.25	3.7	3.8	3.87	5.26
5	admission-to-death shape	1.03	1.05	1.06	1.79	2.22	2.79
6	admission-to-death scale	2.02	3.4	4.41	4.88	5.05	5.45
7	onset-to-detection shape	0.85	0.91	0.92	1.05	1.26	1.54
8	onset-to-detection scale	2.29	3.48	4.38	5.34	5.45	5.58
9	report-to-receipt mean	2.97	3.33	3.37	3.67	4.41	5.59
10	report-to-receipt SD	3.1	3.87	4.43	5.18	5.23	5.28
11	in-care BVD length-of-stay mean (mixture)	7.79	10.31	11.12	11.32	11.67	13.92
12	isolation non-BVD rule-out stay mean	1.64	1.71	1.76	2.82	4.23	6.08
13	confirmation-to-recovery mean	15.78	20.08	21.28	21.62	21.92	25.2

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,
        :background_contact_coef,
        Symbol("cases_state.effort_state.β_asc"),
        Symbol("cases_state.effort_state.σ_eff"),
        :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,
        :incare_cfr, :incare_cfr_modifier, :isolation_death_los_mean,
        :isolation_recovery_los_mean, :abscond_fraction,
        :recovery_probability, :recovered_dispersion, :expected_recovered_T];
    digits = 3);

Show surveillance-parameter summary table

37×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.634	0.794	0.866	0.879	0.881	0.943
2	p_uganda	0.627	0.79	0.878	0.89	0.892	0.943
3	k	2.382	2.995	3.021	3.448	5.212	8.611
4	k_cases	4.343	5.709	7.055	9.71	9.82	9.947
5	k_deaths	6.671	6.949	7.009	9.99	13.474	20.058
6	k_confirmed	1.965	2.511	3.059	4.107	4.134	4.214
7	k_confirmed_deaths	1.114	1.158	1.169	1.418	1.78	2.448
8	dispersion_sd	0.649	0.833	0.931	0.975	0.979	1.209
9	tau_test	0.833	0.854	0.863	0.885	0.927	0.967
10	lambda_bg	13.374	14.341	15.053	17.03	21.212	27.218
11	background_contact_coef	0.529	1.734	2.914	3.64	3.719	4.019
12	cases_state.effort_state.β_asc	0.148	0.328	0.336	0.43	0.942	2.084
13	cases_state.effort_state.σ_eff	0.008	0.008	0.009	0.044	0.108	0.22
14	suspected_positivity	0.295	0.352	0.38	0.408	0.418	0.465
15	test_positivity	0.251	0.26	0.27	0.283	0.295	0.32
16	expected_confirmed_T	1530.61	1590.56	1645.13	1722.86	1806.85	1955.93
17	expected_analysed_T	5539.23	5834.9	5930.72	6105.54	6264.69	6684.24
18	death_ascertainment	0.809	0.853	0.858	0.864	0.904	0.943
19	background_cfr	0.187	0.262	0.327	0.358	0.359	0.397
20	tau_death	0.532	0.654	0.689	0.703	0.729	0.929
21	death_composition	0.394	0.416	0.443	0.485	0.561	0.691
22	death_confirmation	0.297	0.31	0.339	0.388	0.465	0.565
23	expected_confirmed_deaths_T	415.773	429.42	450.538	484.321	528.508	606.251
24	isolation_admission	0.252	0.295	0.309	0.317	0.319	0.368
25	isolation_dispersion	138.994	229.045	270.304	280.502	298.396	520.386
26	expected_isolation_T	611.189	630.533	645.009	671.736	689.019	713.552
27	expected_bed_demand_T	617.783	636.756	652.324	685.718	704.401	745.646
28	bed_capacity	627.922	673.74	689.684	699.749	708.561	751.148
29	bed_shortfall_T	0.0	0.0	0.0	0.0	14.93	76.091
30	incare_cfr	0.368	0.409	0.419	0.425	0.464	0.533
31	incare_cfr_modifier	-0.298	0.023	0.147	0.191	0.228	0.531
32	isolation_death_los_mean	3.171	6.014	7.552	8.193	8.294	11.61
33	isolation_recovery_los_mean	7.85	12.626	13.228	14.136	14.853	19.162
34	abscond_fraction	0.004	0.005	0.005	0.005	0.006	0.006
35	recovery_probability	0.315	0.392	0.463	0.481	0.499	0.537
36	recovered_dispersion	1.947	2.969	4.194	5.175	5.335	6.22
37	expected_recovered_T	248.682	277.484	289.392	303.023	312.43	346.698

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 cover two groups: the dated DRC surveillance streams and the Uganda exports. The latent infection process is not checked here, as it carries no direct observation, and is shown instead as the estimated cumulative trajectories in the joint model estimates figure.

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),
        # Kept so the generator's occupancy-break dimension matches the fitted
        # chain (the offset step on the `[occupancy_break_dates]` days).
        occupancy_break_days = obs.occupancy_break_days,
        recovered_history = _days_only(obs.recovered_history),
        treatment_admissions_history =
        _days_only(obs.treatment_admissions_history),
        treatment_deaths_history = _days_only(obs.treatment_deaths_history),
        treatment_ruleout_history = _days_only(obs.treatment_ruleout_history),
        treatment_absconded_history =
        _days_only(obs.treatment_absconded_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);

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

# 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.
vintage_panels = [
    reported_panel, suspected_daily_panel, isolation_panel, confirmed_panel,
    deaths_panel, suspected_daily_deaths_panel, confirmed_deaths_panel,
    recovered_panel, tests_analysed_panel, tests_analysed_daily_panel,
    admissions_panel, incare_deaths_panel, ruleouts_panel, absconded_panel];
joint_vintage_ppc_fig = plot_vintage_conditional_ppc(vintage_panels);

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

We score each stream's per-vintage conditional predictions against the observed counts. bias is the mean forecast bias over the vintages (negative = under-predicted, positive = over-predicted, zero = the observed counts sit at the predictive median); 50%/90% coverage are the fractions of vintages whose observed count falls inside the central 50% and 90% predictive intervals, which a well-calibrated stream keeps near those nominal levels. Streams with a large bias or coverage far from nominal are the ones the joint fit reproduces less well.

stream_calibration_table = stream_calibration(vintage_panels);

The calibration plot reads the table at a glance: the left panel marks each stream's empirical 50% and 90% coverage against dashed reference lines at the nominal levels, and the right panel marks the mean forecast bias against a dashed line at zero.

Per-stream calibration plot

stream_calibration_fig = plot_stream_calibration(stream_calibration_table);

Per-stream calibration table

14×5 DataFrame

Row	Stream	Vintages	Bias	50% coverage	90% coverage
	String	Int64	Float64	Float64	Float64
1	Suspected cases	9	-0.03	0.56	0.78
2	New suspects/day	28	-0.21	0.61	1.0
3	Patients in isolation	28	-0.18	0.64	0.89
4	Confirmed cases	45	0.06	0.53	0.93
5	Suspected deaths	9	-0.02	0.56	0.78
6	New suspected deaths/day	24	-0.11	0.62	1.0
7	Confirmed deaths	47	-0.02	0.68	0.98
8	Recovered (confirmed)	25	0.06	0.6	0.92
9	Specimens analysed (cumulative)	6	-0.23	0.0	0.67
10	Specimens analysed (24h)	28	-0.3	0.57	0.93
11	Admissions/day	16	-0.3	0.81	0.94
12	In-care deaths/day	16	0.04	0.5	0.81
13	Rule-outs/day	15	-0.19	0.6	0.8
14	Absconded/day	16	0.09	0.5	0.88

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

Posterior correlations and stream totals

The heatmap is the posterior correlation between each pair of headline quantities: the outbreak size ($C_T$), the reproduction number ($R_T$), the outbreak age ($T$), the case-fatality ratio (CFR), the DRC and Uganda ascertainment fractions ($p_{\text{drc}}$, $p_{\text{ug}}$), the non-BVD background rate (${\lambda }_{\text{bg}}$), the fraction tested (${\tau }_{\text{test}}$), and the cut-off total expected for each stream. Blue is positive, red negative.

Posterior correlation heatmap

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

The stream-total plot takes each posterior draw, sums every stream over its own reporting dates, and marks the observed total with a crosshair. The diagonal panels are the predictive spread of each total against the observed value; the off-diagonal panels show whether the totals move together from draw to draw.

Stream totals against observed

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

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	499.0	1144.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;

4×3 DataFrame

Row	Quantity	Central estimate	Narrowest interval
	String	String	String
1	Delay-corrected confirmed CFR	35.2%	31.7–50.4%
2	Structural (infection-based) CFR	38.0%	33.1–51.2%
3	Uncorrected modelled confirmed ratio	27.3%	24.2–35.3%
4	Naive observed confirmed ratio	32.4%	—

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.

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

One-week-ahead forecast summary table

11×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	1716.0	1864.0	1961.0	2196.0	2370.0	2763.0
2	DRC confirmed cases	new this week	155.0	303.0	400.0	635.0	809.0	1202.0
3	DRC confirmed deaths	cumulative by T+7	299.0	486.0	589.0	749.0	863.0	1107.0
4	DRC confirmed deaths	new this week	0.0	0.0	83.0	243.0	357.0	601.0
5	DRC isolation beds	demand at T+7	590.0	769.0	916.0	1292.0	1401.0	1574.0
6	DRC isolation beds	occupancy at T+7	578.0	656.0	688.0	698.0	706.0	734.0
7	DRC isolation admissions	daily at T+7	77.0	106.0	131.0	207.0	226.0	253.0
8	DRC in-care deaths	daily at T+7	13.0	20.0	25.0	35.0	40.0	49.0
9	DRC isolation rule-outs	daily at T+7	37.0	49.0	60.0	102.0	111.0	125.0
10	DRC recovered	cumulative by T+7	245.0	269.0	308.0	421.0	496.0	612.0
11	DRC recovered	new this week	6.0	30.0	69.0	182.0	257.0	373.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);

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

The flow figure projects the daily isolation/treatment flows a week ahead: new admissions, in-care deaths and rule-outs, each grown from its cut-off daily rate and replicated through the isolation dispersion. These are the daily-flow counterparts of the bed-stock forecast above.

One-week-ahead treatment-flow forecast plot

forecast_flows_fig = plot_forecast_flows(forecast);

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)

# 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, "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.