BVDOutbreakSize

Appearance

Estimating the current size of the 2026 DRC Bundibugyo virus outbreak ​

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

Stable Dev Tests codecov Aqua QA Code Style: SciML DOI

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

Data as of: 1 August 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  , where 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. is held flat at the established 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 . 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
julia
# Shared setup: packages, observations, the fit registry and every model fit
# (loaded from the content-addressed cache). See `docs/examples/_setup.jl`.
using BVDOutbreakSize
include(joinpath(pkgdir(BVDOutbreakSize), "docs", "examples", "_setup.jl"))
45

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. The same table breaks the occupancy into dont confirmés (NC+AC) and dont suspects sub-rows, two prevalence sub-stocks that sum to the total each day, which we read as two further census streams splitting the occupancy. 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
julia
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
Rowfielddatevalue
StringDate?Int64
1exported_cases2026-05-233
2exports_deaths2026-05-141
3suspected_deaths2026-05-26246
4suspected_cases2026-05-261077
5confirmed_cases2026-08-013748
6confirmed_deaths2026-08-011657
7specimens_analysed2026-05-28755
8treatment_admissions2026-08-01122
9treatment_deaths2026-08-0125
10treatment_ruleouts2026-08-0186
11treatment_absconded2026-08-014
12genetic_tmrca_bound2026-03-15139
13daily_outbound_travellers (prior mean)missing1871
14daily_outbound_travellers_sd (prior SD)missing200
15source_populationmissing4392200

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
julia
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
72×10 DataFrame
Rowdatesuspected_casessuspected_new_dailypatients_isolatedsuspected_deathsconfirmed_casesconfirmed_deathsrecovered_confirmedspecimens_receivedspecimens_analysed
DateInt64?Int64?Int64?Int64?Int64Int64Int64?Int64?Int64?
12026-05-14missingmissingmissingmissing84missingmissingmissing
22026-05-17missingmissingmissingmissing134missingmissingmissing
32026-05-18516missingmissing131334missingmissingmissing
42026-05-19575missingmissing148514missingmissingmissing
52026-05-20672missingmissing160646missingmissingmissing
62026-05-21745missingmissing175839missingmissingmissing
72026-05-22872missingmissing2049110missingmissingmissing
82026-05-23904missingmissing22010110missing418211
92026-05-24906missingmissing22310510missing431295
102026-05-25998missingmissing23810612missing431295
112026-05-261077missingmissing24612117missing662403
122026-05-27missingmissingmissingmissing12517missing774648
132026-05-28missingmissingmissingmissing21017missing883755
142026-05-29missingmissingmissingmissing26342missingmissingmissing
152026-05-30missingmissingmissingmissing28242missingmissingmissing
162026-05-31missingmissingmissingmissing32148missingmissingmissing
172026-06-01missingmissing173missing34460missingmissingmissing
182026-06-02missingmissing206missing36362missingmissingmissing
192026-06-03missingmissing233missing38164missingmissingmissing
202026-06-04missing153258missing45282missingmissingmissing
212026-06-05missing119267missing48886missingmissingmissing
222026-06-06missing117283missing5159112missingmissing
232026-06-07missing94309missing55010119missingmissing
242026-06-08missing138297missing59811522missingmissing
252026-06-09missing119260missing63512730missingmissing
262026-06-10missing119262missing67613632missingmissing
272026-06-11missing168315missing68913932missingmissing
282026-06-13missing136359missing78218140missingmissing
292026-06-14missing165363missing80819248missingmissing
302026-06-15missing235376missing83719649missingmissing
312026-06-16missing192379missing87520267missingmissing
322026-06-17missing151383missing89623278missingmissing
332026-06-18missing238416missing93324580missingmissing
342026-06-19missing162361missing95624792missingmissing
352026-06-20missing201365missing1003254100missingmissing
362026-06-21missing202371missing1048267112missingmissing
372026-06-22missing131387missing1094277115missingmissing
382026-06-23missing138408missing1118291122missingmissing
392026-06-24missing154385missing1155304138missingmissing
402026-06-25missing265419missing1203321148missingmissing
412026-06-27missing239502missing1274360178missingmissing
422026-06-29missing309609missing1333399189missingmissing
432026-06-30missing301missingmissing1406438208missingmissing
442026-07-01missing150641missing1460452213missingmissing
452026-07-02missing213628missing1502473229missingmissing
462026-07-03missing185missingmissing1528492239missingmissing
472026-07-04missing354missingmissing1561506254missingmissing
482026-07-05missing135646missing1624521273missingmissing
492026-07-06missing237680missing1708580280missingmissing
502026-07-07missing304750missing1759600285missingmissing
512026-07-08missing227764missing1792625295missingmissing
522026-07-09missing284780missing1830648300missingmissing
532026-07-10missing299763missing1873672306missingmissing
542026-07-11missing299753missing1926702318missingmissing
552026-07-12missing268736missing1963719333missingmissing
562026-07-13missing268753missing2011754366missingmissing
572026-07-14missing418736missing2073796377missingmissing
582026-07-15missing389725missing2124828390missingmissing
592026-07-17missing236722missing2267893412missingmissing
602026-07-18missing192724missing2344930466missingmissing
612026-07-19missing252734missing2423967469missingmissing
622026-07-20missing322737missing2473999482missingmissing
632026-07-21missing306738missing25361033506missingmissing
642026-07-22missing318722missing29051269519missingmissing
652026-07-23missing315766missing29731309540missingmissing
662026-07-24missing274755missing30751354556missingmissing
672026-07-25missing340773missing32001405571missingmissing
682026-07-26missing326723missing32621437583missingmissing
692026-07-27missing321733missing33601487597missingmissing
702026-07-30missing374748missing36051587651missingmissing
712026-07-31missing321760missing36741621666missingmissing
722026-08-01missing227690missing37481657708missingmissing

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 ). The generating infection process produces daily infection incidence via the discrete renewal equation

where 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 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:

ParameterExportsDeathsCasesAnalysedConfirmedConf. deathsExport deaths
Reproduction number ●●●●●●●
Generation interval●●●●●●●
Incubation period●●●●●●●
Seed ●●●●●●●
Onset-to-death delay●●●
Case-fatality ratio●●●
Death ascertainment ●●
Background CFR ●●
Onset-to-report delay●●●
Receipt delay●●●
Onset-to-detection delay●
Assay sensitivity / specificity●●
Severity enrichment ●
Death testing fraction ●
Testing fraction ●●
Background rate ●●●●●
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 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 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 at its first knot:

We do not place a prior on directly. We put the prior on the initial growth rate 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 :

The step-size prior keeps weekly changes in the reproduction number moderate. We set the half-normal on 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 .

Daily 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 (the interpolation clamps below the first knot rather than extrapolating):

with the day of knot . 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:

Submodel: rt_walk_model
julia
@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 = RT_INTERVENTION_RAMP,
        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 is a Gamma distribution with a sampled shape and scale , 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 and scale near , 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:

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
julia
@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 on the last day of the seeding window has a prior centred on a single seed:

From that seed we assume the outbreak grew deterministically through an unobserved cryptic exponential phase, doubling times before sustained transmission was established. The cryptic phase grows the seed to infections at the renewal start, the day the renewal takes over, over a duration   with the doubling time. Field epidemiology in Mongbwalu traced a sustained transmission chain back to a death on 25 January 2026 and identified more than 500 suspected cases between mid-January and mid-May (Kupferschmidt, 2026); the genetic TMRCA (Mbala-Kingebeni and others, 2026) is a lower bound on the outbreak age that is consistent with, but does not by itself fix, an origin that early. We place a wide prior on centred so that the implied origin sits at the end of January, at the documented first deaths:

The growth rate carries the prior the genetic source informs. The BEAST X reanalysis of 139 BDBV genomes (Mbala-Kingebeni and others, 2026) reports an Exponential-growth doubling time of 11.7 d (95% HPD 6.8–17.5). We put a log-normal prior on equivalent to a log-normal prior on the doubling time centred on 11.7 d, with its log spread read from that HPD and inflated a little, so the prior is slightly wider in spread than the source but unbiased relative to it:

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 to uses its own generation interval; deriving forward from the shared growth rate under our generation interval is the consistent choice.

Submodel: seed_model
julia
@model function seed_model(; i0_prior = truncated(Normal(0.1, 0.1); lower = 0))
    I0 ~ i0_prior
    return (; I0)
