BVDOutbreakSize

Appearance

Observation models ​

One submodel per reported data stream, mapping the latent infections onto what surveillance recorded. Each takes the latent trajectory and the stream's observations and adds that stream's likelihood.

Index ​

Reference ​

BVDOutbreakSize.break_step_centres Method
julia
break_step_centres(
    window_days,
    increments,
    break_days,
    gross
) -> Tuple{Vector{Int64}, Vector{Float64}}
julia
break_step_centres(window_days, increments, break_days, gross)

Break days that land on a scored window, paired with the data-derived centre for each one's fitted step. The centre is observed increment − gross, the part of that vintage's step the report itself attributes to integrating a harmonised provincial base rather than to the day's own notifications, with gross the printed 24h new-confirmed count. On 22 July 2026 that is 369 − 97 = 272 cases and 236 − 62 = 174 deaths. Centring on the published discrepancy rather than zero is what makes the step identifiable from a single observation. The sampled deviation around the centre then carries only the residual uncertainty about how much of the discrepancy is truly retrospective.

Break days not matching a window are dropped so no inert step is sampled. A missing or short gross is treated as zero, so the centre becomes the whole increment and the entire vintage step is attributed to the artefact. That is conservative rather than neutral, since the de-anchor drops the day's positivity denominator and a step centred near zero over a de-anchored day leaves nothing to absorb the backlog. Supply the printed count so the split is data-derived.

increments may itself be missing, the generator path where no cumulative counts are supplied to difference. The day still gets a step, so a predictive keeps the fitted chain's dimensions, but there is no published discrepancy to centre it on and the centre falls back to zero. Returns (days, centres), both possibly empty.

source
BVDOutbreakSize.composition_shares Method
julia
composition_shares(
    weight::AbstractVector,
    modelled::AbstractMatrix
) -> Any
julia
composition_shares(weight, modelled)

Expected share of each patch at each vintage: the modelled per-patch increments modelled (n_patches × n_vintages) weighted by the per-patch weight and normalised within each vintage. safe_rate floors the modelled increments away from zero so a vintage with no modelled cases in a patch still gives a defined (tiny) share rather than a 0/0.

source
BVDOutbreakSize.composition_split_model Method
julia
composition_split_model(
    obs_increments::Union{Missing, AbstractMatrix{<:Integer}},
    shares::AbstractMatrix,
    totals::AbstractVector{<:Integer},
    rho::Real
) -> Any

Split each vintage's totals across the patches by stick-breaking at the expected shares with the overdispersed Binomial safe_betabinomial at overdispersion rho (see province_composition_model). Patch p is drawn from the cases patches 1 … p−1 left, at its share among the patches not yet allocated, and the last patch takes the remainder. obs_increments is the observed split, or missing to draw one, which is how a forecast splits a national forecast across the provinces. On that path the last patch's remainder is filled in too, so the returned matrix sums to totals in every column.

source
BVDOutbreakSize.confirmed_break_offset Method
julia
confirmed_break_offset(late_days, break_days, b) -> Any
julia
confirmed_break_offset(late_days, break_days, b)

Per-late-window confirmed harmonisation offset from the fitted break steps b on the manually specified break_days (grid day-indices). Unlike cumulative_occupancy_offset this is not cumulative. The confirmed likelihood fits between-vintage increments, so a one-off retrospective base integration inflates exactly one increment and every later vintage differences against the new, higher base. Each entry is therefore the step for that window alone, zero on every non-break window.

Added to the modelled confirmed mean of the break window in confirmed_cases_model and confirmed_deaths_model, so the fit can explain an increment that is mostly reattached backlog without inflating the latent incidence that drives Rt. Each b_j is sampled there around the published discrepancy (break_step_centres), symmetric because a harmonisation can revise down as well as up. Returns a length-length(late_days) vector of eltype(b), zeros when there are no break days.

source
BVDOutbreakSize.confirmed_cases_model Method
julia
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,
    lab_daily_history,
    tests_analysed,
    receipt,
    specimen_intensity,
    severity_enrichment,
    sensitivity,
    specificity,
    overdispersion,
    confirmed_break_days,
    confirmed_break_gross,
    confirmed_break_sd,
    cutoff,
    simulated
) -> Any

Laboratory pipeline likelihood. Two streams driven by the shared renewal onsets and the suspected-case pipeline from reported_cases_model:

  • Specimens analysed. The suspected daily pipeline (p_drc-scaled BVD onset-to-report signal plus the non-BVD background λ_bg) is carried through a sampled report-to-analysed delay and thinned by the tested fraction τ_test, giving the expected daily analysed-specimen volume. Its between-vintage increments are fitted against lab_history with a NegativeBinomial sharing the surveillance dispersion k, identifying τ_test and the delay. The same volume is the denominator proxy in the early and unanchored late windows below. The specimens-received series is not modelled: analysed is the throughput that produces confirmed cases, and received overshoots it by the laboratory backlog.

  • Confirmed positives. The confirmed counts are scored as an overdispersed BetaBinomial of the observed specimens-analysed denominator in each laboratory window (confirmed_positivity_windows, safe_betabinomial), with the composition-linked per-window positivity below and a shared intra-window overdispersion (confirmed_overdispersion_model). Conditioning the positives on the observed denominator, rather than on a modelled count scaled by p_drc · s_test · τ_test, removes the multiplicative ascertainment ridge, so the outbreak size is pinned by the deaths and exports streams while the laboratory positivity tracks the noisy per-vintage data. The early confirmed vintages with no per-vintage analysed denominator (18-23 May) are scored as NegativeBinomial counts against the modelled laboratory volume with the same pooled positivity, so all the confirmed data is used.

The per-window positivity ties the tested BVD share to the suspect-pool composition φ_v = (p_drc·BVD)_v / ((p_drc·BVD)_v + λ_bg_v), severity-upsampled by a decaying enrichment δ0, then mapped to the tested-positive probability through the assay sensitivity and specificity, p = s · q + (1 − spec)(1 − q). The false-positive term (1 − spec)(1 − q) makes the confirmed counts respond to the non-BVD share 1 − q, so the laboratory positivity identifies the background λ_bg rather than absorbing it into a free curve.

The tested fraction τ_test and background rate λ_bg come from reported_cases_model so the suspected and laboratory streams share them. Exposes the per-window positivity, the expected analysed and confirmed totals at the cut-off, and the cut-off positivity as derived quantities.

source
BVDOutbreakSize.confirmed_deaths_model Method
julia
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,
    confirmed_break_days,
    confirmed_break_gross,
    confirmed_break_sd,
    receipt_pmf,
    capacity_start,
    case_analysed_daily,
    case_suspected_daily,
    scaling,
    testing,
    sensitivity,
    specificity,
    cutoff,
    simulated
) -> Any

Laboratory-confirmed-deaths likelihood, the death analogue of the confirmed-case laboratory pipeline (confirmed_cases_model). The death side has no published analysed denominator, so the confirmed-death increments are NegBinomial(k) counts of the modelled death volume, the same modelled-volume route the early and post-lab confirmed-case windows use.

Three pieces:

  • 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 · case_analysed_daily · susp_death / susp_case, with susp_death and susp_case the suspected deaths and cases carried to receipt by the confirmed cases' delay. The case volume carries the laboratory capacity onset, so the death volume inherits it. The scaling (death_testing_scaling_model) is the per-suspect testing-intensity difference between deaths and living suspects, and the death-to-case ratio carries the suspect-pool severity and the suspected-death level. The death-only composer has no case stream and falls back to a death testing fraction (death_testing_fraction_model).

  • Death-pool composition. The BVD share of the suspected deaths at receipt, q_death = bvd_death / (bvd_death + bg_death) per day, from the death series' own BVD and background components. The death background, tied to the case background by cfr_bg (see deaths_model and background_cfr_model), keeps q_death below one.

  • Assay positivity. p = s_test · q_death + (1 − spec)(1 − q_death) with PCR sensitivity s_test, named as confirmed_cases_model names it (test_sensitivity_model), and specificity spec (test_specificity_model), the same form as the confirmed-case positivity, drawn from the same priors as separate death-stream parameters.

