Real-time nowcasting: correcting right-truncation
In real time, the most recent days of a surveillance series are incomplete: a case with reference day
This case study takes a real, fully-reported case series — the daily confirmed COVID-19 cases from Italy's first wave used in the delays example — treats it as the eventual totals, and truncates its recent tail to mimic the real-time snapshot an analyst would have seen mid-outbreak. It then contrasts two fits: a naive one that treats the truncated counts as complete, and one that wraps the observation model in RightTruncate to scale each reference day's expected count by the fraction of its eventual total reported so far. The correction is the EpiNow2-style CDF-scaling nowcast [3], expressed here as a one-line observation modifier.
The idea
The infection pipeline produces RightTruncate conditions the observation error on that down-weighted mean, so the model's
A naive fit drops the
The full-data model
We build the same composed renewal model as the renewal case study: an autoregressive Renewal infection process, observed with a NegativeBinomialError.
using ComposableTuringIDModels, Distributions, Random, Turing
using CSV, DataFrames
Random.seed!(20240625)
latent = AR(
damp_priors = [truncated(Normal(0.8, 0.05), 0, 1)],
init_priors = [Normal(0.0, 0.25)],
ϵ_t = HierarchicalNormal(std_prior = HalfNormal(0.1)))
renewal = Renewal(gen_distribution = Gamma(1.4, 1 / 0.38);
rt = latent, initialisation_prior = Normal(log(1.0), 1.0))
error = NegativeBinomialError(cluster_factor_prior = HalfNormal(0.1))Take a real series and truncate it
We use the fully-reported Italy confirmed-case series as the eventual totals, then impose a reporting delay and truncate at ReportingCDF builds the completeness curve
datapath = joinpath(pkgdir(ComposableTuringIDModels),
"docs", "src", "case-studies", "data", "italy_data.csv")
italy = CSV.read(datapath, DataFrame)
n = 50
eventual = italy.confirm[1:n] # the eventual (complete) totals
reporting_delay = LogNormal(1.6, 0.5) # mean ≈ 5.6 days
cdf_curve = ReportingCDF(reporting_delay)Now form the right-truncated snapshot: thin each reference day's eventual total to the fraction reported by
completeness = as_turing_model(cdf_curve, n)() # F by age, F[1] = age 0
scale = reverse(completeness) # by reference day t = 1..n
observed_so_far = @. rand(Binomial(eventual, scale))
(complete_tail = eventual[(end - 4):end], truncated_tail = observed_so_far[(end - 4):end])(complete_tail = [3599, 3039, 3836, 4204, 3951], truncated_tail = [1497, 752, 355, 64, 0])The last few days are visibly thinned: the most recent day shows only a fraction of its eventual count.
Fit 1 — naive (no truncation correction)
Condition the plain renewal model on the truncated counts as though they were complete.
naive_model = IDModel(renewal, error)
naive_post = as_turing_model(naive_model, observed_so_far, n)
naive_chain = sample(
naive_post, NUTS(0.9), MCMCThreads(), 500, 2; progress = false)┌ Warning: Only a single thread available: MCMC chains are not sampled in parallel
└ @ AbstractMCMC ~/.julia/packages/AbstractMCMC/C1aKp/src/sample.jl:544
┌ Info: Found initial step size
└ ϵ = 0.00625
┌ Warning: There were 29 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/WteH7/src/mcmc/hmc.jl:483
┌ Info: Found initial step size
└ ϵ = 0.003125
┌ Warning: There were 16 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/WteH7/src/mcmc/hmc.jl:483
┌ Warning: There were 45 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/WteH7/src/mcmc/hmc.jl:483Fit 2 — corrected with RightTruncate
Wrap the same error model in RightTruncate with the same reporting delay. Nothing else changes — the infection process and its latent
corrected_obs = RightTruncate(error, reporting_delay)
corrected_model = IDModel(renewal, corrected_obs)
corrected_post = as_turing_model(corrected_model, observed_so_far, n)
corrected_chain = sample(
corrected_post, NUTS(0.9), MCMCThreads(), 500, 2; progress = false)┌ Warning: Only a single thread available: MCMC chains are not sampled in parallel
└ @ AbstractMCMC ~/.julia/packages/AbstractMCMC/C1aKp/src/sample.jl:544
┌ Info: Found initial step size
└ ϵ = 0.0015625
┌ Warning: There were 6 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/WteH7/src/mcmc/hmc.jl:483
┌ Info: Found initial step size
└ ϵ = 0.0125
┌ Warning: There were 10 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/WteH7/src/mcmc/hmc.jl:483
┌ Warning: There were 16 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/WteH7/src/mcmc/hmc.jl:483The right-truncation bias and its correction
The reproduction number returned (re-exported by Turing) recovers the latent RightTruncate fit removes that bias.
using Statistics
function recent_Rt(posterior, chain; window = 7)
# Per-draw R_t = exp(Z_t); average over draws, then over the recent window.
gens = vec(returned(posterior, chain))
Rt_mean = mean(exp.(g.Z_t) for g in gens)
mean(Rt_mean[(end - window + 1):end])
end
complete_post = as_turing_model(naive_model, eventual, n)
complete_chain = sample(
complete_post, NUTS(0.9), MCMCThreads(), 500, 2; progress = false)
R_complete_recent = recent_Rt(complete_post, complete_chain)
R_naive_recent = recent_Rt(naive_post, naive_chain)
R_corrected_recent = recent_Rt(corrected_post, corrected_chain)
(complete = round(R_complete_recent, digits = 2),
naive = round(R_naive_recent, digits = 2),
corrected = round(R_corrected_recent, digits = 2))(complete = 0.9, naive = 0.32, corrected = 0.95)The naive recent-RightTruncate fit, which knows the recent days are incomplete, pulls the recent
Prior versus posterior
Sampling the corrected model with Prior gives a prior draw over the shared renewal parameters. Overlaying it on the posterior with PairPlots.jl confirms the truncation correction still identifies them from the thinned tail.
using CairoMakie, PairPlots
prior_chain = sample(corrected_post, Prior(), 1000; progress = false)
pp_keys = [@varname(damp_AR), @varname(std),
@varname(cluster_factor), @varname(init_incidence)]
pairplot(
PairPlots.Series(corrected_chain[pp_keys]; label = "posterior"),
PairPlots.Series(prior_chain[pp_keys]; label = "prior"))
The correction in a figure
generated_observables. Plotting the posterior median and 95% band for all three fits over time — the complete-data reference, the naive truncated fit, and the RightTruncate-corrected fit — shows the right-truncation artefact and its removal in the recent window (shaded).
Rt_complete = rt_bands(complete_post, eventual, complete_chain)
Rt_naive = rt_bands(naive_post, observed_so_far, naive_chain)
Rt_corrected = rt_bands(corrected_post, observed_so_far, corrected_chain)
ts = 1:n
fig = Figure(; size = (760, 420))
ax = Axis(fig[1, 1]; xlabel = "Reference day",
ylabel = "Reproduction number Rₜ")
vspan!(ax, n - 6, n; color = (:grey, 0.15))
rt_line!(ax, ts, Rt_complete; color = :black, label = "complete (reference)")
rt_line!(ax, ts, Rt_naive; color = :crimson, label = "naive (truncated)")
rt_line!(ax, ts, Rt_corrected; color = :seagreen, label = "corrected")
hlines!(ax, [1.0]; color = :grey, linestyle = :dash)
axislegend(ax; position = :lb)
fig
In the shaded recent window the naive fit (red) dips below the complete-data reference (black) — the artefactual late down-turn — while the RightTruncate correction (green) pulls the recent
Reading the shared parameters
Wrapping the error model in RightTruncate does not touch the renewal process, so the corrected fit recovers the same shared parameters: the autoregressive damping damp_AR[1]), the innovation scale std), the observation overdispersion (cluster_factor), and the initial infections (init_incidence). They keep their flat names (prefixing is disabled throughout the package).
sample returns a FlexiChains chain, which summarystats summarises directly — no conversion step — giving point estimates and their uncertainty alongside the effective sample size and
using MCMCChains
summarystats(corrected_chain)╭─FlexiSummary (9 statistics) ─────────────────────────────────────────────────╮
│ iter collapsed │
│ chain collapsed │
│ ↓ stat = [mean, std, mcse, ess_bulk, ess_tail, rhat, q5, q50, q95] │
│ │
│ Parameters (54) ── AbstractPPL.VarName │
│ Float64 ar_init[1], damp_AR[1], std, ϵ_t[1], ϵ_t[2], ϵ_t[3], ϵ_t[4], │
│ ϵ_t[5], ϵ_t[6], ϵ_t[7], ϵ_t[8], ϵ_t[9], ϵ_t[10], ϵ_t[11], ϵ_t[12], │
│ ϵ_t[13], ϵ_t[14], ϵ_t[15], ϵ_t[16], ϵ_t[17], ϵ_t[18], ϵ_t[19], │
│ ϵ_t[20], ϵ_t[21], ϵ_t[22], ϵ_t[23], ϵ_t[24], ϵ_t[25], ϵ_t[26], │
│ ϵ_t[27], ϵ_t[28], ϵ_t[29], ϵ_t[30], ϵ_t[31], ϵ_t[32], ϵ_t[33], │
│ ϵ_t[34], ϵ_t[35], ϵ_t[36], ϵ_t[37], ϵ_t[38], ϵ_t[39], ϵ_t[40], │
│ ϵ_t[41], ϵ_t[42], ϵ_t[43], ϵ_t[44], ϵ_t[45], ϵ_t[46], ϵ_t[47], │
│ ϵ_t[48], ϵ_t[49], init_incidence, cluster_factor │
│ │
│ Extras (14) │
│ Float64 n_steps, is_accept, acceptance_rate, log_density, │
│ hamiltonian_energy, hamiltonian_energy_error, │
│ max_hamiltonian_energy_error, tree_depth, numerical_error, │
│ step_size, nom_step_size, logprior, loglikelihood, logjoint │
│ │
│ Summary │
│ param mean std mcse ess_bulk ess_tail rhat … │
│ ar_init[1] 0.3573 0.2944 0.0162 332.4600 515.8706 1.0042 … │
│ damp_AR[1] 0.8679 0.0462 0.0031 231.7431 462.5285 1.0043 … │
│ std 0.1784 0.0564 0.0083 47.3079 152.0927 1.0128 … │
│ ϵ_t[1] 2.2763 0.8776 0.0354 629.8363 725.8954 1.0021 … │
│ ϵ_t[2] 0.8195 0.8434 0.0302 793.5643 837.1625 0.9991 … │
│ ϵ_t[3] 0.9742 0.8515 0.0228 1406.0714 757.5137 0.9995 … │
│ ϵ_t[4] 0.2449 0.8961 0.0447 400.2101 445.9117 1.0031 … │
│ ϵ_t[5] 0.4095 0.9068 0.0285 1034.7777 720.8138 1.0009 … │
│ ϵ_t[6] 1.7240 0.9294 0.0559 269.3324 648.2117 1.0032 … │
│ ϵ_t[7] 0.1519 0.7899 0.0293 728.4025 576.7201 1.0016 … │
│ ϵ_t[8] 0.2048 0.8077 0.0298 740.5552 692.5039 0.9996 … │
│ ϵ_t[9] 0.9968 0.8370 0.0559 225.5175 622.6519 1.0044 … │
│ ϵ_t[10] -0.7650 0.8125 0.0353 527.6271 703.7770 1.0011 … │
│ ϵ_t[11] 0.2445 0.7661 0.0185 1688.7120 738.9525 0.9993 … │
│ ϵ_t[12] 0.4914 0.7746 0.0210 1277.8583 838.2441 0.9996 … │
│ ϵ_t[13] 0.5745 0.8042 0.0208 1471.9637 635.1793 0.9994 … │
│ ϵ_t[14] 0.3479 0.7479 0.0221 1186.2822 745.3713 1.0075 … │
│ ϵ_t[15] 0.8730 0.7530 0.0209 1325.3849 646.3871 1.0056 … │
│ ϵ_t[16] 0.3834 0.7910 0.0238 1103.3199 651.5570 1.0008 … │
│ ϵ_t[17] 0.0154 0.7796 0.0185 1821.1285 558.9979 1.0088 … │
│ ϵ_t[18] -0.7029 0.8580 0.0459 343.4261 696.1371 0.9995 … │
│ ϵ_t[19] 1.2192 0.8487 0.0449 368.2001 720.5786 1.0011 … │
│ ϵ_t[20] 0.4927 0.7602 0.0191 1545.0323 519.2485 1.0009 … │
│ ϵ_t[21] 0.0227 0.7671 0.0221 1209.2308 759.6198 1.0000 … │
│ ϵ_t[22] 0.3044 0.8302 0.0270 951.6558 498.4579 1.0001 … │
│ ϵ_t[23] -0.2695 0.7837 0.0269 828.7602 684.2067 0.9994 … │
│ ϵ_t[24] 0.3824 0.7506 0.0253 871.3351 717.3066 0.9995 … │
│ ϵ_t[25] -0.1181 0.7373 0.0253 858.0032 642.8424 1.0012 … │
│ ϵ_t[26] 0.3358 0.7253 0.0192 1437.2281 787.0653 1.0001 … │
│ ϵ_t[27] 0.5253 0.7636 0.0218 1225.6534 580.2325 1.0016 … │
│ ϵ_t[28] 0.1756 0.7338 0.0213 1200.4961 763.5188 1.0032 … │
│ ϵ_t[29] -0.0925 0.7808 0.0201 1396.9931 566.2257 1.0014 … │
│ ϵ_t[30] -0.5992 0.7719 0.0221 1213.4942 674.8806 0.9995 … │
│ ϵ_t[31] -0.5576 0.7480 0.0203 1375.5078 696.9227 1.0015 … │
│ ϵ_t[32] -0.0696 0.7644 0.0246 1042.7024 765.0976 0.9992 … │
│ ϵ_t[33] 0.0460 0.7822 0.0228 1191.0075 650.5753 1.0003 … │
│ ϵ_t[34] 0.2158 0.8140 0.0294 735.4550 594.0936 1.0015 … │
│ ϵ_t[35] -0.1422 0.7139 0.0228 968.2825 841.4702 1.0025 … │
│ ϵ_t[36] -0.2210 0.7408 0.0212 1222.0209 635.2568 1.0010 … │
│ ϵ_t[37] -0.5646 0.7225 0.0198 1308.6856 763.2700 1.0028 … │
│ ϵ_t[38] -0.6812 0.6953 0.0192 1320.7332 615.7265 0.9994 … │
│ ϵ_t[39] -0.0934 0.7332 0.0191 1466.0512 749.0589 1.0144 … │
│ ϵ_t[40] 0.3103 0.7472 0.0234 1020.4765 626.6874 1.0028 … │
│ ϵ_t[41] 0.0578 0.7138 0.0222 1042.6943 663.0072 0.9998 … │
│ ϵ_t[42] -0.0128 0.7400 0.0246 909.2509 656.7158 0.9996 … │
│ ϵ_t[43] -0.0233 0.7223 0.0212 1149.7252 719.7487 0.9992 … │
│ ϵ_t[44] -0.2866 0.6997 0.0201 1194.2275 647.1665 0.9991 … │
│ ϵ_t[45] -0.5153 0.7563 0.0237 1001.5204 714.8612 1.0060 … │
│ ϵ_t[46] -0.1828 0.7680 0.0197 1515.6962 663.1914 1.0003 … │
│ ϵ_t[47] 0.6374 0.7385 0.0195 1441.1355 621.6066 1.0050 … │
│ ϵ_t[48] 0.6620 0.8890 0.0221 1663.5732 767.7530 1.0034 … │
│ ϵ_t[49] -0.0968 0.9859 0.0321 932.6397 593.8112 0.9995 … │
│ init_incide… 3.1903 0.2255 0.0097 559.6340 627.0170 1.0001 … │
│ cluster_fac… 0.1745 0.0337 0.0042 68.0632 42.6788 1.0115 … │
╰──────────────────────────────────────────────────────────────────────────────╯From the marginal to the joint
RightTruncate corrects right-truncation by conditioning on each reference day's observed-so-far total — the marginal of the full reference-day × reporting-delay structure. When the delay structure itself is of interest (e.g. reporting that drifts over the outbreak), the ReportTriangle observation model keeps the full 2D reporting triangle and scores it cell by cell; its observed row-sums reconcile with this CDF-scaling to machine precision. The marginal correction here is the cheaper, released-code first step.
References
- S. Abbott, J. Hellewell, R. N. Thompson, K. Sherratt, H. P. Gibbs, N. I. Bosse, J. D. Munday, S. Meakin, E. L. Doughty, J. Y. Chun and others. Estimating the time-varying reproduction number of SARS-CoV-2 using national and subnational case counts. Wellcome Open Research 5, 112 (2020).