end
Submodel: exponential_growth_model
julia
@model function exponential_growth_model(;
        r_prior = LogNormal(log(log(2) / M_PRIOR_DOUBLING_DAYS), 0.28),
        m_prior = truncated(Normal(5.0, 4.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   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 (the cryptic duration plus the observed window, defined in the infection process below):

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
julia
@model function genetic_seeding_model(T::Real,
        tmrca_days::Union{Missing, Real}; tmrca_days_sd::Real = 16.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

The grid days before the renewal start are filled by the cryptic exponential curve at rate ending at , 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:

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 divided by that rate.

Submodel: infection_model
julia
@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,  ; (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:

Every delay is discretised to a daily PMF over lags 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 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

which is then renormalised over lags .

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
julia
@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
julia
@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:

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 -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:

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:

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:

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 deaths in cases (; CDC outbreak history), with confidence bands spanning roughly -. The companion Bundibugyo virus (BDBV) reanalysis reports a baseline of ( CrI -) for non-healthcare-worker (non-HCW) confirmed cases. Based on this we use a prior of

with mean and interval roughly -. The mean matches the CDC   figure and the corrected central CFR in the 20 May report (McCabe and others, May 2026).

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

The prior density, with the CDC 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 scale in non-centred log form:

so    per stream, with setting the pooling (  collapses to one shared dispersion). The population value   is the headline dispersion.

Submodel: pooled_dispersion_model
julia
@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 and share a logit-scale hyperprior with mean and pooling strength , centred on a reporting fraction of , reflecting the active case-finding of a declared Ebola response rather than baseline passive surveillance:

The cases likelihood uses ; the two Uganda-side likelihoods use .

Submodel: pooled_ascertainment_model
julia
@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 is the share of suspected cases routed to the laboratory. A truly BVD specimen tests positive with the assay sensitivity , and a non-BVD specimen tests positive with the false-positive rate   from the assay specificity. We assume that more severe cases, more likely to be Ebola, are preferentially tested, through an enrichment factor 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 of suspected deaths reach the laboratory, and they confirm at the assay positivity       built from the same sensitivity and specificity but the death-pool BVD share . 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:

The non-BVD background rate enters the suspected-case stream and is described with it below; the suspected deaths carry a death ascertainment    and a non-BVD death background tied to the case background by a background CFR  .

Submodel: test_positivity_model
julia
@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
julia
@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
julia
@model function test_sensitivity_model(;
        sensitivity_prior = Beta(38.0, 2.0))
    s_test ~ sensitivity_prior
    return (; s_test)
end
Submodel: test_specificity_model
julia
@model function test_specificity_model(; specificity_prior = Beta(60.0, 2.0))
    spec ~ specificity_prior
    return (; spec)
end
Submodel: severity_enrichment_model
julia
@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
julia
@model function death_testing_fraction_model(; fraction_prior = Beta(5.0, 2.0))
    τ_death ~ fraction_prior
    return (; τ_death)
end
Submodel: death_ascertainment_model
julia
@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
julia
@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 is a sample mean across roughly - point-of-entry-weeks. We use a Normal prior centred on with SD ( CV), truncated at zero, covering point-of-entry variation and the sitrep sampling uncertainty; the source population is kept fixed (census):

Submodel: traveller_volume_model
julia
@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 and scaled by the DRC ascertainment . The convolution of a daily series with a delay PMF is the lagged sum

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

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 that follows a lognormal random walk on weekly knots around a baseline , linearly interpolated to the daily grid,

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

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

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" on report day , 4-7 June). This is a genuine daily incidence, not a cumulative total, so it is scored against the modelled daily suspected count on that day directly (a single-day mean, not a between-vintage sum) with a NegBinomial sharing :

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 (see deaths_model).

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

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

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

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

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

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

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

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

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

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

    return (; p_drc, λ_bg = λ_bg_base, τ_test, report_pmf,
        report_mean = report_state.mean, report_sd = report_state.sd,
        bvd_reports_daily,
        reports_daily, expected_reports, positivity, bg_daily, bg_sigma,
        bg_total)
end
Treatment-centre flow ​

The treatment-centre stream models the daily patient flow through the isolation/treatment centres: the occupied-bed count ("Patients en isolement"), the daily admissions, and the daily discharges split by outcome — in-care deaths, rule-outs and absconded — read from the situation-report Tableau 6 patient-movement table. Two parallel processes act on each patient. A clinical course governs how long a patient occupies a bed and how they leave it, and so sets the total occupancy and every discharge flow; a laboratory label runs alongside it and only relabels a patient from suspected to confirmed, carving the suspect/confirmed split of the census. Death is clinical and happens under either label, so a true case may die before its test confirms it. Separating the two keeps the operational churn in the suspected pool out of the part of the occupancy the infection estimate leans on.

A proportion of the reported suspects need a bed. These admissions split into a BVD true-case inflow, admitted at a severity-skewed rate     above the base rate, and a non-BVD background inflow at the base rate,

each carried through a short suspected-to-admission delay that captures triage, transport and the wait for a bed. A patient then leaves by one of four routes. A BVD true case dies at the in-care case-fatality ratio over the admission-to-death stay or recovers over the longer admission-to-recovery stay; a non-case is ruled out by a negative test over the rule-out stay or absconds. 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 takes the report-to-receipt laboratory turnaround.

Occupancy is the running balance of these latent events rather than a length-of-stay convolution: each day the bed stock is yesterday's stock plus the day's admissions less the day's deaths, recoveries, rule-outs and absconds. The abscond outflow drains the suspected pool at a small daily fraction of the previous day's suspected occupancy,

so the total bed demand is the forward balance

with   . The stay lives entirely in the discharge flows, each an admission stream convolved with its outcome density, so one inflow and one set of outcome timings generate the bed stock, the discharge flows and the demand together. The occupancy accumulates admissions over the stay, so it is a smooth integral of the infection signal: a high, sustained bed stock informs the reproduction number, and anything in the occupancy that is not infection, such as an overnight reclassification of who is counted, is modelled rather than left to bend the transmission estimate.

In-care deaths combine the two labels: a true case who dies before its test returns is a suspected death and one who dies after is a confirmed death, and the report records the two together, so the death flow is the in-care fatality applied to the BVD inflow over the admission-to-death stay, scored against the combined deaths directly and never gated by confirmation. The in-care fatality is a sampled log-odds modifier on the infection case-fatality ratio,

identified by the in-care death flow relative to admissions and occupancy. It is a fatality conditional on admission rather than a causal treatment effect, sitting below the infection case-fatality ratio where care lowers mortality, and it is reported with and the overall length of stay (the death/recovery mixture mean).

The laboratory label carves the census into a confirmed and a suspected sub-stock. Confirmation relabels a true case already in a bed at the daily hazard   : the community confirmation hazard   — the share of suspects routed to the laboratory times the day's positivity — borrowed from the confirmed-case pipeline rather than re-estimated, so the in-care confirmed stock is a subset of the total confirmed by construction. The confirmed-in-care stock is tracked by admission cohort. Each true-case admission carries two clocks from the day it enters a bed, a confirmation clock and a clinical-stay clock, and it counts toward the confirmed census only once it has been confirmed and while it is still in a bed,

with the cumulative confirmation probability of a cohort admitted on day ,

a cumulative product along the cohort's age rather than a fixed distribution because the hazard is time-varying, and the clinical-stay survival,

the probability an admitted case is still in a bed after days, the discharge-complement of the death/recovery mixture built from the same admission-to-death and admission-to-recovery stays the discharge flows use. Because    reconstructs the occupied true-case stock, the confirmed-and-present cohort is a subset of it and    holds by construction.

Cohort tracking is needed because deaths and recoveries are observed combined across the two labels, so the data do not say which departing patients had already been confirmed. Carrying the confirmation and stay clocks separately excludes cases that die before their test returns from the confirmed pool, which a single pool-average discharge rate would over-attribute. The suspected sub-stock is the remainder   , holding the not-yet-confirmed BVD occupancy together with the non-case occupancy awaiting rule-out, and the abscond outflow drains this suspected stock at the daily fraction . Recoveries among the confirmed (the published recovery total) are the confirmed subset of recoveries and are modelled as a separate confirmed-recovery stream (below).

Capacity enters only as a censored observation; the latent demand is never capped, because the demand is the quantity of interest. The bed capacity 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 divided by the reported occupancy rate (about rising to beds over 9–13 June). The occupied beds are scored as the latent demand right-censored at the recorded implied capacity, so demand above a saturated capacity is left uncensored, and the daily admissions are right-censored at the recorded free-bed headroom, the implied capacity less the previous day's observed occupancy. Both censoring bounds are fixed recorded data, which keeps the admissions censor stable where a bound tied to the modelled, wandering capacity would drift. The occupancy and a flow stream are scored as

with each flow mean the matching modelled event series — the admissions, the in-care deaths, the rule-outs and the absconds — all sharing the treatment dispersion , and the implied capacity carried by a NegBinomial of its own. Occupancy below capacity identifies the demand directly; demand above a saturated capacity is only partially identified, since the occupancy reveals that demand was at least the beds filled but not how much more, so the bed shortfall above capacity is informed by the demand model and its priors rather than measured. Bed demand is the uncapped diagnostic, and 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. Out of sample, where no occupancy is observed, the forecast caps admissions at the modelled free beds; this modelled bound is safe because the forecast is a forward simulation outside the likelihood.

The fitted occupancy series is the all-patients column from 1 June (SitRep 018) onward, and from 13 June the report adds a two-row breakdown into confirmed and suspected beds that sums to the total each day. The total-occupancy term is the backbone present from 1 June; the daily flows and the confirmed/suspected census add likelihood on the days they exist, scored per day as either the total or the split so the total and its parts are never both counted on one day. The early window, with only the total occupancy reported, fits the backbone alone while the latent admissions still drive the stock, and the split is scored only where the borrowed confirmation hazard is non-zero, that is, where the laboratory pipeline of the full model supplies it.

One reporting artefact is modelled, on identified days only: an overnight reclassification of the total, where the published start-of-day in-bed count is differenced against the previous report day's occupancy, and a day whose gap exceeds a threshold is flagged as a break day. One step is fitted per flagged day, with a prior centred on that day's observed gap but free to move, so the fit can attribute part of a gap to genuine change in demand; the steps accumulate into a persistent additive offset on the modelled total occupancy, carried forward to every later day, that absorbs the overnight gap without bending the reproduction number to chase it. The split does not change the occupancy before 13 June, since no breakdown is published there and the total backbone carries that window.

Submodel: treatment_flow_model
julia
@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[]),
        ## Tableau 6 occupancy split (`dont confirmés` / `dont suspects`): two
        ## census series scored in place of the total occupancy on the days they
        ## are present, only when the confirmation hazard is non-zero. Empty by
        ## default, leaving the total-occupancy likelihood alone.
        confirmed_incare_history = (; days = Int[], counts = Int[]),
        suspect_incare_history = (; days = Int[], counts = Int[]),
        ## Daily in-care confirmation hazard `τ_test · p_pos` borrowed from the
        ## lab pipeline ([`confirmed_cases_model`](@ref)). `nothing` (standalone,
        ## no lab stream) gives a zero hazard, so the confirmed sub-stock stays
        ## empty and the suspect sub-stock carries the whole occupancy.
        conf_hazard_daily::Union{Nothing, AbstractVector} = nothing,
        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),
        ## In-care confirmation-rate modifier prior (log scale): γ_conf scales
        ## the borrowed community hazard to the effective in-care rate
        ## ρ = exp(γ_conf). Centred on zero (ρ = 1) so the census sets the split.
        incare_confirm_log_prior = Normal(0.0, 0.5),
        ## 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)),
        ## Opt-in occupancy reclassification-break days (grid indices). A level
        ## step is fitted into the modelled total at each, absorbing a
        ## measurement-basis discontinuity in the isolation series. Empty (the
        ## default) is a no-op. See `cumulative_occupancy_offset`.
        occupancy_break_days::AbstractVector{<:Integer} = Int[],
        ## Prior sd of each occupancy break step (beds), centred on zero.
        occupancy_break_sd::Real = 25.0)
    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 (loss-to-follow-up) drains the suspect pool at a small daily
    ## fraction κ of the previous-day suspect occupancy.
    abscond_frac ~ abscond_prior
    κ = abscond_frac

    ## Opt-in cumulative occupancy reclassification-break offset Δ(t) added to
    ## the modelled total occupancy, absorbing a between-report measurement-
    ## basis discontinuity in the observed isolation series (the `au-lit-J-1`
    ## versus `Fin-J` reclassification) on the manually supplied
    ## `occupancy_break_days`. A level step is fitted at each break day, sampled
    ## non-centred and centred on zero, so the fit partitions it into reporting
    ## artifact vs real demand. Applied additively to the modelled total only;
    ## demand (the diagnostic) stays the un-offset latent stock. Empty (the
    ## default) → no sampled step, Δ = 0, a no-op. Only break days on or before
    ## an observed occupancy day can move the likelihood; later ones are dropped
    ## so no inert step is sampled. See `cumulative_occupancy_offset`.
    iso_last = isempty(isolation_history.days) ? 0 :
               maximum(Int.(isolation_history.days))
    brk_days = [Int(d) for d in occupancy_break_days if Int(d) <= iso_last]
    if isempty(brk_days)
        b = Float64[]
    else
        occupancy_step ~ product_distribution(
            fill(Normal(0, 1), length(brk_days)))
        b = occupancy_break_sd .* occupancy_step
    end
    ## Cumulative reclassification offset Δ over the 1:n grid.
    break_grid_days = brk_days
    occ_break_offset = cumulative_occupancy_offset(1:n, break_grid_days, b)

    ## Admission inflow through the suspected→admission delay, split into BVD
    ## true-case (`p_iso_bvd`) and non-BVD (`p_iso`) inflows. Uncapped latent
    ## demand; capacity enters only as a censoring bound below.
    A_bvd = convolve_delay(p_iso_bvd .* p_drc .* bvd_reports_daily,
        adm_delay_state.pmf)
    A_bg = convolve_delay(p_iso .* bg_daily, adm_delay_state.pmf)
    if eltype(A_bvd) === Any
        A_bvd = convert(Vector{eltype(C)}, A_bvd)
        A_bg = convert(Vector{eltype(C)}, A_bg)
    end

    ## Label-independent clinical discharge events. Deaths and recoveries split
    ## `A_bvd` by `CFR_iso`; rule-outs discharge `A_bg`.
    dpmf = death_los_state.pmf
    rpmf = recovery_los_state.pmf
    deaths_daily = convolve_delay(CFR_iso .* A_bvd, dpmf)
    recover_daily = convolve_delay((one(CFR_iso) - CFR_iso) .* A_bvd, rpmf)
    ruleout_daily = convolve_delay(A_bg, ruleout_los_state.pmf)
    admit_daily = A_bvd .+ A_bg

    ## Community confirmation hazard `τ_test · p_pos` borrowed from the lab
    ## pipeline; `nothing` (standalone) gives a zero hazard.
    borrowed_hazard = if conf_hazard_daily === nothing
        zeros(eltype(A_bvd), n)
    elseif eltype(conf_hazard_daily) === Any
        convert(Vector{eltype(C)}, conf_hazard_daily)
    else
        conf_hazard_daily
    end

    ## The split census is identified only when the borrowed hazard is non-zero
    ## (a lab stream supplies it). With a structural zero the confirmed sub-stock
    ## is empty, so the split likelihood no-ops and those days stay on the total.
    split_active = any(>(zero(eltype(borrowed_hazard))), borrowed_hazard)

    ## In-care confirmation-rate modifier ρ = exp(γ_conf) on the borrowed hazard.
    ## Sampled only when the hazard is non-zero, so no unidentified dimension is
    ## added when the split is absent (matching the `k_external` pattern).
    if split_active
        incare_confirm_log ~ incare_confirm_log_prior
    else
        incare_confirm_log = zero(eltype(borrowed_hazard))
    end
    ρ_conf = exp(incare_confirm_log)
    conf_hazard = ρ_conf .* borrowed_hazard

    ## Forward running-balance occupancy: total demand and the BVD/non-case
    ## stocks. The scored abscond flow is recomputed below off the two-clock
    ## suspect stock.
    acc = accumulate_occupancy(A_bvd, A_bg, deaths_daily, recover_daily,
        ruleout_daily, κ, conf_hazard)
    demand = acc.demand
    O_bvd = acc.O_bvd

    ## Two-clock confirmed-in-care sub-stock: cohort-tracked confirmed-and-present
    ## prevalence, exact in the fast-death tail where the running balance's
    ## proportional drain is only mean-field. Demand and `O_bvd` stay as
    ## `accumulate_occupancy` built them; `O_conf` (and `O_susp = D − O_conf`) is
    ## replaced.
    S_clin = clinical_stay_survival(dpmf, rpmf, CFR_iso)
    O_conf_raw = two_clock_confirmed(A_bvd, conf_hazard, S_clin)
    ## `O_conf ≤ O_bvd` holds by construction; clamp into `[0, O_bvd]` as a guard
    ## under any prior draw. The suspect sub-stock is the demand remainder.
    O_conf = map((c, b) -> clamp(c, zero(eltype(demand)), b), O_conf_raw, O_bvd)
    O_susp = map((d, c) -> max(d - c, zero(eltype(demand))), demand, O_conf)
    ## Abscond outflow off the two-clock suspect stock, `κ · O_susp(t-1)`; day 1
    ## has no prior stock.
    abscond_daily = [t == 1 ? zero(eltype(demand)) : κ * O_susp[t - 1]
                     for t in 1:n]
    if eltype(demand) === Any
        demand = convert(Vector{eltype(C)}, demand)
        O_conf = convert(Vector{eltype(C)}, O_conf)
        O_susp = convert(Vector{eltype(C)}, O_susp)
        abscond_daily = convert(Vector{eltype(C)}, abscond_daily)
    end

    ## Add the reclassification offset Δ(t) to the modelled total only; demand
    ## (the diagnostic) stays the un-offset latent stock.
    occ_offset = eltype(occ_break_offset) === Any ?
                 convert(Vector{eltype(C)}, occ_break_offset) : occ_break_offset
    ## Broadcasts, not `map(1:n) do t`: the `do`-block builds an anonymous
    ## closure over `demand`/`occ_offset`/`O_conf` whose reverse-mode shadow
    ## Enzyme cannot construct; an elementwise broadcast creates no closure
    ## and is bit-identical under Mooncake.
    occ_obs_total = demand .+ occ_offset

    ## Confirmed and suspect census means, summing to the offset total.
    conf_split = copy(O_conf)
    susp_split = max.(occ_obs_total .- conf_split, zero(eltype(demand)))

    ## Occupancy likelihood: NegativeBinomial around the latent demand,
    ## right-censored at the implied-capacity bound. Days with a published split
    ## are scored as the two sub-stocks below, not the total, so the total and
    ## its parts are never both scored on one day. When the overlay is
    ## unidentified the split is unscored and those days stay on the total.
    split_days = split_active ? Set(Int.(confirmed_incare_history.days)) :
                 Set{Int}()
    iso_all_days = isolation_history.days
    iso_keep = [!(Int(d) in split_days) for d in iso_all_days]
    iso_days = iso_all_days[iso_keep]
    iso_obs = isempty(isolation_history.counts) ? missing :
              collect(Int.(isolation_history.counts))[iso_keep]
    iso_means = [occ_obs_total[clamp(Int(d), 1, n)] for d in iso_days]
    iso_ceil = censoring_cap(iso_days, iso_obs, capacity_history)
    isolation ~ to_submodel(
        censored_occupancy_model(iso_means, iso_ceil, iso_obs, k))

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

    ## Split likelihoods: on days with a published split, score the two sub-stock
    ## censuses instead of the total. Guarded by `split_active`, so they no-op
    ## when the hazard is structurally zero.
    ci_days = split_active ? confirmed_incare_history.days : Int[]
    ci_obs = (!split_active || isempty(confirmed_incare_history.counts)) ?
             missing : collect(Int.(confirmed_incare_history.counts))
    confirmed_incare_obs ~ to_submodel(vintage_increments_model(
        [conf_split[clamp(Int(d), 1, n)] for d in ci_days], ci_obs, k))
    si_days = split_active ? suspect_incare_history.days : Int[]
    si_obs = (!split_active || isempty(suspect_incare_history.counts)) ?
             missing : collect(Int.(suspect_incare_history.counts))
    suspect_incare_obs ~ to_submodel(vintage_increments_model(
        [max(susp_split[clamp(Int(d), 1, n)], zero(eltype(susp_split)))
         for d in si_days], si_obs, k))
    ## Confirmed-in-care deaths, attributed as the death flow times the confirmed
    ## share of the BVD stock. Exposed for the report, not separately scored.
    bvd_stock = acc.O_bvd
    conf_share = [bvd_stock[t] > zero(eltype(demand)) ?
                  O_conf[t] / bvd_stock[t] : zero(eltype(demand)) for t in 1:n]
    confirmed_incare_deaths_daily = deaths_daily .* conf_share

    ## Optional daily Tableau 6 flow likelihoods, each a no-op on empty history.
    ## Admissions are right-censored at the recorded free-bed headroom.
    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))
    adm_means = [admit_daily[clamp(Int(d), 1, n)] for d in adm_h_days]
    adm_ceil = admission_headroom(adm_h_days, adm_h_obs, capacity_history,
        isolation_history)
    admissions ~ to_submodel(
        censored_occupancy_model(adm_means, adm_ceil, 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. Occupancy is the censored stock — capacity
    ## bounds it at the cut-off, so the reported occupancy is the demand capped
    ## at the fixed bed count; bed demand is the uncapped latent stock.
    z0 = zero(eltype(C))
    dem_T = isempty(demand) ? z0 : demand[end]
    occ_T = min(dem_T, C_T)
    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])
    ## Unmet demand: the uncapped demand above the censored occupancy.
    bed_shortfall := safe_rate(max(dem_T - occ_T, z0))
    bed_utilisation := safe_rate(occ_T) / safe_rate(C_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
    ## Cut-off occupancy split: the two sub-stock prevalences at the grid end,
    ## and the confirmed share of the occupied stock.
    conf_incare_T = isempty(conf_split) ? z0 : conf_split[end]
    susp_incare_T = isempty(susp_split) ? z0 : max(susp_split[end], z0)
    expected_confirmed_incare := safe_rate(conf_incare_T)
    expected_suspect_incare := safe_rate(susp_incare_T)
    incare_confirmed_share := safe_rate(conf_incare_T) / safe_rate(dem_T)
    ## In-care confirmation-rate modifier ρ (raw; can exceed one, so reported
    ## directly). ρ < 1 means occupied suspects are confirmed slower than the
    ## borrowed community hazard, held for repeated exclusion testing.
    incare_confirm_modifier := ρ_conf
    ## Cut-off cumulative occupancy reclassification offset (fitted; can be
    ## negative, so reported raw rather than through `safe_rate`). Reports how
    ## much of the observed reclassification the model absorbed as a reporting
    ## artifact, the rest carried by real demand. The per-day steps `b` and the
    ## grid offset `occ_break_offset` are returned for the report.
    occupancy_break := isempty(occ_break_offset) ? z0 : last(occ_break_offset)

    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 = min.(demand, C), isolation,
        deaths_daily, recover_daily, ruleout_daily, admit_daily,
        abscond_daily,
        break_steps = b, break_offset = occ_break_offset,
        break_grid_days,
        occupancy_break = isempty(occ_break_offset) ? z0 :
                          last(occ_break_offset),
        confirmed_incare = conf_split, suspect_incare = susp_split,
        confirmed_incare_deaths_daily, incare_confirm_modifier = ρ_conf,
        expected_confirmed_incare = safe_rate(conf_incare_T),
        expected_suspect_incare = safe_rate(susp_incare_T),
        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 , 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 , the death analogue of the case ascertainment , 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 applied to the per-day non-BVD suspected-case background , 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

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

Submodel: deaths_model
julia
@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 (  plus the non-BVD background ) carried through the report-to-analysed delay and thinned by the testing fraction (the share of suspected cases routed to the laboratory),

This analysed volume is gated to zero before the testing onset: no specimens are analysed before the laboratory existed, so 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 :

The confirmed positives in each laboratory window are scored as a Binomial of the observed specimens-analysed denominator with a per-window tested-positive probability , and where no analysed count is observed (the early and unanchored windows) the modelled volume 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 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:

The tested BVD share raises by the decaying severity enrichment . A truly BVD specimen then tests positive with the sensitivity , and a non-BVD specimen with the false-positive rate  , so the false-positive term carries the non-BVD share and the laboratory data identify the background:

with the cumulative modelled laboratory volume at window , 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 , the daily modelled volume summed over the window, with the same composition-linked positivity, so all the confirmed data are used:

Submodel: lab_delay_model (receipt delay)
julia
@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
julia
@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,
        ## Opt-in retrospective harmonisation-break days (grid day-indices):
        ## days whose cumulative confirmed step is mostly a provincial base
        ## integration rather than 24h notifications. De-anchored from the
        ## positivity denominator and given a fitted level step in the
        ## modelled mean. Empty (the default) → no-op. See issue #484.
        confirmed_break_days::AbstractVector{<:Integer} = Int[],
        ## Printed 24h new-confirmed counts on each break day. The step is
        ## centred on `observed increment − gross`, the part of the vintage
        ## step the report attributes to base integration, so the magnitude
        ## comes from published data rather than a prior guess. Empty or
        ## all-zero centres on the WHOLE increment, attributing all of it to
        ## the artefact — conservative, not neutral.
        confirmed_break_gross::AbstractVector{<:Integer} = Int[],
        ## Residual uncertainty about how much of that discrepancy is truly
        ## retrospective rather than coincident same-day incidence — NOT the
        ## harmonisation magnitude, which `confirmed_break_gross` supplies.
        ## Zero pins the step at the published discrepancy and samples no
        ## parameter, taking the printed count as exact.
        confirmed_break_sd::Real = 25.0)
    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, confirmed_break_days)
    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.
        ## The per-window pool composition φ, severity upshift, assay
        ## sensitivity/specificity transform and (0,1) guards live in
        ## `composition_positivity` (a plain function with an explicit loop,
        ## so its working variables are not boxed closure captures under
        ## Enzyme).
        dscale = max(convert(Tt, decay_scale), one(Tt))
        p_pos = composition_positivity(window_days, bvd_window, bg_window,
            c_window, δ0, dscale, s_test, spec, lo, hi)
    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` is bound unconditionally to one array allocation site
    ## (a `bin_increments(...)[2:end]`) rather than a two-branch ternary.
    ## ternary. The two-allocation ternary compiles to a pointer that is one
    ## of two array allocations depending on a branch, which Enzyme's
    ## `nodecayed_phis!` LLVM pass cannot trace (an `EnzymeInternalError`; see
    ## #445). With no late/early window days the edge is the singleton
    ## `[start]`, so `bin_increments(...)[2:end]` is empty — same value as the
    ## former `similar(analysed_daily, 0)` branch. Mooncake is unaffected.
    early_edges = n_early > 0 ?
                  vcat(windows.early_start, windows.early_days) :
                  [windows.early_start]
    early_volume = bin_increments(analysed_daily, early_edges)[2:end]
    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]
    ## Unconditional single-allocation binding (see `early_volume` above):
    ## avoid the two-allocation `? bin_increments : similar` pointer-PHI that
    ## Enzyme's `nodecayed_phis!` pass cannot trace. Empty when `n_late = 0`.
    late_edges = n_late > 0 ? vcat(windows.late_start, windows.late_days) :
                 [windows.late_start]
    late_volume = bin_increments(analysed_daily, late_edges)[2:end]
    ## Opt-in retrospective harmonisation step. On a declared break day the
    ## observed increment is mostly a provincial base integration, so a single
    ## level step is fitted into that window's modelled mean and the fit
    ## partitions the increment into reporting artefact vs real incidence with
    ## uncertainty. Sampled non-centred around the published discrepancy
    ## (`break_step_centres`: the vintage increment minus the printed 24h
    ## count), with `confirmed_break_sd` the residual uncertainty about how much
    ## of that discrepancy is truly retrospective. Symmetric, because that
    ## residual can fall either side. Only break days that land
    ## on a late window can move the likelihood; others are dropped so no inert
    ## step is sampled. Same conditional-`b`-then-pure-function shape as the
    ## occupancy break (`cumulative_occupancy_offset`), which keeps the
    ## allocation single-site for the AD backends. Empty → Δ = 0, a no-op.
    late_day_ints = Int.(windows.late_days)
    conf_brk_days, conf_brk_centre = break_step_centres(late_day_ints,
        windows.late_increments, confirmed_break_days, confirmed_break_gross)
    if isempty(conf_brk_days)
        cb = Float64[]
    elseif iszero(confirmed_break_sd)
        ## `confirmed_break_sd = 0` is the DETERMINISTIC correction: the report
        ## publishes the 24h count, so the artefact size is taken as exact and
        ## the increment is fitted against the gross count with no sampled
        ## parameter. Two fewer parameters, and no ridge between a step and the
        ## ascertainment it trades off against.
        cb = conf_brk_centre
    else
        confirmed_step ~ product_distribution(
            fill(Normal(0, 1), length(conf_brk_days)))
        cb = conf_brk_centre .+ confirmed_break_sd .* confirmed_step
    end
    late_break_offset = confirmed_break_offset(windows.late_days,
        conf_brk_days, cb)
    late_mean = late_p .* late_volume .+ late_break_offset
    ## 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))

    ## Plain `=`, not `:=`: these derived quantities are surfaced through the
    ## returned NamedTuple and re-tracked at the joint level (`joint.jl`
    ## `expected_confirmed_T`/`expected_analysed_T`/`test_positivity`), so the
    ## submodel-level `:=` tracking is redundant. `:=` also builds a
    ## DynamicPPL tracking closure that captures `p_pos` (assigned in both the
    ## `:composition` and `else` positivity branches, so boxed in a
    ## `Core.Box`), and Enzyme's `nodecayed_phis!` LLVM pass cannot
    ## differentiate through the resulting boxed pointer-PHI (see #445).
    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)
    ## Unconditional broadcasts (single allocation site each): with
    ## `n_late = 0` every operand (`amask`, `late_den_a`, `late_volume`,
    ## `late_mean`) is empty, so the `ifelse.` result is empty too — the same
    ## value the former `? … : similar(…, 0)` ternary produced, without the
    ## two-allocation pointer-PHI that Enzyme's `nodecayed_phis!` pass cannot
    ## trace (see #445). Mooncake is unaffected.
    late_den = ifelse.(amask, late_den_a, late_volume)
    late_expected = ifelse.(amask, late_p .* late_den_a, late_mean)
    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: per-window tested-positive
    ## probability expanded onto the daily grid times the modelled analysed
    ## volume. In predict mode `p_pos` can widen to `Vector{Any}`, so pin it to
    ## the analysed volume's element type before expanding.
    p_pos_daily = p_pos
    if eltype(p_pos_daily) === Any
        p_pos_daily = convert(Vector{eltype(analysed_daily)}, p_pos_daily)
    end
    ## Per-day positivity `p_pos_grid`, exposed with `τ_test` so the treatment
    ## model can form the in-care confirmation hazard `τ_test · p_pos_grid[t]`.
    p_pos_grid = expand_vintage_rate(p_pos_daily, window_days, n)
    confirmed_daily = p_pos_grid .* analysed_daily

    return (; τ_test, bg_daily, p_pos, p_pos_grid, windows, analysed_daily,
        confirmed_daily,
        receipt_pmf = receipt_state.pmf,
        receipt_mean = receipt_state.mean, receipt_sd = receipt_state.sd,
        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 at the per-day suspected death-to-case ratio, times a testing-intensity scaling,

with and the modelled suspected-death and suspected-case series and 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

with and the BVD and non-BVD components of the suspected deaths (both at receipt) and , the same assay sensitivity and specificity as the confirmed cases. The false-positive term    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,

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

The death analysed volume inherits the laboratory capacity onset from the case volume , so is zero before the first confirmed-case vintage: no deaths are confirmed before the laboratory existed.

Submodel: confirmed_deaths_model
julia
@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[]),
        ## Opt-in retrospective harmonisation-break days, as in
        ## `confirmed_cases_model`: the 22 July 2026 base integration steps the
        ## confirmed-death cumulative by +236 against a printed 24h count of
        ## +62, and this stream fits between-vintage increments too, so the
        ## backlog would otherwise be read as one day of confirmed deaths.
        ## Empty (the default) → no-op. See issue #484.
        confirmed_break_days::AbstractVector{<:Integer} = Int[],
        confirmed_break_gross::AbstractVector{<:Integer} = Int[],
        ## As on the cases path: the residual uncertainty around the published
        ## discrepancy, with zero pinning the step and sampling no parameter.
        confirmed_break_sd::Real = 25.0,
        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]
            v = den > lo ? sc * case_analysed_daily[t] * susp_death[t] / den :
                zero(sc)
            ## Cap the volume at the suspected-death pool so confirmed deaths
            ## stay a subset of suspected and the realised τ_death ≤ 1.
            min(v, susp_death[t])
        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)
    ## Retrospective harmonisation step, mirroring `confirmed_cases_model`:
    ## centred on the published discrepancy (vintage increment minus the
    ## printed 24h death count) and added to that window's modelled mean, so
    ## the increment likelihood does not attribute reattached deaths to the
    ## day. Sampled whenever a break day lands on a vintage, as on the cases
    ## path, so a posterior predictive carries the same dimensions as the
    ## fitted chain and replicates the break instead of leaving the vintage as
    ## an outlier.
    cd_brk_days, cd_brk_centre = break_step_centres(vobs.days,
        vobs.obs_increments, confirmed_break_days, confirmed_break_gross)
    if isempty(cd_brk_days)
        cdb = Float64[]
    elseif iszero(confirmed_break_sd)
        ## Deterministic correction, as on the cases path: the printed 24h death
        ## count is taken as exact and no step parameter is sampled.
        cdb = cd_brk_centre
    else
        cdeath_step ~ product_distribution(
            fill(Normal(0, 1), length(cd_brk_days)))
        cdb = cd_brk_centre .+ confirmed_break_sd .* cdeath_step
    end
    modelled_inc = modelled_inc .+
                   confirmed_break_offset(vobs.days, cd_brk_days, cdb)
    ## The cut-off scalar is the generator gate, as in `confirmed_cases_model`:
    ## nulling it leaves `predict` to resample the increments while the dated
    ## history still supplies the vintage grid and the break-step centres. That
    ## keeps the published discrepancy available to a predictive, which
    ## differencing an emptied history cannot do.
    cdeath_obs = ismissing(confirmed_deaths) ? missing : vobs.obs_increments
    cdeath_increments ~ to_submodel(
        vintage_increments_model(modelled_inc, cdeath_obs, 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 (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 and convolved with a sampled confirmation-to-recovery delay ,

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   adjusted on the log-odds scale by a sampled offset  , since the confirmed cases are a slightly different population from the one the CFR is defined over,

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 :

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
julia
@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   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

Then the daily export intensity and its running sum are

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 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:

Submodel: exports_model
julia
@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 , because a death among an exported case would be reported whether or not the case itself was ascertained as an import. Writing it    , the daily export-death intensity is

Its running sum is the cumulative export-death intensity:

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 :

Submodel: exports_deaths_model
julia
@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 and the pooled ascertainment fractions, threading to the suspected-case, laboratory and confirmed-death likelihoods and 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
julia
@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
julia
@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
julia
@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
julia
@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,
        confirmed_break_days::AbstractVector{<:Integer} = Int[],
        confirmed_break_gross_cases::AbstractVector{<:Integer} = Int[],
        confirmed_break_sd::Real = 25.0,
        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, confirmed_break_days,
        confirmed_break_gross = confirmed_break_gross_cases,
        confirmed_break_sd,
        positivity_link = confirmed_positivity_link))
    ## Cut-off expected confirmed count, aliased under the same un-prefixed
    ## name [`bvd_joint`](@ref) uses so both fit kinds carry one key and the
    ## confirmed stream can be forecast from this fit
    ## ([`forecast_stream`](@ref)). It is aliased HERE, at the composer
    ## level, rather than `:=`-tracked inside
    ## [`confirmed_cases_model`](@ref): that submodel deliberately keeps its
    ## derived quantities on plain `=` because a `:=` there builds a
    ## DynamicPPL tracking closure over the boxed `p_pos` that Enzyme's
    ## `nodecayed_phis!` pass cannot differentiate through (see #445, #453).
    ## The alias below closes over the submodel's returned NamedTuple, which
    ## is assigned once and not boxed — the same pattern the joint already
    ## uses for this quantity.
    expected_confirmed_T := confirmed_state.expected_confirmed
end
Composer: exports-deaths-only fit
julia
@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
julia
@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[]),
        treatment_confirmed_incare_history = (; days = Int[], counts = Int[]),
        treatment_suspect_incare_history = (; days = Int[], counts = Int[]),
        occupancy_break_days::AbstractVector{<:Integer} = Int[],
        confirmed_break_days::AbstractVector{<:Integer} = Int[],
        confirmed_break_gross_cases::AbstractVector{<:Integer} = Int[],
        confirmed_break_gross_deaths::AbstractVector{<:Integer} = Int[],
        confirmed_break_sd::Real = 25.0,
        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,
        confirmed_positivity_link::Symbol = :composition,
        genetic = nothing,
        onset_to_sample = nejm_onset_to_sample(),
        tmrca_days::Union{Missing, Real} = missing,
        tmrca_days_sd::Real = 16.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`.
    ## The pooling SD is sampled only when the random effect is active (the
    ## tilde must stay gated), but `σ_rw_shared` and `bg_onset` are bound
    ## unconditionally to plain values so the closure below captures
    ## non-conditional, write-once variables. Conditionally-scoped captures
    ## get boxed in a `Base.RefValue`, which Enzyme's reverse mode cannot
    ## differentiate through (Mooncake tolerates it); binding them in both
    ## paths keeps the capture type-stable and the closure un-boxed.
    if background_re
        bg_pool ~ to_submodel(background_pooling_model())
        σ_rw_shared = bg_pool.σ_bg
    else
        σ_rw_shared = 0.0
    end
    ## 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. Bound unconditionally (unused when the effect is off).
    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)
    ## The closure is built unconditionally (one concrete closure type, not
    ## closure-or-`Nothing`); `background_re` then selects the closure or the
    ## `nothing` sentinel. Identical behaviour to the gated form: with
    ## `background_re = false` the closure is never passed, so the unused
    ## `σ_rw_shared = 0` never enters the log-density.
    make_case_bg = nn -> background_walk_model(nn, σ_rw_shared;
        onset = bg_onset)
    case_bg_re = background_re ? make_case_bg : nothing

    ## Cases first so the suspected-case background `bg_daily` is available to
    ## the deaths stream (which scales it by `cfr_bg` for the death background)
    ## and to the laboratory pipeline.
    cases_state ~ to_submodel(
        cases(reported_history, reported_cases, onsets, k_cases, p_drc;
        suspected_daily_history, background_re = case_bg_re))
    deaths_state ~ to_submodel(
        deaths(deaths_history, total_deaths, onsets, k_deaths;
        suspected_daily_deaths_history, case_bg_daily = cases_state.bg_daily))
    confirmed_state ~ to_submodel(
        confirmed(confirmed_history, confirmed_cases, onsets, k_confirmed,
        p_drc, cases_state.bg_daily, cases_state.τ_test,
        cases_state.bvd_reports_daily;
        lab_history, lab_daily_history,
        tests_analysed, confirmed_break_days,
        confirmed_break_gross = confirmed_break_gross_cases,
        confirmed_break_sd,
        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,
        confirmed_break_days,
        confirmed_break_gross = confirmed_break_gross_deaths,
        confirmed_break_sd,
        case_analysed_daily = confirmed_state.analysed_daily,
        case_suspected_daily = cases_state.reports_daily))
    ## Treatment-centre patient flow ([`treatment_flow_model`](@ref)): occupancy
    ## plus the in-care outcome flows, with the in-care fatality CFR_iso
    ## identified by the in-care death flow. The occupancy split borrows the
    ## in-care confirmation hazard `τ_test · p_pos` from the confirmed pipeline to
    ## carve the occupied true-case stock into confirmed and suspect sub-stocks,
    ## scored against the Tableau 6 `dont confirmés` / `dont suspects` census. The
    ## known DHIS2 harmonisation days carry the overnight total reporting break.
    conf_hazard_daily = confirmed_state.τ_test .* confirmed_state.p_pos_grid
    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,
        confirmed_incare_history = treatment_confirmed_incare_history,
        suspect_incare_history = treatment_suspect_incare_history,
        occupancy_break_days = occupancy_break_days,
        conf_hazard_daily = conf_hazard_daily,
        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)
    ## External onset-to-sample constraint on the confirmed sampling delay
    ## (grounded on the NEJM DRC 2026 cohort by default, see
    ## [`nejm_onset_to_sample`](@ref)). The onset→report and report→receipt legs
    ## convolve to the confirmed onset-to-sample delay, so its continuous mean is
    ## the sum of the two legs' means and its continuous SD the root-sum of their
    ## variances; both are exposed here. The cohort's reported (continuous) mean
    ## and median are fitted to these as soft Normal observations
    ## ([`onset_to_sample_logweight`](@ref)), grounding the otherwise-
    ## unidentified receipt (lab-turnaround) leg without touching either prior.
    ## The term only exists on the confirmed report⊕receipt path, so single-
    ## stream and isolation composers carry none; passing `nothing` drops it.
    onset_to_sample_mean := cases_state.report_mean +
                            confirmed_state.receipt_mean
    onset_to_sample_sd := sqrt(cases_state.report_sd^2 +
                               confirmed_state.receipt_sd^2)
    if onset_to_sample !== nothing
        @addlogprob! onset_to_sample_logweight(cases_state.report_mean,
            cases_state.report_sd, confirmed_state.receipt_mean,
            confirmed_state.receipt_sd, onset_to_sample)
    end
    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 occupancy split: the confirmed-in-care and suspect-in-care
    ## sub-stock prevalences carved from the occupied true-case stock by the
    ## confirmation overlay.
    expected_confirmed_incare_T := treatment_state.expected_confirmed_incare
    expected_suspect_incare_T := treatment_state.expected_suspect_incare
    ## 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 outcome mixture: `isolation_bvd_los_mean` is 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
    ## In-care confirmation-rate modifier ρ on the borrowed community
    ## confirmation hazard, identified by the confirmed/suspected-in-care split.
    incare_confirm_modifier := treatment_state.incare_confirm_modifier
    expected_recovered_T := recovered_state.expected_recovered
    recovery_probability := recovered_state.p_recover
    recovery_delay_mean := recovered_state.recovery_delay_mean
    recovered_dispersion := recovered_state.k_recovered
    tau_test := cases_state.τ_test
    lambda_bg := cases_state.λ_bg
    bg_sigma := cases_state.bg_sigma
    background_total := cases_state.bg_total
    death_ascertainment := deaths_state.p_death
    background_cfr := deaths_state.cfr_bg
    lambda_bg_death := deaths_state.λ_bg_death
    bg_death_sigma := deaths_state.bg_death_sigma
    background_death_total := deaths_state.bg_death_total
    tau_death := confirmed_deaths_state.τ_death
    death_testing_scaling := confirmed_deaths_state.scaling
    suspected_positivity := cases_state.positivity
    test_positivity := confirmed_state.p_positive
    death_composition := confirmed_deaths_state.q_death
    death_confirmation := confirmed_deaths_state.p_death_conf