Returns the cut-off realised death testing intensity τ_death, specimens per suspected death rather than a probability, which exceeds one where the laboratory analyses more specimens than there are suspected deaths. The death-only composer, with no case volume to scale from, falls back to a sampled fraction that is bounded by one (death_testing_scaling_model). Also returns the testing scaling, the cut-off death-pool composition q_death, the confirmation positivity and the expected confirmed-death count.

source
BVDOutbreakSize.confirmed_positivity_windows Function
julia
confirmed_positivity_windows(
    confirmed_history,
    lab_history
) -> Union{NamedTuple{(:obs_days, :obs_positives, :obs_analysed, :early_days, :early_increments, :early_start, :late_days, :late_increments, :late_analysed, :late_start), <:Tuple{Vector{Int64}, Vector{Int64}, Vector{Int64}, Any, Any, Any, Vector{Int64}, Vector{Int64}, Vector{Int64}, Int64}}, NamedTuple{(:obs_days, :obs_positives, :obs_analysed, :early_days, :early_increments, :early_start, :late_days, :late_increments, :late_analysed, :late_start), <:Tuple{Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Any, Vector{Int64}, Vector{Int64}, Vector{Int64}, Any}}}
confirmed_positivity_windows(
    confirmed_history,
    lab_history,
    lab_daily_history
) -> Union{NamedTuple{(:obs_days, :obs_positives, :obs_analysed, :early_days, :early_increments, :early_start, :late_days, :late_increments, :late_analysed, :late_start), <:Tuple{Vector{Int64}, Vector{Int64}, Vector{Int64}, Any, Any, Any, Vector{Int64}, Vector{Int64}, Vector{Int64}, Int64}}, NamedTuple{(:obs_days, :obs_positives, :obs_analysed, :early_days, :early_increments, :early_start, :late_days, :late_increments, :late_analysed, :late_start), <:Tuple{Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Any, Vector{Int64}, Vector{Int64}, Vector{Int64}, Any}}}
confirmed_positivity_windows(
    confirmed_history,
    lab_history,
    lab_daily_history,
    confirmed_break_days
) -> Union{NamedTuple{(:obs_days, :obs_positives, :obs_analysed, :early_days, :early_increments, :early_start, :late_days, :late_increments, :late_analysed, :late_start), <:Tuple{Vector{Int64}, Vector{Int64}, Vector{Int64}, Any, Any, Any, Vector{Int64}, Vector{Int64}, Vector{Int64}, Int64}}, NamedTuple{(:obs_days, :obs_positives, :obs_analysed, :early_days, :early_increments, :early_start, :late_days, :late_increments, :late_analysed, :late_start), <:Tuple{Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Vector{Int64}, Any, Vector{Int64}, Vector{Int64}, Vector{Int64}, Any}}}