end

Model fitting and evaluation ​

Prior predictive check ​

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

Sample the joint prior
julia
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
RowQuantityLower 90%Lower 60%Lower 30%Upper 30%Upper 60%Upper 90%
StringFloat64Float64Float64Float64Float64Float64
1C_T3291.024125.098348.01.3336e67.70433e62.07528e8

Pair plot of the prior over the latent quantities.

Prior pair plot
julia
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
10×4 DataFrame
Rowfitmax_rhatmin_ess_bulkdivergences
StringFloat64Float64Int64
1joint1.014327.031
2exports1.006302.00
3deaths (DRC)1.015300.08
4cases (DRC)1.019379.035
5confirmed (DRC)1.016301.014
6confirmed deaths (DRC)1.01515.038
7isolation (DRC)1.018391.051
8frozen (1wk back)1.018246.013
9delay sensitivity1.016261.05
10clock sensitivity (ExpGrowth)1.008460.01

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,    , where 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:

with the cumulative confirmed deaths, the modelled daily confirmed-case incidence, and   the residual delay between a confirmed case and its confirmed death. is the onset-to-death-confirmation lag (onset-to-death convolved with the report-to-receipt laboratory delay) and 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 DRC observation streams as forecast targets: the reported cases and suspected deaths, the laboratory-confirmed cases and confirmed deaths, the isolation/treatment beds and the 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 reported case and suspected death streams are no longer published, so their forecasts extend the last published cumulative total rather than a still-growing series. Exports are not forecast, 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.

Each release now saves its forecast as an asset so it can later be scored against what is observed. Earlier releases showed a forecast but did not store it, so those forecasts are reconstructed. scripts/backfill_forecasts.jl checks out each past release at its own tag, re-runs that release's own model code on that release's data through its own fit, and writes the forecast in the same archive schema, so a reconstructed forecast is the release's own output rather than a current-code approximation. The reconstruction re-resolves dependencies at current versions, since the release manifests were not pinned, so it reproduces the release's model and data but not its exact solver build. Reconstruction reaches back to the first release that carried any forecast, so it covers the whole release history, but the streams available grow over time as the model and data did. The renewal releases from v1.4.0 reconstruct the incident case and death streams, extending to all four streams from v1.6.0 once the recovered and isolation series entered the data. The integral-era v1.3.0 reconstructs the confirmed case and death streams. The earliest releases, v1.0.0 to v1.2.0, forecast the reported case, suspected death, and export streams, and are reconstructed from each tag's own inline model code. Reconstructed forecasts are published as a separate backfill release and scored alongside the stored ones by the forecast scoring across releases section.

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)
julia
# The frozen re-fits are defined in the fit registry (`docs/fits/registry.jl`) and
# loaded through the cache in the setup block above.

Forecast scoring against a persistence baseline ​