Align the confirmed-case counts onto the laboratory windows, splitting them into three non-overlapping groups so all the confirmed data is used:

  • Observed windows, where an analysed-specimen increment is available (the laboratory series only differences cleanly from its second vintage onward. The first cumulative value is the baseline). Each window's analysed increment is the Binomial denominator and the matching confirmed increment the positives. A zero analysed increment (the 24-25 May analysis stall) or a window whose positives exceed its denominator is merged forward, so every observed window has obs_analysed > 0 and 0 ≤ obs_positives ≤ obs_analysed.

  • Early windows, the confirmed vintages up to and including the first laboratory date (18-23 May), which have no per-vintage analysed denominator. Their confirmed increments are returned so the model can score them against the modelled laboratory volume, with the per-window positivity partially pooled with the observed windows.

  • Late windows, the confirmed vintages strictly after the last laboratory date (INSP's confirmed-only format, with no national analysed-specimen denominators). Like the early windows, their confirmed increments are scored against the modelled laboratory volume with the pooled positivity. Their day grid carries a late_start day (the last laboratory day) so the model bins the modelled volume over each late window's own day range, pinned there rather than from day 0, which would double-count the observed window volume. A late day that publishes a 24h analysed count (lab_daily_history) is flagged in late_analysed so the model can score it as a Binomial of that observed denominator, anchoring its positivity to data. A day listed in confirmed_break_days is the exception: its denominator is dropped (late_analysed set to 0) so it stays a modelled-volume window. On such a day INSP has retrospectively integrated a harmonised provincial base, so most of the confirmed increment is reattached backlog rather than positives out of that day's specimens, and pairing the two would score a wildly inflated positivity (22 July 2026: 369 of 414 analysed, 89% implied against the 23% actually reported). The increment is still scored, with the backlog absorbed by a fitted level step in confirmed_cases_model (see confirmed_break_offset).

The three groups partition the confirmed counts at the first and last laboratory dates, so no confirmed case is counted twice. Returns (; obs_days, obs_positives, obs_analysed, early_days, early_increments, late_days, late_increments, late_analysed, late_start) of grid day-indices and per-window counts. late_analysed[i] is the observed 24h analysed denominator for late day i (0 when none was published). The observed and late groups are empty when no laboratory history is present and every confirmed vintage becomes an early window.

early_start and late_start are grid days on every path but two. An empty confirmed_history returns 0 for both, and a confirmed_history with no lab_history returns 0 for late_start alone. Those two are the shortcut returns, and 0 there is a sentinel rather than a day. Elsewhere each is a real day whether or not its own window list is empty, since early_start is the first confirmed vintage and late_start the last laboratory day, neither of which depends on how many windows fell out.

An empty group's start is still read: confirmed_cases_model passes it to bin_increments as a one-element edge list. The single bin that comes back is dropped, so the value cannot reach the likelihood, but a caller reading the sentinel as a grid day would index on 0.

Pure integer bookkeeping on the observed data, so it carries no gradient.

source
BVDOutbreakSize.deaths_model Method
julia
deaths_model(
    deaths_history,
    total_deaths::Union{Missing, Integer},
    onsets::AbstractVector,
    k::Real;
    suspected_daily_deaths_history,
    cfr,
    ascertainment,
    case_bg_daily,
    background_cfr,
    onset_to_death,
    cutoff,
    simulated
) -> Any

DRC suspected-deaths likelihood, per-vintage time series. Convolves the daily onsets with the sampled onset-to-death delay, scales by the CFR and the death ascertainment p_death, and fits the between-vintage increments with a NegativeBinomial sharing the surveillance dispersion k (surveillance_dispersion_model). The death history ends at the cut-off, so the cut-off total is the final increment and is not scored separately. The onset-to-death prior is the convolution of the two atomic line-list components (onset-to-admission and admission-to-death, from the bdbv-linelist-analysis submodule), implied mean ≈ 12.8 d, SD ≈ 7.0 d.

A fatal BVD infection is counted as a suspected death only when ascertained, the death analogue of the suspected-case ascertainment p_drc (see death_ascertainment_model). Its prior is informative and high (centre ≈ 0.9), since a death is more reliably reported than a living suspect.

Non-BVD background suspected deaths are added on top. The joint passes the per-day non-BVD suspected-case background case_bg_daily and a background CFR submodel (background_cfr_model), so the background deaths are cfr_bg times that case background, lagged by the onset-to-death delay. The death background therefore inherits its level and time profile from the identified case background rather than carrying a second free, outbreak-size-degenerate rate. With no case_bg_daily the stream is pure BVD, which is what the deaths-only composer fits.

An optional suspected_daily_deaths_history adds the 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. It continues the suspected-death signal where the cumulative deaths_history stops, and its days fall strictly after that series ends, so the two likelihoods cover disjoint days. Empty by default.

Returns the cut-off expected count, the daily death series (total and the BVD and background components), the onset-to-death PMF, the CFR, the death ascertainment and the background CFR for reuse by confirmed_deaths_model and exports_deaths_model.

source
BVDOutbreakSize.expand_vintage_rate Method
julia
expand_vintage_rate(
    rate::AbstractVector,
    days::AbstractVector{<:Integer},
    n::Integer
) -> Any

Expand a per-vintage rate vector rate onto a length-n daily grid, assigning each day the rate of the vintage window it falls in. The windows are delimited by the ascending day indices days (1-based into the grid). Day t ≤ days[1] takes rate[1], a day in (days[i-1], days[i]] takes rate[i], and any day beyond the last vintage takes the last rate (a flat carry-forward of the final window). When days is empty the whole grid takes rate[1] if present, else zero, so a scalar background is recovered. Pure and AD-transparent. The element type follows rate. Used to put the per-window laboratory positivity on the daily grid (confirmed_cases_model).

source
BVDOutbreakSize.exports_deaths_model Method
julia
exports_deaths_model(
    exports_deaths::Union{Missing, Integer},
    travelled_prevalence::AbstractVector,
    CFR::Real,
    od_pmf::AbstractVector,
    incubation_pmf::AbstractVector;
    export_death_days,
    pre_death_exports,
    cutoff,
    simulated
) -> Any

Deaths-among-exported-cases likelihood, dated per-day. Deaths accrue among the travelled at-risk person-time q · prevalence from exports_model, before the export-case ascertainment p_uganda: a death among an exported case would be reported whether or not the case was ascertained as an import, so the death model is not thinned by p_uganda (unlike the export-case count). The day-t expected export-death increment is the CFR-scaled convolution of that travelled prevalence with the infection→death PMF,

so the cumulative export-death intensity Λ_d(t) is its running sum, the discrete analogue of the integral model's ∫ C(s)·S_det·F_death ds. The onset-to-death PMF od_pmf shared from deaths_model is convolved with the incubation PMF to give the infection→death distribution f_d, keyed to infection like the prevalence.

The observed Uganda export deaths are a dated series export_death_days (grid day-indices, ascending), fitted as an inhomogeneous Poisson exactly as exports_model fits the imports: one Poisson term per death day (the increment of Λ_d between consecutive death-day edges), a pre_death_exports ~ Poisson(Λ_d(δ₁−1)) zero term bounding the first death day δ₁, and the last_offset truncation that stops the clock at the last observed death day. With an empty export_death_days the model falls back to the cut-off cumulative Poisson exports_deaths ~ Poisson(Λ_d(n)).

source
BVDOutbreakSize.exports_model Method
julia
exports_model(
    exported_cases::Union{Missing, Integer},
    infections::AbstractVector,
    p_uganda::Real;
    export_case_days,
    pre_detection_exports,
    incubation_pmf,
    source_population,
    traveller,
    onset_to_detection,
    cutoff,
    simulated
)

Uganda exports likelihood (geographic spread). The exports stream is travel-gated, so the at-risk clock starts at infection: a traveller moves and is exported during incubation (pre-symptomatic) and stays at risk of being exported and detected abroad only until the infection→detection delay has elapsed. The expected detected exports by the cut-off therefore accumulate the per-capita travel rate q = daily_travellers / source_population over the at-risk person-time, not over single-day onset events:

with C(s) the cumulative infections and detected(s) the cumulative infections that have already completed the infection→detection delay by day s. The infection→detection delay is the sampled onset-to-detection delay convolved with the shared incubation PMF, so incubation sits inside it, keyed to infection like C(s). With f_det that delay's PMF, the prevalence is the infections convolved with its survival P(delay > τ) = 1 − Σ_{u ≤ τ} f_det(u). Summing the daily at-risk prevalence is the discrete person-time integral. Summing q · onsets instead would charge each case only a single day of travel risk. The onset-to-detection prior is centred on the Ebola onset-to-hospitalisation delay (mean 5.0 d, SD 4.7 d; WHO Ebola Response Team 2014, NEJM), the delay from symptom onset to detection at a point of entry abroad.

Dated per-day likelihood

The observed Uganda imports are a dated series, export_case_days giving the grid day-index of each detection (sorted ascending), fitted as an inhomogeneous Poisson process rather than a single cumulative count. The cumulative export intensity is Λ(t) = sum(export_prevalence[1:t]), so the per-day expected export increment Λ(t) − Λ(t−1) is the day-t at-risk person-time export_prevalence[t]. The likelihood places one Poisson term per import at its detection day, the increment between consecutive detection-day edges (via bin_increments on the daily prevalence series), with two extra terms:

  • pre_detection_exports ~ Poisson(Λ(d₁−1)) observed at zero, the first-detection timing bound: no export is expected before the earliest detection day d₁. The first import's increment is measured from Λ(d₁−1), so the pre-detection weight and the import increments partition Λ(t_last) and the model still conditions on the same total as a cumulative single-Poisson would.

  • last_offset truncation: the travel-gated export clock stops at the last observed import t_last = day of the most recent detection. Days after the last import carry no informative zero (cross-border movement shifts over the outbreak and the most recent days are right-truncated by reporting lag), so prevalence past t_last does not accrue. With an empty export_case_days the model falls back to the cut-off cumulative Poisson exported_cases ~ Poisson(Λ(T)).

Returns the expected cumulative count at t_last, the per-capita travel rate and the daily at-risk prevalence for reuse by exports_deaths_model.

source
BVDOutbreakSize.onset_ascertainment_model Method
julia
onset_ascertainment_model(
    anchor_series::AbstractVector,
    grid_start::Integer,
    grid_end::Integer;
    beta_prior,
    pooling_prior,
    week
) -> Any

Ascertainment level over the onset-date grid [grid_start, grid_end]: a logit-scale offset and slow random walk on top of a delay-weighted anchor series anchor_series (onset_report_anchor_series),

with alpha given by onset_report_ascertainment. beta_prior = Normal(0, 0.75) is a zero-centred logit-scale departure from the anchor, so the triangle's ascertainment defaults to the anchor series' own level and can depart by roughly a factor of two either way. pooling_prior = truncated(Normal(0, 0.1); lower = 0) is deliberately tight: omega shares the onset-date axis with rt_walk_model, so a flat ascertainment level is the default the data has to argue away from. Weekly knots, linearly interpolated (knot_days/interpolate_knots), first knot pinned at zero, mirroring onset_report_hazard_model's calendar walk exactly.

Returns (; alpha, β, σ_a, z_a, ω), alpha and ω length nt = max(grid_end - grid_start + 1, 1).

source
BVDOutbreakSize.onset_nowcast Method
julia
onset_nowcast(
    observed::Real,
    onsets_u::Real,
    δ::Integer,
    logit_h0::AbstractVector,
    γ::AbstractVector,
    u::Integer,
    grid_start::Integer,
    α::Real;
    until
) -> Any
julia
onset_nowcast(observed, onsets_u, δ, logit_h0, γ, u, grid_start, α; until)

Reported count for onset date u at a later delay than the observed count already printed for it at delay δ:

with n_u the modelled symptom onsets, α the ascertainment level and F the cumulative reported proportion (onset_report_F). until is the target delay δ*, the default nothing targeting the eventual total F(u, δ*) = α. Only the reporting between the two delays is modelled, so the estimate stays anchored on the count already printed and closes on it when the two delays leave nothing outstanding.

Its spread is parameter uncertainty in the onsets and the delay curve. The count it targets is latent, so put it through the stream's own bar measurement error before comparing it with a digitised reading. Pure, top-level, allocation-free.

source
BVDOutbreakSize.onset_report_F Method
julia
onset_report_F(
    δ::Integer,
    logit_h0::AbstractVector,
    γ::AbstractVector,
    u::Integer,
    grid_start::Integer,
    α::Real
) -> Any
julia
onset_report_F(δ, logit_h0, γ, u, grid_start, α)

Cumulative reported proportion F(u, δ) = α * G(u, δ) of onset date u's eventual symptom-onset cases reported within δ days, with α the ascertainment level for onset date u (onset_report_ascertainment) and G the normalised delay CDF (onset_report_G). F(u, D-1) = α exactly, so the hazard's delay shape and its ascertainment level are two separate factors rather than one asymptote. Pure, top-level, allocation-free.

source
BVDOutbreakSize.onset_report_G Method
julia
onset_report_G(
    δ::Integer,
    logit_h0::AbstractVector,
    γ::AbstractVector,
    u::Integer,
    grid_start::Integer
) -> Any
julia
onset_report_G(δ, logit_h0, γ, u, grid_start)

Normalised delay CDF G(u, δ) = cdf(u, δ) / cdf(u, D-1), both terms being onset_report_cdf_extrapolated values off one walk, so G(u, D-1) = 1 by construction and G is a proper delay distribution rather than an asymptote that drifts with the hazard level. Built on the extrapolated cdf (calendar index clamped rather than assumed in-range) because the denominator always reaches D - 1 days ahead of u, which for an onset date within D - 1 days of grid_end runs past γ's fitted support even though every scored δ itself stays in range. The calendar effect is held flat at the walk's nearest edge there. The denominator is guarded with safe_rate, and G(u, δ) = 0 for δ < 0 since the numerator already returns 0.

G(u, D-1) = 1 holds except in one corner: if every h(j, ·) is small enough that each 1 - h rounds to 1, both terms underflow to exactly 0 and G is 0 rather than 1. That needs all D baseline hazards below about 1e-16 at once, roughly 30 prior standard deviations away, and it degrades quietly rather than producing NaN. onset_report_ascertainment clamps against the one consequence that would spread, an anchor of exactly zero. Pure, top-level, allocation-free.

source
BVDOutbreakSize.onset_report_anchor Method
julia
onset_report_anchor(
    logit_h0::AbstractVector,
    γ::AbstractVector,
    u::Integer,
    grid_start::Integer,
    a::AbstractVector
) -> Any
julia
onset_report_anchor(logit_h0, γ, u, grid_start, a)

Delay-weighted average of the calendar-indexed daily ascertainment series a over onset date u's reporting window, anchor(u) = Σ_d g(u, d) * a[clamp(u + d, 1, length(a))], with g(u, d) = G(u, d) - G(u, d-1) (onset_report_G) the normalised delay PMF, d = 0 … D-1. Since Σ_d g(u, d) = 1, anchor(u) lies within [minimum(a), maximum(a)] and a constant a gives back that constant exactly. Both hold wherever G reaches one, so they inherit the underflow corner onset_report_G documents, where the weights sum to zero and the anchor with them. a is indexed on the calendar/report axis and clamped at both ends, which is what reconciles the onset-indexed anchor(u) with a report-indexed series such as the confirmed pipeline's daily ascertainment (see onset_reporting_model). Pure and top-level, over a one-column onset_report_cdf_table.

source
BVDOutbreakSize.onset_report_anchor_series Method
julia
onset_report_anchor_series(
    cdf_table::AbstractMatrix,
    u_lo::Integer,
    a::AbstractVector
) -> Any
julia
onset_report_anchor_series(cdf_table, u_lo, a)

onset_report_anchor for every onset date the delay-CDF table onset_report_cdf_table spans, its column k being onset date u_lo + k - 1. This is the form onset_reporting_model calls, off the table it shares with onset_report_moments. Pure, top-level, one indexed loop over the table and no further hazard evaluation.

source
BVDOutbreakSize.onset_report_anchor_series Method
julia
onset_report_anchor_series(
    logit_h0::AbstractVector,
    γ::AbstractVector,
    grid_start::Integer,
    grid_end::Integer,
    a::AbstractVector
) -> Vector
julia
onset_report_anchor_series(logit_h0, γ, grid_start, grid_end, a)

onset_report_anchor evaluated for every onset date u in grid_start:grid_end, the onset-date grid the ascertainment walk spans (see onset_reporting_model), over a delay-CDF table (onset_report_cdf_table) built here. Returns an empty vector when grid_end < grid_start.

source
BVDOutbreakSize.onset_report_ascertainment Method
julia
onset_report_ascertainment(
    anchor_series::AbstractVector,
    β::Real,
    ω::AbstractVector
) -> Any
julia
onset_report_ascertainment(anchor_series, β, ω)

Ascertainment level α(u) = logistic(logit(anchor(u)) + β + ω(u)) for every onset date u the ascertainment walk spans, anchor_series[k] being the delay-weighted anchor at onset date grid_start + k - 1 (onset_report_anchor_series), β a sampled logit-scale offset and ω the calendar-time random walk over the same onset dates (onset_ascertainment_model). α is a level rather than an asymptote needing D days of walk to become observable, so it is returned over the walk's full span with no restriction.

The anchor is clamped into (0, 1) before the logit, the same guard composition_positivity applies to its own probability. Exactly 0 or 1 is not reachable through the confirmed pipeline, whose positivity is already clamped, but it is reachable if every hazard underflows so that onset_report_anchor's weights sum to zero. logit(0) is -Inf, and while the forward value stays finite the gradient is NaN, which would poison the whole log-density rather than this stream's part of it. Pure, top-level, elementwise broadcast.

source
BVDOutbreakSize.onset_report_cdf Method
julia
onset_report_cdf(
    δ::Integer,
    logit_h0::AbstractVector,
    γ::AbstractVector,
    u::Integer,
    grid_start::Integer
) -> Any
julia
onset_report_cdf(δ, logit_h0, γ, u, grid_start)

Un-normalised reported proportion of onset date u's eventual symptom-onset cases reported within δ days, under the discrete hazard h(d, t) = logistic(logit_h0[d+1] + γ[t - grid_start + 1]) (delay d, report calendar day t = u + d, indexed into the calendar-time random-walk vector γ which starts at grid day grid_start):

This is the delay shape only: it is not the reported proportion F the model scores, which also carries ascertainment (see onset_report_G and onset_report_F). δ < 0 returns exactly 0, the right- truncation case, and is the building block both G and the delay-weighted onset_report_anchor rely on. δ is capped at D-1, so the returned value is constant for every δ >= D-1.

Pure, top-level and allocation-free (a single indexed @inbounds loop, no closures), so it composes AD-transparently into the vectorised moment and total functions below, avoiding the Enzyme boxing a map(...) do or a closure over logit_h0/γ inside an @model body hits.

source
BVDOutbreakSize.onset_report_cdf_extrapolated Method
julia
onset_report_cdf_extrapolated(
    δ::Integer,
    logit_h0::AbstractVector,
    γ::AbstractVector,
    u::Integer,
    grid_start::Integer
) -> Any
julia
onset_report_cdf_extrapolated(δ, logit_h0, γ, u, grid_start)

Like onset_report_cdf, but the calendar-time index into γ is clamped to [1, length(γ)] rather than assumed in-range, so an onset date u outside [grid_start, grid_start + length(γ) - 1] (a calendar day the fitted γ walk has no support for) still returns a well-defined value: the calendar effect at that unseen report day is held flat at the walk's nearest known edge (γ[1] if the whole delay window falls before grid_start, γ[end] if after). A strict extension: it agrees exactly with onset_report_cdf whenever every index touched is in-range.

Used by onset_report_expected_total to extend the cut-off total to onset dates the digitised triangle never covers (days 1:grid_start-1). The increment likelihood in onset_reporting_model never needs it, since every scored onset date is inside [grid_start, grid_end] by construction (grid_start = minimum(onset_days)). Pure, top-level, allocation-free.

source
BVDOutbreakSize.onset_report_cdf_table Method
julia
onset_report_cdf_table(
    logit_h0::AbstractVector,
    γ::AbstractVector,
    grid_start::Integer,
    u_lo::Integer,
    u_hi::Integer
) -> Matrix
julia
onset_report_cdf_table(logit_h0, γ, grid_start, u_lo, u_hi)

onset_report_cdf_extrapolated at every delay d = 0 … D-1 for every onset date u in u_lo:u_hi, as a D × (u_hi - u_lo + 1) matrix: table[d + 1, k] is the value at delay d and onset date u_lo + k - 1. A column is one pass of the survival recurrence, so it holds the normalised delay CDF's numerators alongside its denominator table[D, k]. Returns a D × 0 matrix when u_hi < u_lo.

onset_reporting_model builds one table over its own onset-date grid, which onset_report_anchor_series and onset_report_moments then both read, so the reporting hazard is evaluated once per (delay, onset date) cell for the whole stream rather than once per use. Pure, top-level, single indexed loop (see onset_report_cdf for the AD-safety rationale).

source
BVDOutbreakSize.onset_report_expected_total Method
julia
onset_report_expected_total(
    onsets::AbstractVector,
    logit_h0::AbstractVector,
    γ::AbstractVector,
    grid_start::Integer,
    alpha::AbstractVector,
    as_of::Integer
) -> Any
julia
onset_report_expected_total(onsets, logit_h0, γ, grid_start, alpha, as_of)

Expected reported symptom-onset total as of grid day as_of, Σ_u onsets[u] · F(u, as_of - u) for u in 1:as_of (clamped to 1:length(onsets)), the onset-stream analogue of every other stream's expected_*_T cut-off deterministic. Callers pass the model cut-off n for as_of, so the total is directly comparable to expected_confirmed_T/expected_deaths_T, which sum their full 1:n series. Passing the triangle's own last report day instead would give a total anchored a few days earlier than every sibling.

γ and alpha both span only the digitised triangle's own grid (see onset_reporting_model), which starts after grid day 1 and ends at or before as_of, so both the oldest and the most recent terms need a calendar effect and an ascertainment level the fit has no estimate for. onset_report_F already holds both flat at their nearest fitted edge (see onset_report_G), so no separate extrapolated form is needed here.

Safe for any as_of and any γ/alpha length, including the degenerate length(γ) < D case, because both indices are clamped rather than assumed in range. An empty delay support (D = 0) gives zero. Pure, top-level, single indexed loop.

source
BVDOutbreakSize.onset_report_hazard_model Method
julia
onset_report_hazard_model(
    grid_start::Integer,
    grid_end::Integer;
    D,
    baseline_prior,
    pooling_prior,
    walk_sigma_prior,
    week,
    walk_start,
    basis
) -> Any

Discrete symptom-onset reporting-delay hazard, nonparametric over the delay and drifting over calendar time. Two non-centred random effects:

  • a baseline logit hazard over the delay dimension d = 0 … D-1, a partially-pooled non-centred random effect over delay whose deviations sum to zero, so η0 is the mean logit hazard:

    with Q the sum-to-zero basis (sum_to_zero_basis).

  • a calendar-time random walk on report date, weekly knots linearly interpolated to the daily grid (rt_walk_model's non-centred cumulative-sum walk, same construction):

The walk is indexed on the report-date grid [grid_start, grid_end], not the onset/infection-date axis rt_walk_model already carries a smooth calendar-time effect on. Indexing this walk on onset date would leave it unidentified against the reproduction-number walk, both free to explain the same onset-date-indexed rise and fall. Report date lags onset date by the delay itself (mean ≈6 d) plus the ≈6 d incubation period the onset series is already convolved from, so the two walks act on different, if overlapping, calendar windows.

The walk holds at zero on every day before walk_start and only moves from there. It changes the reporting delay, and a delay is seen only between snapshots, so before the first snapshot a calendar shift cannot be told apart from the baseline hazard. onset_reporting_model sets walk_start to the first snapshot's report day. The grid, and so the indexing of γ, is unchanged, so every consumer reads it as before. The default walk_start = grid_start starts the walk on the grid's first day.

The default baseline_prior = Normal(logit(0.13), 0.7) targets a median onset-to-report delay of roughly 5 days under a constant-hazard approximation ((1 - 0.13)^5 ≈ 0.5), close to the ≈6-day median a line-list re-analysis found. The wide SD leaves the per-delay random effect free to depart substantially, so the fitted delay is correspondingly uncertain (that analysis put its 7-day reporting fraction at 43-68%, wide because the digitisation noise is comparable to the between-vintage increments the estimate rests on).

walk_sigma_prior is a half-normal with SD 0.3, sized so a drift the scored triangle can show is reachable rather than extreme. A 14% shift in a snapshot's printed level needs a calendar shift of roughly 0.32 at delay 8 and 0.58 at delay 12 on the logit hazard, holding the baseline at its prior median, and after the six weekly knots the scored window spans γ has SD E[σ_γ]·√6 = 0.586, putting that shift at 0.5-1σ. It stays concentrated at zero, so a flat reporting profile is the default the data has to argue away from. Widening it further risks the walk explaining recent onset-date structure that belongs to R_t, so any change here should report R_t over the final fortnight and C_T either side.

Returns (; logit_h0, γ, grid_start, η0, σ_h0, σ_γ), with γ length max(grid_end - grid_start + 1, 1), indexed from grid_start.

source
BVDOutbreakSize.onset_report_moments Method
julia
onset_report_moments(
    cdf_table::AbstractMatrix,
    u_lo::Integer,
    onsets::AbstractVector,
    grid_start::Integer,
    alpha::AbstractVector,
    onset_idx::AbstractVector{<:Integer},
    cur_report_idx::AbstractVector{<:Integer},
    prev_report_idx::AbstractVector{<:Integer}
) -> NamedTuple{(:means, :level_cur, :level_prev), <:Tuple{Any, Any, Any}}
julia
onset_report_moments(cdf_table, u_lo, onsets, grid_start, alpha,
    onset_idx, cur_report_idx, prev_report_idx)

onset_report_moments off a delay-CDF table (onset_report_cdf_table) whose column k is onset date u_lo + k - 1, which is the form onset_reporting_model calls off the table it shares with onset_report_anchor_series. The table must span every onset_idx; the method above builds one over extrema(onset_idx), and the model's own grid starts at minimum(onset_days) and ends at or after maximum(report_days). Pure, top-level, one indexed loop over the table and no further hazard evaluation.

source
BVDOutbreakSize.onset_report_moments Method
julia
onset_report_moments(
    onsets::AbstractVector,
    logit_h0::AbstractVector,
    γ::AbstractVector,
    grid_start::Integer,
    alpha::AbstractVector,
    onset_idx::AbstractVector{<:Integer},
    cur_report_idx::AbstractVector{<:Integer},
    prev_report_idx::AbstractVector{<:Integer}
) -> NamedTuple{(:means, :level_cur, :level_prev), <:Tuple{Any, Any, Any}}
julia
onset_report_moments(onsets, logit_h0, γ, grid_start, alpha, onset_idx,
    cur_report_idx, prev_report_idx)

Per-cell modelled increment mean and the two modelled cumulative levels it is the difference of, for a batch of (onset day, current report day, previous report day) triples (the increment cells load_onset_curve builds). For cell i,

with δ = report_idx - u_i and onset_report_F supplying F, alpha indexed at onset date u_i (clamped into 1:length(alpha), the grid_start:grid_end ascertainment grid) for its α argument (so δ_prev < 0 at the sentinel prev_report_idx = 0, the virtual empty predecessor of an onset date's first print, contributes ℓ_prev = 0 with no special-casing). An onset_idx outside 1:length(onsets) contributes a zero rate rather than indexing out of bounds. Returns (; means, level_cur, level_prev), each a length-length(onset_idx) vector, with level_cur/level_prev the two cumulative levels each increment differences. The hazard enters through a delay-CDF table (onset_report_cdf_table) built here over extrema(onset_idx). Pure and top-level (see onset_report_cdf for the AD-safety rationale).

source
BVDOutbreakSize.onset_report_scale Method
julia
onset_report_scale(μ::Real, τ::Real, reads::Integer) -> Any
julia
onset_report_scale(μ, τ, reads)

Scalar form of onset_report_scales's per-cell formula, for one increment mean μ over reads digitised bars with the read SD τ. The vector method calls this, so the two cannot drift apart. The onset forecast (onset_forecast_model) calls it directly to give a future reporting increment the same observation scale the likelihood gives a scored cell. See onset_report_scales for what each term means.

source
BVDOutbreakSize.onset_report_scales Method
julia
onset_report_scales(
    means::AbstractVector,
    τ::Real,
    reads::AbstractVector{<:Integer}
) -> Any
julia
onset_report_scales(means, τ, reads)

Per-cell observation scale for the reporting-triangle increment likelihood, the square root of a variance built from three sources.

  • Counting variation of the cases the cell actually reports. The cell is a count of newly reported cases, so it carries its own sampling variation of about its mean, means[i], on top of any reading error. This term is what makes the scale correct for the very first snapshot's cells, which are differenced against an empty predecessor and so score a level rather than a correction (see load_onset_curve): a level of 40 cases has counting variation of about sqrt(40) ≈ 6, far larger than the reading error below, and scoring it on reading error alone would let 28 level cells dominate the joint likelihood.

  • Rounding of each read. A digitised bar is an integer, so every read carries the 1/12 variance of rounding to the nearest count. The term is structural rather than measured, and it is what keeps the read SD below from collapsing to zero on the settled cells whose residual is exactly zero.

  • Read error on each digitised bar, the fitted read SD τ of onset_reporting_model.

reads[i] is the number of bars cell i differences: 1 for a level cell off the first scored snapshot and 2 for a correction between two snapshots. The read variances add, so

with μ_i = means[i] the modelled increment and r_i = reads[i]. The counting term cancels for a genuine correction between two snapshots only to the extent that the two reads share the same realised cases: the newly reported cases in between are a fresh count, and μ_i is exactly their expected number, so the same formula covers both cell kinds without a branch.

The magnitude entering the scale is the modelled increment (means from onset_report_moments), never the raw observed count: feeding the likelihood's own noisy observation back into its variance would bias towards overconfidence on cells that happen to undershoot. Pure, top-level, single indexed loop.

The counting term is Poisson-like, with no separate overdispersion parameter.

source
BVDOutbreakSize.onset_reporting_model Method
julia
onset_reporting_model(
    onset_curve_history,
    onsets::AbstractVector;
    hazard,
    ascertainment,
    anchor,
    D,
    read_sd_prior,
    ν
) -> Any

Symptom-onset reporting-triangle observation model: fits the between-vintage increments of the digitised reporting triangle (load_onset_curve) against the shared latent daily onset series onsets, through a nonparametric delay hazard (onset_report_hazard_model) and an explicit ascertainment level (onset_ascertainment_model). See onset_report_F for how the two combine and onset_report_cdf for how right truncation enters.

Close in spirit to how the R package baselinenowcast treats a reporting triangle: a triangle of between-vintage increments, not a single total column. It differs from reported_cases_model and every other stream here, and from EpiNow2, which nowcast a single evolving total. The same right-truncation mechanism (F(u, δ) = 0 for δ < 0) is applied to corrections between snapshots rather than to a level, so a case already scored in an earlier snapshot is never counted again.

What the triangle separates. Four time-varying multiplicative objects act on or near the same latent series: the reproduction-number walk (rt_walk_model) on the infection/onset axis, this stream's calendar walk γ on the report axis, the delay shape's baseline logit_h0, and the ascertainment walk ω on the onset axis. A change in the onset series moves a column of scored cells, a change in γ moves a row, and a change in logit_h0 moves a diagonal band. Column, row and band are distinguishable once there is more than one snapshot, which is the structural reason this stream is worth fitting, and the reason γ is indexed on the report day rather than the onset day.

Ascertainment is anchored on the confirmed pipeline's own daily ascertainment (p_drc · τ_test · p_pos_grid, delay-weighted onto the onset axis by onset_report_anchor) rather than left free: bvd_joint passes that series in as anchor, onsets_only_model falls back to a constant 0.15 anchor (no confirmed pipeline to borrow from). β and ω let the fitted level depart from that anchor, with ω's tight prior making a flat departure the default the data has to argue away from (see onset_ascertainment_model).

Three things stay genuinely weak. First, logit_h0 and alpha are pinned by levels, not by corrections, since corrections constrain only differences of F. What breaks the tie is each onset date's level cell, its first print differenced against an empty predecessor (see load_onset_curve), together with the onset series being pinned by the other streams. A single-stream onsets_only_model fit has neither, so its ascertainment and C_T stay close to prior-driven. Second, the hazard at the shortest delays is barely observed, since published figures stop short of the report date, so those hazards rest on partial pooling to η0 rather than data. Third, a falling ascertainment and a slowing delay shape both suppress recent bars within one snapshot. Right truncation self-corrects for the delay explanation but not a genuine ascertainment fall, so the vintage structure separates them only as well as the under-ten snapshots covering the affected dates allow.

Whether γ, ω and rt_walk_model pull on each other in practice has not been checked. All three act on overlapping calendar windows and are least constrained over the final fortnight, so any change here should report R_t over the final fortnight and C_T either side.

The alive/dead split the raw triangle carries is not modelled separately: only confirmed_total is fitted, since the confirmed-death stream already carries that split from other data.

onset_curve_history is the load_onset_curve return shape (; onset_days, report_days, prev_report_days, increments). The default empty history makes every loop here a no-op, the degrade-gracefully path for a missing input file. increments may be missing to sample instead of condition (the predictive-generator path).

The observation scale (onset_report_scales) is built from counting variation, the rounding variance of each integer read and a fitted read SD τ ~ read_sd_prior, one for every read of a digitised bar: a correction cell carries two reads' error and a first-snapshot level cell one read's. One count is about 2.9 pixels on the published figures, so a read is a rounding plus an outline pixel, of order one count. The prior is centred on that scale and the data set the value.

The likelihood is Student-t with fixed degrees of freedom ν (default 4). With only a few hundred cells ν is weakly identified, so it is not sampled, as in lab_delay_model. The heavy tail lets the frequently negative measured increments score as large-but-plausible residuals rather than breaking a count likelihood.

Returns (; increments, modelled, scales, logit_h0, γ, grid_start, grid_end, alpha, τ, η0, σ_h0, σ_γ, β, σ_a, ν) with modelled the per-cell increment means the likelihood scores, scales the per-cell observation scales it scores them with, grid_end the report-date grid day the calendar walk was built up to (max(report_days), or grid_start when the history is empty), and the hyperparameters re-exposed at this level for the pairs-plot summary.

source
BVDOutbreakSize.province_composition_model Method
julia
province_composition_model(
    obs_increments::Union{Missing, AbstractMatrix{<:Integer}},
    modelled_confirmed::AbstractMatrix;
    rho_prior,
    ascertainment_sd_prior,
    severity_sd_prior,
    ascertainment_offset_prior,
    basis
) -> Any

Per-province composition of the confirmed cases, from the Tableau 1 spatial tables.

Why a composition and not per-province counts

The per-province confirmed totals are an exact partition of the national ones: at every one of the 17 shared vintages, Ituri + Nord-Kivu + Sud-Kivu equals the national confirmed count confirmed_cases_model already scores. Fitting the per-province counts with their own count likelihood would put the same observations into the joint density twice.

The information the spatial tables add is purely spatial, how the confirmed cases divide between provinces, so the likelihood factorises as

where the national confirmed stream supplies the total term P(N) and this model supplies only the conditional composition term. The vintage totals are conditioned on, never scored, so nothing is counted twice.

Likelihood

The composition is scored by stick-breaking over the patches (stick_breaking_loglik) with the overdispersed Binomial safe_betabinomial, the sequential form of a Dirichlet-multinomial. For each vintage the observed count in patch p is drawn from the cases not yet allocated to patches 1 … p−1, at the conditional share implied by the modelled per-patch confirmed increments:

with the last patch taking the remainder. ρ is an overdispersion shared across provinces and vintages, absorbing the extra-Binomial variation in how cases are attributed to provinces (reporting lags between the provincial and national tables, reassignment of cases between health zones). ρ → 0 recovers a plain Multinomial split.

Arguments

  • obs_increments: (n_patches × n_vintages) observed per-province new-confirmed counts, or missing for the predictive path.

  • modelled_confirmed: (n_patches × n_vintages) modelled expected per-province confirmed increments, binned to the same vintages.

  • rho_prior: prior on the composition overdispersion.

  • ascertainment_sd_prior: prior on the pooling scale of the per-patch ascertainment contrast. nothing turns the contrast off, so the shares are the modelled split alone; the laboratory composition passes that.

severity_sd_prior optionally adds a second sum-to-zero log multiplier on the shares, partially pooled toward equality. The death composition uses it for the per-province case-fatality ratio; the case composition leaves it nothing and samples nothing. A composition identifies only the product of the two multipliers, so their priors are what separates them.

Returns (; shares, rho, obs_increments, province_ascertainment, ascertainment_sd, province_severity, severity_sd) where shares[p, i] is the modelled expected share of patch p at vintage i.

source
BVDOutbreakSize.recovered_model Method
julia
recovered_model(
    recovered_history,
    recovered_total::Union{Missing, Integer},
    confirmed_daily::AbstractVector,
    CFR::Real;
    recovery,
    dispersion,
    confirmation_to_recovery,
    k_external,
    cutoff,
    simulated
) -> Any

DRC recovered-among-confirmed likelihood ("cumul guéris"), an incidence (scaled-convolution) stream: the renewal analogue of the convolution-and- scaling secondary-observation model of EpiNow2 (Abbott et al., 2020). Recoveries are the survivors among laboratory-confirmed cases: the modelled daily confirmed-case incidence confirmed_daily (from confirmed_cases_model) is scaled by the recovery probability p_recover (the confirmed-case survival fraction, recovery_probability_model) and convolved with a sampled confirmation-to-recovery delay,

The recovery fraction p_recover is grounded on the case-fatality ratio CFR, since a recovered case is one that did not die, adjusted for the confirmed population by a sampled log-odds offset (see recovery_probability_model). A case is taken to be confirmed before it is recorded as recovered, so the recovery follows confirmation by the confirmation-to-recovery delay. A positive result could in principle return after a patient has already recovered and been discharged, but "cumul guéris" is assumed to count confirmed cases carefully recorded as recovered.

The cumulative recovered series ends at the cut-off, so its between-vintage increments are fitted with a NegativeBinomial whose dispersion is sampled here rather than shared with the other streams. The convolution right-censors recoveries that have not yet resolved by the cut-off, so a small observed recovered count is consistent with a high eventual survival fraction and a long recovery delay. Empty by default. A missing cut-off total leaves the increments missing (the predictive-generator path). Returns the recovery probability, the recovery-delay mean, the dispersion, the daily recovered series and the cut-off total.

source
BVDOutbreakSize.reported_cases_model Method
julia
reported_cases_model(
    reported_history,
    reported_cases::Union{Missing, Integer},
    onsets::AbstractVector,
    k::Real,
    p_drc::Real;
    suspected_daily_history,
    positivity,
    background_re,
    onset_to_report,
    cutoff,
    simulated
) -> Any

DRC suspected (reported) cases likelihood, per-vintage time series. The expected daily suspected cases are a BVD-driven onset-to-report convolution scaled by the DRC ascertainment fraction p_drc, plus an additive non-BVD background rate λ_bg per day (so a suspected case need not be a true BVD infection). Reads the modelled cumulative suspected cases at each vintage day off the daily series and fits the increments with a NegativeBinomial sharing k.

An optional suspected_daily_history adds the daily new-suspect inflow ("nouveaux cas suspects du jour"), per-day counts scored against the modelled daily suspected series at each report day. It continues the suspected signal where the cumulative reported_history stops, once INSP began reclassifying suspects downward and the cumulative total fell. It shares the suspect pipeline and dispersion with the cumulative stream, and its days fall strictly after that series ends, so the two likelihoods cover disjoint days. Empty by default.

The background and testing fraction are sampled by an injected test_positivity_model, and the onset-to-report delay is injected, defaulting to a weakly-informative prior on the onset-to-notification delay (mean 4.5 d, SD 3.6 d), consistent with Ebola surveillance reporting delays.

Returns the onset-to-report PMF and the BVD onset-to-report daily series (unit ascertainment, no background) so confirmed_cases_model can reuse the same report kernel, the sampled background rate and testing fraction, and the implied per-suspected positivity (the BVD share of the expected suspected total) as a derived quantity for comparison with the sitrep.

source
BVDOutbreakSize.treatment_flow_model Method
julia
treatment_flow_model(
    isolation_history,
    bvd_reports_daily::AbstractVector,
    bg_daily::AbstractVector,
    p_drc::Real,
    CFR::Real;
    capacity_history,
    admissions_history,
    deaths_history,
    ruleout_history,
    absconded_history,
    confirmed_incare_history,
    suspect_incare_history,
    conf_hazard_daily,
    defaults,
    admission,
    severity,
    capacity,
    dispersion,
    k_external,
    cfr_modifier_prior,
    abscond_prior,
    incare_confirm_log_prior,
    admission_delay,
    death_los,
    recovery_los,
    ruleout_los,
    occupancy_break_days,
    occupancy_break_sd,
    bvd_reports_matrix,
    background_split,
    patch_ascertainment,
    province_isolation,
    province_capacity,
    province_admissions,
    patch_capacity,
    province_split_rho_prior,
    cutoff,
    simulated
) -> Any

DRC treatment-centre patient-flow likelihood. Occupancy is built as a running balance of latent admission, discharge and abscond events rather than a convolution. The "Patients en isolement" figure is the daily occupied-bed count, fitted as latent bed demand right-censored at the implied-capacity bound.

Admissions are the reported suspects (reported_cases_model) carried through a short suspected→admission delay and split into a BVD true-case inflow A_bvd = p_iso_bvd · p_drc · bvd_reports at the severity-skewed rate p_iso_bvd (isolation_severity_model) and a non-BVD inflow A_bg = p_iso · bg_daily at the base rate (isolation_admission_model).

The clinical layer sets the total occupancy and discharge flows and is label-independent, so a true case can die before confirmation. A BVD true case leaves by death (weight CFR_iso) on the admission→death stay or recovery (weight 1 − CFR_iso) on the admission→recovery stay. A non-BVD admission rules out on the rule-out stay. Absconds drain the suspect pool at κ · O_susp(t−1). The total latent demand is the running balance

with A = A_bvd + A_bg. The in-care deaths flow is suspect and confirmed combined (the Tableau 6 décédés row), scored directly and not gated by confirmation. The in-care fatality CFR_iso = logistic(logit(CFR) + β_iso) adjusts the infection CFR to the admitted population by a sampled modifier β_iso. It is conditional on admission, not a causal treatment effect.

A label overlay carves the occupied stock into confirmed and suspect sub-stocks without removing anyone from the total. Confirmation relabels an occupied true case at the daily hazard ρ · τ_test · p_pos, the community hazard borrowed from the lab pipeline (confirmed_cases_model) scaled by a sampled in-care modifier ρ = exp(γ_conf) the confirmed/suspect census identifies. The suspect sub-stock is the remainder O_susp(t) = D(t) − O_conf(t), so a case that dies before confirmation is a suspect death in the combined deaths flow.

Capacity enters only as a fixed, data-derived censoring bound, the latent demand staying uncapped. The occupancy likelihood is a NegativeBinomial around the demand, right-censored at the recorded implied-capacity series (censoring_cap, censored_occupancy_model). The capacity walk C(t) (bed_capacity_walk_model) carries the implied-capacity likelihood. It is an ingredient of the beds, not the beds reported: the cut-off beds, occupancy and shortfall are built by patch and summed (cutoff_occupancy), each patch's beds its modelled capacity floored at its recorded beds.

The daily in-care outcome flows (deaths, rule-outs, admissions, absconds) are each an optional NegativeBinomial stream scored against the Tableau 6 patient-movement counts, a no-op when their history is empty. The two sub-stocks are scored against the Tableau 6 census breakdown (dont confirmés / dont suspects) where published, in place of the total on those days. An opt-in occupancy offset Δ(t) on the supplied occupancy_break_days (cumulative_occupancy_offset) absorbs a between-report measurement-basis discontinuity in the isolation series. Empty (the default) is a no-op.

With per-patch reports (bvd_reports_matrix) and the province rows province_isolation, province_capacity and province_admissions, the province occupancy, bed and 24h admission counts are each scored as a split of the printed sum of the provinces present that day (province_split_model). The occupancy split is on the uncapped per-patch demand, the national demand shared out by each patch's stock of admissions through the stays. The bed split is on each patch's daily share of the national walk, centred on its cumulative admissions (patch_capacity_share_model). The admission split is on each patch's modelled admissions. The national likelihoods are kept, so the splits add only the spatial signal.

Exposes the cut-off occupancy, bed demand and shortfall, the utilisation, the BVD share of demand, CFR_iso and β_iso, the length-of-stay, the two sub-stock prevalences and their stay means, the in-care fraction and the daily series for forecasting and replication.

source
BVDOutbreakSize.BetaBinomialVector Type
julia
struct BetaBinomialVector{I<:(AbstractVector{<:Integer}), P<:(AbstractVector{<:Real}), R<:Real} <: Distributions.Distribution{Distributions.Multivariate, Distributions.Discrete}

Independent safe_betabinomial counts, entry i out of trials[i] with mean probability p[i] and the shared overdispersion ρ. logpdf sums the entries' terms in one loop.

Fields

  • trials::AbstractVector{<:Integer}: Per-entry trial counts.

  • p::AbstractVector{<:Real}: Per-entry mean probabilities.

  • ρ::Real: Shared overdispersion.

source
BVDOutbreakSize.CensoredNegBinomialVector Type
julia
struct CensoredNegBinomialVector{T<:Real, V<:(AbstractVector{<:Real}), U<:(AbstractVector{<:Real})} <: Distributions.Distribution{Distributions.Multivariate, Distributions.Discrete}

Independent right-censored safe_nbinomial counts, entry i about μ[i] and censored at upper[i], both through safe_rate, with the shared dispersion k. Built by censored(NegBinomialVector(k, μ); upper). A count at its ceiling scores the censored tail nbinomial_logtail.

Fields

  • k::Real: Shared dispersion.

  • μ::AbstractVector{<:Real}: Per-entry means.

  • upper::AbstractVector{<:Real}: Per-entry censoring ceilings.

source
BVDOutbreakSize.NegBinomialVector Type
julia
struct NegBinomialVector{T<:Real, V<:(AbstractVector{<:Real})} <: Distributions.Distribution{Distributions.Multivariate, Distributions.Discrete}

Independent safe_nbinomial counts, entry i about the safe_rate of μ[i] with the shared dispersion k. logpdf sums the entries' terms in one loop.

Fields

  • k::Real: Shared dispersion.

  • μ::AbstractVector{<:Real}: Per-entry means.

source
BVDOutbreakSize.SafeNegBinomial Type
julia
struct SafeNegBinomial{T<:Real, M<:Real} <: Distributions.Distribution{Distributions.Univariate, Distributions.Discrete}
julia
SafeNegBinomial(k, μ)

safe_nbinomial(k, μ) whose draws saturate at typemax(Int) rather than throwing InexactError past it, as SafePoisson does. logpdf is the safe_nbinomial one.

Fields

  • k::Real: Dispersion.

  • μ::Real: Mean.

source
BVDOutbreakSize.SafePoisson Type
julia
struct SafePoisson{T<:Real} <: Distributions.Distribution{Distributions.Univariate, Distributions.Discrete}
julia
SafePoisson(λ)

Poisson with mean λ whose draws saturate at typemax(Int) rather than throwing InexactError past it. logpdf is the Poisson one. The draw follows SafePoisson in ComposableTuringIDModels (https://github.com/EpiAware/ComposableTuringIDModels.jl, src/utils/SafePoisson.jl).

Fields

  • λ::Real
source
BVDOutbreakSize.SplitCountVector Type
julia
struct SplitCountVector{B<:BVDOutbreakSize.BetaBinomialVector, N<:BVDOutbreakSize.NegBinomialVector} <: Distributions.Distribution{Distributions.Multivariate, Distributions.Discrete}

Independent counts split on their trial counts: entry i with trials[i] > 0 is a safe_betabinomial of trials[i] with mean probability p[i] and overdispersion ρ, and every other entry a safe_nbinomial about μ[i] with dispersion k. logpdf is a BetaBinomialVector over the first group plus a NegBinomialVector over the second.

Fields

  • anchored::Vector{Int64}: Entries with a trial count.

  • unanchored::Vector{Int64}: Entries without one.

  • binomial::BVDOutbreakSize.BetaBinomialVector: Distribution of the anchored entries.

  • negbinomial::BVDOutbreakSize.NegBinomialVector: Distribution of the unanchored entries.

source
BVDOutbreakSize.StudentTVector Type
julia
struct StudentTVector{V<:(AbstractVector{<:Real}), S<:(AbstractVector{<:Real}), T<:Real} <: Distributions.Distribution{Distributions.Multivariate, Distributions.Continuous}

Independent safe_studentt values, entry i about μ[i] with scale σ[i] and the shared degrees of freedom ν. logpdf sums the cells' terms in one loop, with the normalising constant evaluated once.

Fields

  • μ::AbstractVector{<:Real}: Per-entry locations.

  • σ::AbstractVector{<:Real}: Per-entry scales.

  • ν::Real: Shared degrees of freedom.

source
BVDOutbreakSize.safe_studentt Method
julia
safe_studentt(μ::Real, σ::Real, ν::Real) -> Any
julia
safe_studentt(μ, σ, ν)

NaN / Inf-safe location-scale Student-t distribution μ + σ · Tν, built from Distributions.TDist(ν) via the affine-combination operators. σ is floored away from zero and non-finite values, mirroring safe_nbinomial's domain guard. A non-positive or non-finite ν falls back to 4, the caller's own default, rather than to the smallest value TDist accepts: TDist(1) is Cauchy, so a floor at the domain edge would turn a bad degrees-of-freedom argument into a likelihood with no mean or variance.

Used by onset_reporting_model to score the reporting-triangle increments, which are frequently negative (a later scan reads fewer cases at some onset date than an earlier one, from digitisation noise rather than a real reporting reversal) and so cannot take a count distribution.

source