Every forecast above, and every stored forecast from a past release, is scored with the continuous ranked probability score, on the count scale and on a log scale that stops the largest counts dominating. We report the score split into the predictive spread and the cost of reading high or low, the share of observations inside the 50% and 90% predictive intervals, and a bias running from when every forecast sits below the observation to when every one sits above it. Relative skill between fits and is

each mean taken over the forecasts both fits scored, so a comparator that happens to score zero on one forecast cannot send the ratio to infinity. A value below one beats the comparator.

The count streams are running cumulative totals, so each is scored on its increment over the forecast window rather than on the level it reaches, while bed occupancy is a level and is scored as one. A stream is scored only where its own reporting covers the window, from the day it was first reported to the day it was last updated, since outside that period a cumulative total that has not moved is the absence of a series rather than an observed zero. On the two confirmed streams a reported step that is mostly a retrospective integration of harmonised provincial records has that backfill removed from both the target and the baseline, so such a window is scored on transmission rather than on transmission plus an integration.

Each forecast is also compared against a persistence baseline built from the same stream. Write the vintages recorded by the day the forecast was made as dates    carrying cumulative values , let be the value at the last vintage on or before , and let be the day the forecast was made and the horizon in days. The baseline centres on the occupancy reached, or on the increment over the preceding window of the same length,

and takes its spread from the record's own first differences, each rescaled to a one-day step and entered with both signs,

Under a driftless walk of per-day variance a change over days has variance , so dividing by puts vintages recorded at different spacings on a common one-day scale, and holding both signs makes the pool mean zero so a record that only rises does not give the walk a direction. One predictive draw iterates the walk to the horizon,

so before the floor it has mean and variance for  , and the interval widens with the square root of the horizon. This follows the COVID-19 Forecast Hub baseline, except that the centre for a count stream pools the whole window rather than the single most recent increment, since these vintages are sparse and irregularly spaced. Fewer than three recorded differences leaves no usable pool and the baseline falls back to a Poisson draw around the centre.

The baseline reads only the vintages recorded by the day the forecast was made, from the archived data snapshot the release itself was built on, so no later correction or backfill reaches it. For a forecast made at a release's own cut-off that snapshot is the one the forecast was made from and the guarantee is exact. The frozen re-fits below forecast from fixed historical cut-offs reused across later releases, so their snapshot can post-date the day the forecast was made by weeks, and a correction landing in between is already in it. Closing that would need a snapshot archived per frozen cut-off, which does not exist. The earliest releases archived their cut-off totals without the dated vintage record, leaving their baseline no history to draw on, so skill against it is not informative there.

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
julia
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% 6907–8028, 60% 6375–8993, 90% 5702–11219 infections to date, reported and unreported.

  • Against the 3748 laboratory-confirmed cases by the cut-off that is roughly 1.5–3.0× 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-15–2026-03-22, 60% 2026-03-11–2026-03-25, 90% 2026-03-01–2026-03-27, an elapsed age to the cut-off of 30% 132–139, 60% 129–143, 90% 127–153 days.

  • Growth rate and doubling time: the initial growth rate is estimated to have been 30% 0.054–0.064, 60% 0.049–0.07, 90% 0.04–0.082 per day, an initial doubling time of 30% 10.8–12.8, 60% 9.9–14.2, 90% 8.4–17.4 days. The latest growth rate is estimated to be 30% 0.004–0.013, 60% -0.003–0.019, 90% -0.017–0.036 per day, a latest doubling time of 30% 25.4–60.4, 60% -47.0–99.1, 90% -306.1–301.2 days.

  • Reproduction number: the initial reproduction number is estimated to have been 30% 1.76–2.16, 60% 1.6–2.41, 90% 1.36–3.1 and the latest to be 30% 1.04–1.17, 60% 0.97–1.28, 90% 0.81–1.5.

  • Case-fatality ratio: the case-fatality ratio is estimated to be 30% 0.46–0.51, 60% 0.42–0.54, 90% 0.37–0.59.

  • 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.09, the outbreak age by -0.87 and the doubling time by 1.56; the largest move is in the doubling time.

Joint model estimates ​

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

Cumulative infection count summary table
julia
cumulative_cases_summary = summary_table(
    chn_joint, [:C_T]; digits = 0);
1×7 DataFrame
RowQuantityLower 90%Lower 60%Lower 30%Upper 30%Upper 60%Upper 90%
StringFloat64Float64Float64Float64Float64Float64
1C_T5702.06375.06907.08028.08993.011219.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
julia
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
julia
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
julia
infection_summary = summary_table(chn_joint,
    [:r, :doubling_time, :T, :R_T, :CFR, :C_T]; digits = 2);
6×7 DataFrame
RowQuantityLower 90%Lower 60%Lower 30%Upper 30%Upper 60%Upper 90%
StringFloat64Float64Float64Float64Float64Float64
1r-0.02-0.00.00.010.020.04
2doubling_time-306.07-46.9725.3760.4299.14301.23
3T126.64129.35132.14138.67143.32152.58
4R_T0.810.971.041.171.281.5
5CFR0.370.420.460.510.540.59
6C_T5702.26375.326906.978027.78993.311219.5
Infection-parameter pair plot (prior overlaid)
julia
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
julia
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
RowQuantityLower 90%Lower 60%Lower 30%Upper 30%Upper 60%Upper 90%
StringFloat64Float64Float64Float64Float64Float64
1generation interval shape1.432.022.322.773.063.64
2generation interval scale3.014.224.935.966.697.76
3incubation period mean5.596.06.226.596.837.31
4incubation period SD2.152.763.153.744.074.69
Infection-delay pair plot (prior overlaid)
julia
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
julia
# `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 = RT_INTERVENTION_RAMP);

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
julia
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
RowStreamLower 90%Lower 60%Lower 30%Upper 30%Upper 60%Upper 90%
StringFloat64Float64Float64Float64Float64Float64
1Rt multiplier exp(effect)0.450.580.650.790.860.96

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
julia
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
RowQuantityLower 90%Lower 60%Lower 30%Upper 30%Upper 60%Upper 90%
StringFloat64Float64Float64Float64Float64Float64
1onset-to-report shape0.630.830.931.141.261.46
2onset-to-report scale1.412.132.583.493.914.67
3onset-to-admission shape0.70.961.091.321.461.72
4onset-to-admission scale1.932.833.344.34.95.76
5admission-to-death shape1.491.872.092.512.783.19
6admission-to-death scale2.363.263.864.785.376.61
7onset-to-detection shape0.740.981.121.321.431.65
8onset-to-detection scale2.082.953.474.44.855.71
9report-to-receipt mean2.53.153.454.064.455.0
10report-to-receipt SD2.943.513.84.394.715.34
11in-care BVD length-of-stay mean (mixture)8.449.229.6910.6411.212.25
12isolation non-BVD rule-out stay mean6.667.257.648.138.499.11
13confirmation-to-recovery mean12.7115.4417.1320.0121.4223.73
Observation-delay pair plot (prior overlaid)
julia
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: is the population-level dispersion, , , and 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 (, ) 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
julia
surveillance_summary = summary_table(chn_joint,
    [:p_drc, :p_uganda, :k, :k_cases, :k_deaths, :k_confirmed,
        :k_confirmed_deaths, :dispersion_sd, :tau_test, :lambda_bg,
        :suspected_positivity, :test_positivity, :expected_confirmed_T,
        :expected_analysed_T, :death_ascertainment, :background_cfr,
        :tau_death, :death_composition,
        :death_confirmation, :expected_confirmed_deaths_T,
        :isolation_admission, :isolation_dispersion, :expected_isolation_T,
        :expected_bed_demand_T, :bed_capacity, :bed_shortfall_T,
        :incare_cfr, :incare_cfr_modifier, :incare_confirm_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
35×7 DataFrame
RowQuantityLower 90%Lower 60%Lower 30%Upper 30%Upper 60%Upper 90%
StringFloat64Float64Float64Float64Float64Float64
1p_drc0.4940.6230.6990.820.8710.945
2p_uganda0.4780.6460.7190.8290.8820.945
3k2.9234.0294.9316.77.7359.785
4k_cases6.3517.5858.4289.84510.85212.536
5k_deaths6.1547.8218.9311.09912.71116.131
6k_confirmed3.153.7334.0844.6995.1896.11
7k_confirmed_deaths2.2692.7643.0643.764.2375.124
8dispersion_sd0.4020.50.5610.690.790.977
9tau_test0.8760.9140.9310.9590.9710.987
10lambda_bg14.07717.34919.33223.47625.88429.837
11suspected_positivity0.240.260.2710.2880.2980.321
12test_positivity0.2640.2780.2840.2970.3050.319
13expected_confirmed_T3653.213834.073924.024093.764197.924396.69
14expected_analysed_T13468.714011.414259.714760.015039.015559.8
15death_ascertainment0.7760.8430.8740.9130.9320.954
16background_cfr0.1840.2480.2830.340.3760.437
17tau_death0.5360.6250.6820.8010.891.0
18death_composition0.2690.3190.3590.4390.4960.604
19death_confirmation0.280.3190.360.4270.4840.578
20expected_confirmed_deaths_T1240.61321.051366.951441.261486.171581.29
21isolation_admission0.2870.3030.3110.3250.3330.351
22isolation_dispersion51.95661.49567.49381.19590.959115.822
23expected_isolation_T688.167722.922743.912779.809798.313836.398
24expected_bed_demand_T688.167722.922743.912779.809798.313836.398
25bed_capacity953.583993.9761021.191069.251097.381160.52
26bed_shortfall_T0.00.00.00.00.00.0
27incare_cfr0.4190.4510.470.5070.530.582
28incare_cfr_modifier-0.441-0.217-0.0730.1630.3060.536
29incare_confirm_modifier0.6110.8060.9291.2551.5142.116
30isolation_death_los_mean2.5114.4235.6787.6818.94611.337
31isolation_recovery_los_mean8.42810.93612.20614.67415.97118.546
32abscond_fraction0.010.0110.0110.0120.0130.014
33recovery_probability0.2390.2710.2910.330.3520.395
34recovered_dispersion2.2782.8623.2063.894.3585.492
35expected_recovered_T668.903716.778746.926796.32833.644904.006
Surveillance-parameter pair plot (prior overlaid)
julia
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
julia
# 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),
        treatment_confirmed_incare_history =
        _days_only(obs.treatment_confirmed_incare_history),
        treatment_suspect_incare_history =
        _days_only(obs.treatment_suspect_incare_history),
        confirmed_history = obs.confirmed_history,
        # Counts kept, like the confirmed cases above: the cut-off scalar
        # (`confirmed_deaths = missing`) is this stream's generator gate, so
        # `predict` still resamples the increments while the dated history
        # supplies both the vintage grid and the published break discrepancy
        # the step is centred on. Differencing an emptied history cannot
        # recover that discrepancy, which would leave the harmonised vintage
        # replicated as a day of real deaths.
        confirmed_deaths_history = obs.confirmed_deaths_history,
        lab_history = obs.lab_history,
        lab_daily_history = obs.lab_daily_history,
        # Kept, like the occupancy break above, so the generator's confirmed
        # break dimension matches the fitted chain (the level step and the
        # de-anchored positivity denominator on the
        # `[confirmed_break_dates]` days). Without them the harmonised
        # vintage is replicated as though its whole increment were one day of
        # incidence, so 22 July plots as a gross outlier against a chain that
        # fitted it as mostly backlog, and the `confirmed_step` columns go
        # unused.
        confirmed_break_days = obs.confirmed_break_days,
        confirmed_break_gross_cases = obs.confirmed_break_gross_cases,
        confirmed_break_gross_deaths = obs.confirmed_break_gross_deaths,
        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.
# The treatment model scores the total occupancy only on the days without a
# published confirmed/suspect split (a per-day total-or-split switch): on the
# split days the two sub-stock census panels carry the fit instead, so the
# `isolation.obs` predictive holds only the non-split days. Drop the split
# days from the panel's dates and observed counts to match that length.
_iso_split_days = Set(Int.(obs.treatment_confirmed_incare_history.days))
_iso_keep = [!(Int(d) in _iso_split_days) for d in obs.isolation_history.days]
isolation_panel = (;
    title = "Patients in isolation",
    dates = _vintage_dates(obs.isolation_history.days[_iso_keep]),
    replicates = _vintage_replicates(
        pp_joint, @varname(isolation.obs)),
    observed = obs.isolation_history.counts[_iso_keep],
    colour = :darkorange, cumulative = false);
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.obs)),
    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);

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

# 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,
    confirmed_incare_panel, suspect_incare_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
julia
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.

julia
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
julia
stream_calibration_fig = plot_stream_calibration(stream_calibration_table);
Per-stream calibration table
16×5 DataFrame
RowStreamVintagesBias50% coverage90% coverage
StringInt64Float64Float64Float64
1Suspected cases90.070.440.67
2New suspects/day530.010.70.98
3Patients in isolation160.110.560.94
4Confirmed cases700.050.590.93
5Suspected deaths9-0.070.560.78
6New suspected deaths/day32-0.070.620.97
7Confirmed deaths720.040.60.93
8Recovered (confirmed)510.050.510.92
9Specimens analysed (cumulative)6-0.170.00.67
10Specimens analysed (24h)53-0.150.740.94
11Admissions/day41-0.010.510.85
12In-care deaths/day410.130.540.85
13Rule-outs/day400.030.280.68
14Absconded/day410.070.390.88
15Confirmed in care37-0.030.951.0
16Suspects in care37-0.080.811.0

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

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

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

Posterior correlations and stream totals ​

The heatmap is the posterior correlation between each pair of headline quantities: the outbreak size (), the reproduction number (), the outbreak age (), the case-fatality ratio (CFR), the DRC and Uganda ascertainment fractions (, ), the non-BVD background rate (), the fraction tested (), and the cut-off total expected for each stream. Blue is positive, red negative.

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

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

Counterfactual: lower bound under no further transmission ​

The committed future deaths if transmission stopped at the report date, defined in the methods counterfactual.

Project no-onward deaths and summarise
julia
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
RowStreamLower 90%Lower 60%Lower 30%Upper 30%Upper 60%Upper 90%
StringFloat64Float64Float64Float64Float64Float64
1no-onward total246.0246.0246.0246.0246.0608.0

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

No-onward projected-deaths plot
julia
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
julia
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;
julia
4×3 DataFrame
RowQuantityCentral estimateNarrowest interval
StringStringString
1Delay-corrected confirmed CFR45.0%37.4–54.7%
2Structural (infection-based) CFR48.1%36.9–59.2%
3Uncorrected modelled confirmed ratio35.0%30.2–40.7%
4Naive observed confirmed ratio44.2%—

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
julia
confirmed_cfr_fig = plot_confirmed_cfr(confirmed_cfr);

One-week-ahead forecast results ​

The cumulative and new expected counts by   from the no-change projection defined in the methods one-week-ahead forecast. The summary table reports the confirmed case and death streams, the recovered total and the isolation-bed levels and daily flows; the observed-forecast plot below additionally shows the reported cases and suspected deaths, so every projected stream appears.

Generate the one-week-ahead forecast
julia
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
RowStreamQuantityLower 90%Lower 60%Lower 30%Upper 30%Upper 60%Upper 90%
StringStringFloat64Float64Float64Float64Float64Float64
1DRC confirmed casescumulative by T+73981.04125.04235.04439.04588.04920.0
2DRC confirmed casesnew this week233.0377.0487.0691.0840.01172.0
3DRC confirmed deathscumulative by T+7292.0403.0504.0789.01021.01848.0
4DRC confirmed deathsnew this week0.00.00.00.00.0191.0
5DRC isolation bedsdemand at T+7569.0679.0755.0878.0962.01142.0
6DRC isolation bedsoccupancy at T+7569.0679.0755.0878.0952.01039.0
7DRC isolation admissionsdaily at T+790.0109.0120.0142.0157.0189.0
8DRC in-care deathsdaily at T+716.021.024.029.033.041.0
9DRC isolation rule-outsdaily at T+755.067.074.087.096.0115.0
10DRC recoveredcumulative by T+7714.0730.0744.0792.0838.0981.0
11DRC recoverednew this week6.022.036.084.0130.0273.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
julia
forecast_latent_fig = plot_forecast_latent(forecast);

The observed figure shows the new count each reported stream adds over the horizon: reported cases, suspected deaths, laboratory-confirmed cases, confirmed deaths and recovered, one panel per stream the forecast carries.

One-week-ahead observed forecast plot
julia
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
julia
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
julia
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, the one- to four-week-ahead forecasts of the observed streams (so each release records the forecast it made, for later scoring), and a copy of the input observations.toml so the exact data that produced each result is recorded alongside it.

Write outputs to output/
julia
# 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]
# Raw walk base `log_R0`, the renewal reproduction-number walk's starting
# point on the log scale, kept unexponentiated so downstream scoring takes
# its own exp. This is a distinct quantity from `r0`, the growth-clock
# initial rate, so both columns are published side by side.
log_R0_draws = vec(Array(chn_joint[Symbol("rt_state.log_R0")]))
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
)
posterior_draws[!, Symbol("rt_state.log_R0")] = log_R0_draws
posterior_draws = posterior_draws[1:10:end, :]
CSV.write(joinpath(output_dir, "posterior_draws.csv"), posterior_draws);

# One- to four-week-ahead forecasts of the observed streams, saved as a
# release asset so each release records the forecast it made and it can later
# be scored against what is observed. Only the incident and level quantities
# are archived (see `forecast_archive`), thinned to keep the asset compact.
forecast_horizons = (7, 14, 21, 28)
forecast_runs = [(h,
                     forecast_reported(chn_joint; horizon = h,
                         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))
                 for h in forecast_horizons]
CSV.write(joinpath(output_dir, "forecast.csv"),
    forecast_archive(forecast_runs; made_date = obs.cutoff, thin = 5));

# The same one- to four-week-ahead forecast made from each FROZEN joint
# re-fit (the McCabe-matched cut-offs, the Chamla anchor and the one-week-back
# validation fit), each stamped with its OWN cut-off as the made date. This
# gives historical forecast evaluation using the current model at past data
# cut-offs, scored against what has since been observed, without
# reconstructing old release tags. Each frozen fit uses its own frozen
# observations for the cut-off counts. The May cut-offs predate the isolation
# and recovered streams, so those are simply absent for them; the per-stream
# guard in `forecast_archive` skips a stream a fit does not carry.
frozen_forecast_fits = unique(f -> f.o.cutoff,
    [frozen_results; frozen_by_cutoff[chamla_cutoff]; frozen_lastweek])
frozen_forecast_archive = DataFrame(made_date = Date[], horizon = Int[],
    target_date = Date[], stream = String[], draw = Int[], value = Float64[])
for f in frozen_forecast_fits
    runs = [(h,
                forecast_reported(f.chn; horizon = h,
                    obs_cases = f.o.reported_cases,
                    obs_deaths = f.o.total_deaths,
                    obs_confirmed = f.o.confirmed_cases,
                    obs_confirmed_deaths = f.o.confirmed_deaths,
                    obs_recovered = f.o.recovered_cases))
            for h in forecast_horizons]
    append!(frozen_forecast_archive,
        forecast_archive(runs; made_date = f.o.cutoff, thin = 5))
end
CSV.write(joinpath(output_dir, "forecast_frozen.csv"),
    frozen_forecast_archive);

# Per-fit release assets: the reproduction number, outbreak size and forecasts
# for every fit rather than the joint alone, so a release records what each
# dataset implies on its own and can later be scored against the joint. The
# single-stream fits walk Rt from day 1 while the joint walks from
# `RT_WALK_LEAD` days before the first situation report, so each fit carries
# the starts its own fit used. Each single-stream fit forecasts only the
# dataset it observes; the joint forecasts every shared stream. Recovered has
# no single-stream fit, so it stays a joint-only stream in `forecast.csv`.
stream_thin = 5
_rt_walk_start_joint = clamp(_BREAKPOINT - RT_WALK_LEAD, _rt_start_plot, obs.n)
# Observed bed occupancy at the cut-off, the level the isolation forecast
# anchors on.
_iso_at_cutoff = isempty(obs.isolation_history.counts) ? 0 :
                 obs.isolation_history.counts[end]
stream_fits = [
    (; fit = "joint", chn = chn_joint, rt_start = _rt_start_plot,
        rt_walk_start = _rt_walk_start_joint,
        streams = [(:reported_cases, "reported cases", obs.reported_cases),
            (:suspected_deaths, "suspected deaths", obs.total_deaths),
            (:confirmed_cases, "confirmed cases", obs.confirmed_cases),
            (:confirmed_deaths, "confirmed deaths", obs.confirmed_deaths),
            (:isolation_beds, "isolation beds", _iso_at_cutoff),
            (:exports, "exports", obs.exported_cases)]),
    (; fit = "cases", chn = chn_cases, rt_start = 1, rt_walk_start = 1,
        streams = [(:reported_cases, "reported cases", obs.reported_cases)]),
    (; fit = "deaths", chn = chn_deaths, rt_start = 1, rt_walk_start = 1,
        streams = [(:suspected_deaths, "suspected deaths", obs.total_deaths)]),
    (; fit = "confirmed", chn = chn_confirmed, rt_start = 1, rt_walk_start = 1,
        streams = [(:confirmed_cases, "confirmed cases", obs.confirmed_cases)]),
    (; fit = "confirmed_deaths", chn = chn_confirmed_deaths, rt_start = 1,
        rt_walk_start = 1,
        streams = [(:confirmed_deaths, "confirmed deaths",
            obs.confirmed_deaths)]),
    (; fit = "treatment", chn = chn_treatment, rt_start = 1,
        rt_walk_start = 1,
        streams = [(:isolation_beds, "isolation beds", _iso_at_cutoff)]),
    (; fit = "exports", chn = chn_exports, rt_start = 1, rt_walk_start = 1,
        streams = [(:exports, "exports", obs.exported_cases)])
]

# Cut-off reproduction number per fit. The joint exposes it as `R_T`; the
# single-stream composers do not (the alias lives in `bvd_joint`, not the
# shared latent submodel), so theirs is rebuilt from the walk parameters every
# chain carries and read at the cut-off, the last day of the reconstructed
# path. `ramp` matches the model's 21-day intervention scale-up.
function _fit_rt_draws(f)
    f.fit == "joint" && return vec(Array(f.chn[:R_T]))
    rt = reconstruct_rt(f.chn; n = obs.n, breakpoint = _BREAKPOINT,
        rt_start = f.rt_start, rt_walk_start = f.rt_walk_start, ramp = RT_INTERVENTION_RAMP)
    return Float64[rt[i, obs.n] for i in axes(rt, 1)]
end

# Basic reproduction number per fit, the renewal walk's starting value, a
# single distribution per fit rather than a daily series. Every current
# single-stream composer shares the joint's latent submodel and so carries
# its own walk base, but `r0_walk_draws` probes rather than assumes, so a
# single-stream model built without its own renewal walk drops out of this
# quantity instead of breaking the release.
_stream_quantities = [(f.fit, _fit_rt_draws(f), vec(Array(f.chn[:C_T])),
                          r0_walk_draws(f.chn)) for f in stream_fits]

# One row per fit and quantity, with the median and the 30/60/90% credible
# bounds the report's tables use.
function _stream_estimate_row(fit, quantity, draws)
    s = posterior_summary(draws)
    return (fit = fit, quantity = quantity, median = quantile(draws, 0.5),
        lo30 = s.lo30, hi30 = s.hi30, lo60 = s.lo60, hi60 = s.hi60,
        lo90 = s.lo90, hi90 = s.hi90)
end
stream_estimates = DataFrame([_stream_estimate_row(fit, q, d)
                              for (fit, rt, ct, r0) in _stream_quantities
                              for (q, d) in (("R_T", rt), ("C_T", ct),
                                      ("R0", r0))
                              if !isnothing(d)])
CSV.write(joinpath(output_dir, "stream_estimates.csv"), stream_estimates);

# Thinned reproduction-number and outbreak-size draws per fit, so downstream
# scoring can recompute its own summaries rather than reuse the intervals.
stream_draws = DataFrame([(fit = fit, quantity = q, draw = d, value = v)
                          for (fit, rt, ct, r0) in _stream_quantities
                          for (q, vals) in (("R_T", rt), ("C_T", ct),
                                  ("R0", r0))
                          if !isnothing(vals)
                          for (d, v) in enumerate(vals[1:stream_thin:end])])
CSV.write(joinpath(output_dir, "stream_draws.csv"), stream_draws);

# Per-fit forecasts of each fit's own observed stream, in the `forecast.csv`
# long schema plus the fit that made them. Rebuilding a single-stream fit's
# cut-off growth rate needs the grid length and the breakpoint, which are data
# rather than chain contents, so both are passed.
stream_forecasts = DataFrame(made_date = Date[], horizon = Int[],
    target_date = Date[], stream = String[], draw = Int[], value = Float64[],
    fit = String[])
for f in stream_fits, (stream, label, obs_value) in f.streams,
    h in forecast_horizons
    _vals = forecast_stream(f.chn, stream; horizon = h,
        obs_value = obs_value, n = obs.n, breakpoint = _BREAKPOINT,
        rt_start = f.rt_start, rt_walk_start = f.rt_walk_start)
    for (d, i) in enumerate(1:stream_thin:length(_vals))
        push!(stream_forecasts, (obs.cutoff, h, obs.cutoff + Day(h), label,
            d, Float64(_vals[i]), f.fit))
    end
end

# Confirmed/suspect ward-bed occupancy for the joint, taken from the joint's
# own `forecast_reported` runs (partitioned by the cut-off confirmed share in
# `forecast_reported`, not re-derived here) so the ward beds are scored on the
# same footing as the total occupancy in the preferred `stream_forecasts.csv`
# asset. `forecast_stream` cannot project these — the split is not a growable
# stream but a partition of the total — so they are read from the archive
# runs. Dormant until the chain carries the confirmed in-care split:
# `forecast_reported` emits these columns only when the confirmed in-care
# prevalence `expected_confirmed_incare_T` is present, so the guard skips
# them otherwise.
for (h, fc) in forecast_runs,
    (col, label) in ((:suspect_occupancy, "isolation beds (suspected)"),
        (:confirmed_occupancy, "treatment beds"))

    col in propertynames(fc) || continue
    _wvals = fc[!, col]
    for (d, i) in enumerate(1:stream_thin:length(_wvals))
        push!(stream_forecasts, (obs.cutoff, h, obs.cutoff + Day(h), label,
            d, Float64(_wvals[i]), "joint"))
    end
end
CSV.write(joinpath(output_dir, "stream_forecasts.csv"), stream_forecasts);

# 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
julia
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.