Real-time nowcasting: right-truncation vs the reporting triangle

In real time the most recent days of a surveillance series are incomplete. A case with reference day is only counted after a reporting delay, so on day we have seen a fraction of the cases that will eventually be attributed to recent days. Fitting a renewal model to this right-truncated tail without correction produces the familiar artefact of real-time estimation, an apparent late down-turn in that is really the not-yet-reported counts [3].

This package corrects right-truncation two ways, and this tutorial fits and compares both on the same simulated data. RightTruncate is the marginal correction. It keeps only each reference day's observed-so-far total and scales the expected count by the reporting-delay CDF, the EpiNow2-style CDF-scaling nowcast. ReportTriangle is the joint correction. It keeps the full reference-day × reporting-delay triangle and scores it cell by cell, the epinowcast-style nowcast. The two share one reporting-delay kernel, so the triangle's observed row-sums are exactly the observed-so-far totals the marginal correction conditions on. Both leave the renewal infection process and its latent untouched. Both keep the model's as the eventual total, so the nowcast is read back out.

The two corrections

The infection pipeline produces  , the expected eventual total for reference day . Write for the reporting-delay PMF (delay  ) and     for its CDF. A reference day of age    has reported a fraction   of its eventual total.

Summing the joint model's observed cells over recovers the marginal mean    , so the two are consistent by construction. A naive fit drops the delay factor (equivalently  ).

The renewal model

Both corrections wrap the same composed renewal model as the renewal tutorial, an autoregressive process folded into a Renewal infection process observed with a NegativeBinomialError.

using ComposableTuringIDModels, Distributions, Random, Turing, Mooncake
using ADTypes: AutoMooncake
using CSV, DataFrames
Random.seed!(20240625)

adt = AutoMooncake(; config = nothing)

latent = AR(
    damp = [truncated(Normal(0.8, 0.05), 0, 1)],
    init = [Normal(0.0, 0.25)],
    ϵ_t = HierarchicalNormal(std = HalfNormal(0.1)))
renewal = Renewal(; generation_time = Gamma(1.4, 1 / 0.38),
    rt = latent, initialisation = Normal(log(1.0), 1.0))
error = NegativeBinomialError(cluster_factor = HalfNormal(0.1))

Simulate a reporting triangle

We take the fully-reported Italy confirmed-case series as the eventual totals, split each day's total across reporting delays with the delay PMF, then mask the cells not yet reported at  . ReportingPMF discretises the reporting delay with the same CensoredDistributions.jl path the rest of the package uses.

datapath = joinpath(pkgdir(ComposableTuringIDModels),
    "docs", "src", "tutorials", "data", "italy_data.csv")
italy = CSV.read(datapath, DataFrame)
n = 45
eventual = italy.confirm[1:n]                 # eventual (complete) totals

reporting_delay = LogNormal(1.6, 0.5)         # mean ≈ 5.6 days
delay_pmf = ReportingPMF(reporting_delay; D = 10)
pmf = delay_pmf.pmf                            # delays 0 … Dmax
Dmax = length(pmf) - 1

Split each reference day's eventual total across delays with a multinomial draw, then keep only the cells reported by   (  ). The observed-so-far total per reference day is the row-sum of the reported cells, the marginal the CDF-scaling correction sees.

full_triangle = reduce(vcat,
    (permutedims(rand(Multinomial(eventual[t], pmf))) for t in 1:n))
mask = [t + d <= n for t in 1:n, d in 0:Dmax]
reported_triangle = full_triangle .* mask
observed_so_far = vec(sum(reported_triangle, dims = 2))
(complete_tail = eventual[(end - 4):end],
    truncated_tail = observed_so_far[(end - 4):end])

(complete_tail = [4782, 4668, 4585, 4805, 4316], truncated_tail = [2245, 1249, 432, 58, 0])

The recent days are visibly thinned, the most recent day showing only a fraction of its eventual count.

Visualise the data

The reporting triangle is the native data object. Plotting it with AlgebraOfGraphics, the not-yet-reported cells blanked, shows the staircase of missing counts in the recent corner.

using AlgebraOfGraphics, CairoMakie, DataFrames

tri_df = DataFrame(
    reference = repeat(1:n, outer = Dmax + 1),
    delay = repeat(0:Dmax, inner = n),
    count = vec(full_triangle),
    reported = vec(mask))
tri_df.shown = ifelse.(tri_df.reported, Float64.(tri_df.count), NaN)

draw(data(tri_df) *
     mapping(:reference, :delay, :shown => "Reported cases") * visual(Heatmap);
    axis = (xlabel = "Reference day", ylabel = "Reporting delay (days)"))

The same truncation reads as a shortfall in the tail when the observed-so-far row-sums are drawn against the eventual totals.

comp_df = DataFrame(
    reference = repeat(1:n, 2),
    count = vcat(eventual, observed_so_far),
    series = repeat(["eventual total", "observed so far"], inner = n))

draw(data(comp_df) * mapping(:reference, :count, color = :series) *
     visual(Lines); axis = (xlabel = "Reference day", ylabel = "Cases"))

Three fits

We fit the plain renewal model three ways with NUTS, drawing two chains with 1000 warmup iterations each so the cross-chain diagnostic is available, differentiating with Mooncake (see Automatic differentiation backend). The naive fit treats the truncated totals as complete, and motivates the problem. RightTruncate applies the marginal correction to the same observed-so-far totals. ReportTriangle applies the joint correction to the reporting triangle, built through the shared define_y_t hook.

A fourth fit takes the plain model to the complete (untruncated) series, what the analyst would eventually see, as a reference.

The corrections differ from the plain model only in their observation model, so each is derived from it with Accessors.@set rather than rebuilt.

using Accessors

naive_model = IDModel(renewal, error)
rt_obs = RightTruncate(error, ReportingCDF(reporting_delay; D = 10))
rt_model = @set naive_model.observation_model = rt_obs
tri_obs = ReportTriangle(error, delay_pmf)
tri_model = @set naive_model.observation_model = tri_obs

IDModel
├─ infection: Renewal
│  └─ rt: AR
│     └─ ϵ_t: HierarchicalNormal
└─ observation: ReportTriangle
   ├─ error: NegativeBinomialError
   └─ delay: ReportingPMF

A fit is a model and the data it is scored against. The four are a table of pairs, so the sampler call is written once.

specs = (
    naive = (model = naive_model, data = observed_so_far),
    right_truncate = (model = rt_model, data = observed_so_far),
    report_triangle = (model = tri_model,
        data = define_y_t(tri_obs, reported_triangle, eventual)),
    complete = (model = naive_model, data = eventual))

fits = map(specs) do spec
    post = as_turing_model(spec.model, spec.data, n)
    chain = sample(post, NUTS(1000, 0.95; adtype = adt), MCMCThreads(), 250, 2;
        progress = false)
    (; post, chain)
end

naive_post, naive_chain = fits.naive
rt_post, rt_chain = fits.right_truncate
tri_post, tri_chain = fits.report_triangle
complete_post, complete_chain = fits.complete

┌ Info: Found initial step size
└   ϵ = 0.00625
┌ Info: Found initial step size
└   ϵ = 0.0125
┌ Info: Found initial step size
└   ϵ = 0.025
┌ Info: Found initial step size
└   ϵ = 0.0125
┌ Info: Found initial step size
└   ϵ = 0.00625
┌ Info: Found initial step size
└   ϵ = 0.003125
┌ Info: Found initial step size
└   ϵ = 0.0125
┌ Info: Found initial step size
└   ϵ = 0.0125
┌ Warning: There were 1 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4sXd9/src/mcmc/hmc.jl:483
┌ Warning: There were 2 divergent transitions. Consider reparameterising your model or using a smaller step size. For adaptive samplers such as NUTS and HMCDA, consider increasing `target_accept`.
└ @ Turing.Inference ~/.julia/packages/Turing/4sXd9/src/mcmc/hmc.jl:483
┌ Warning: There were 3 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/4sXd9/src/mcmc/hmc.jl:483

Recent Rt

  is a generated quantity, recovered per draw by re-running the fitted model over the chain with generated_observables. Averaging it over the most recent window shows the bias and its removal. Right-truncation biases the naive recent downward, so the naive value sits below the two corrections, which lift it back up towards the complete-data reference. The exact numbers carry Monte Carlo noise on this short run, but the direction repeats.

using Statistics

function recent_Rt(post, chain; window = 7)
    gens = vec(generated_observables(post, nothing, chain).generated)
    Rt_mean = mean(exp.(g.Z_t) for g in gens)
    round(mean(Rt_mean[(end - window + 1):end]), digits = 2)
end

map(fit -> recent_Rt(fit.post, fit.chain), fits)

(naive = 0.3, right_truncate = 0.9, report_triangle = 0.88, complete = 0.91)

Posterior prediction, nowcast, and Rt

Three views of the fits share the reference-day axis with the recent window shaded. The panel plots the bands of all four fits. The nowcast panel plots the reconstructed eventual totals   for the two corrections against the true eventual totals and the truncated observed-so-far. The posterior prediction panel plots the RightTruncate model's posterior-predictive observed-so-far, recovered with predict on the same model with the observations set to missing, against the data it was fit to.

ts = 1:n
Rt_complete = gen_bands(complete_post, complete_chain, g -> exp.(g.Z_t))
Rt_naive = gen_bands(naive_post, naive_chain, g -> exp.(g.Z_t))
Rt_rt = gen_bands(rt_post, rt_chain, g -> exp.(g.Z_t))
Rt_tri = gen_bands(tri_post, tri_chain, g -> exp.(g.Z_t))

now_rt = gen_bands(rt_post, rt_chain, g -> g.I_t)
now_tri = gen_bands(tri_post, tri_chain, g -> g.I_t)

rt_pred = predict(as_turing_model(rt_model, fill(missing, n), n), rt_chain)
yt = predictive_bands(rt_pred, n)

fig = Figure(; size = (780, 820))
ax1 = Axis(fig[1, 1]; ylabel = "Reproduction number Rₜ")
vspan!(ax1, n - 6, n; color = (:grey, 0.15))
ci_ribbon!(ax1, ts, Rt_complete; color = :black, label = "complete (reference)")
ci_ribbon!(ax1, ts, Rt_naive; color = :crimson, label = "naive")
ci_ribbon!(ax1, ts, Rt_rt; color = :seagreen, label = "right-truncate")
ci_ribbon!(ax1, ts, Rt_tri; color = :steelblue, label = "report-triangle")
hlines!(ax1, [1.0]; color = :grey, linestyle = :dash)
axislegend(ax1; position = :lb, nbanks = 2)

ax2 = Axis(fig[2, 1]; ylabel = "Eventual cases (nowcast)")
vspan!(ax2, n - 6, n; color = (:grey, 0.15))
ci_ribbon!(ax2, ts, now_rt; color = :seagreen, label = "right-truncate")
ci_ribbon!(ax2, ts, now_tri; color = :steelblue, label = "report-triangle")
lines!(ax2, ts, Float64.(eventual); color = :black, linewidth = 2,
    label = "true eventual")
scatter!(ax2, ts, observed_so_far; color = :grey, markersize = 6,
    label = "observed so far")
axislegend(ax2; position = :lt)

ax3 = Axis(fig[3, 1]; xlabel = "Reference day", ylabel = "Observed-so-far")
vspan!(ax3, n - 6, n; color = (:grey, 0.15))
ci_ribbon!(ax3, ts, yt; color = :teal, label = "posterior predictive")
scatter!(ax3, ts, observed_so_far; color = :black, markersize = 6,
    label = "observed")
axislegend(ax3; position = :lt)
fig

In the shaded recent window the naive (crimson) dips below the complete-data reference (black), while the two corrections lift the recent back up off that spurious decline. The nowcast panel shows the same story on the count scale, the corrections lifting the reconstructed eventual totals above the truncated observed-so-far towards the true eventual line. The posterior-predictive panel confirms the RightTruncate fit reproduces the observed-so-far series it saw.

Shared parameters

Neither correction touches the renewal process, so both recover the same shared parameters, the autoregressive damping (damp[1]), the innovation scale (std), the observation overdispersion (cluster_factor) and the initial infections (init_incidence).

using MCMCChains
summarystats(tri_chain)

╭─FlexiSummary (9 statistics) ─────────────────────────────────────────────────╮
│   iter    collapsed                                                          │
│   chain   collapsed                                                          │
│ ↓ stat  = [mean, std, mcse, ess_bulk, ess_tail, rhat, q5, q50, q95]          │
│                                                                              │
│ Parameters (50) ── AbstractPPL.VarName                                       │
│  Float64  init[1], damp, ρ, 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], 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  …       │
│        init[1]   0.1086  0.2931  0.0400   53.9299  110.1585  1.0092  …       │
│           damp   0.7846  0.0452  0.0136   11.9950   24.3737  1.1169  …       │
│              ρ   0.7846  0.0452  0.0136   11.9950   24.3737  1.1169  …       │
│            std   0.3144  0.0293  0.0080   13.5512   51.3326  1.0402  …       │
│         ϵ_t[1]   3.0080  0.7570  0.1182   40.5647   75.7999  0.9998  …       │
│         ϵ_t[2]  -0.0795  0.6563  0.1023   39.6271   98.0348  1.0412  …       │
│         ϵ_t[3]   1.2543  0.4957  0.0583   71.7512  121.9873  1.0076  …       │
│         ϵ_t[4]  -0.3354  0.4293  0.0661   42.4113   92.7750  1.0198  …       │
│         ϵ_t[5]  -0.4061  0.4375  0.0531   66.1227  179.1219  1.0068  …       │
│         ϵ_t[6]   2.9800  0.4335  0.0714   37.4075  154.6013  1.0334  …       │
│         ϵ_t[7]  -0.5158  0.3969  0.0527   53.5915   73.9288  1.0119  …       │
│         ϵ_t[8]  -0.2303  0.3081  0.0333   85.6160   81.0039  1.0241  …       │
│         ϵ_t[9]   2.3312  0.3164  0.0911   12.4592  168.7149  1.0854  …       │
│        ϵ_t[10]  -1.8666  0.3157  0.0555   33.8225   49.1867  1.0071  …       │
│        ϵ_t[11]   0.6802  0.2158  0.0285   59.0324  139.4831  1.0463  …       │
│        ϵ_t[12]   0.5810  0.2226  0.0303   55.9988  117.9778  1.0044  …       │
│        ϵ_t[13]   0.7059  0.2330  0.0370   45.5873   40.6738  1.0369  …       │
│        ϵ_t[14]  -0.1081  0.1874  0.0291   43.4217   90.8820  1.0218  …       │
│        ϵ_t[15]   1.3177  0.1887  0.0386   24.7272   93.7628  1.0656  …       │
│        ϵ_t[16]   0.3412  0.1891  0.0444   17.8535   59.8187  1.0962  …       │
│        ϵ_t[17]   0.3383  0.1630  0.0235   49.9377   92.1203  1.0421  …       │
│        ϵ_t[18]  -2.1157  0.2543  0.0638   16.1399   57.0774  1.0260  …       │
│        ϵ_t[19]   2.5832  0.2597  0.0648   16.2388  117.1820  1.0219  …       │
│        ϵ_t[20]   0.2857  0.1482  0.0313   22.4290   36.4198  1.0774  …       │
│        ϵ_t[21]  -0.3213  0.1449  0.0318   22.8188   36.0683  1.0328  …       │
│        ϵ_t[22]   0.8797  0.1225  0.0500    6.1554  108.9433  1.1397  …       │
│        ϵ_t[23]  -0.8323  0.1471  0.0322   21.3275   50.6233  1.0371  …       │
│        ϵ_t[24]   1.0097  0.1270  0.0535    5.6866   50.5214  1.1513  …       │
│        ϵ_t[25]  -0.5082  0.1214  0.0261   23.5860   29.8464  1.0087  …       │
│        ϵ_t[26]   0.4902  0.0865  0.0153   32.6622  143.3287  1.0360  …       │
│        ϵ_t[27]   0.6835  0.1032  0.0341    9.7906   74.7083  1.1108  …       │
│        ϵ_t[28]   0.2746  0.0906  0.0199   21.0484  143.9406  1.0755  …       │
│        ϵ_t[29]   0.1741  0.0953  0.0187   26.2251   58.4054  1.0562  …       │
│        ϵ_t[30]  -0.6275  0.1059  0.0222   24.7746   50.2483  1.0118  …       │
│        ϵ_t[31]  -0.5618  0.0915  0.0185   23.9872  137.6375  1.0026  …       │
│        ϵ_t[32]   0.2610  0.0703  0.0107   43.4444  100.4579  1.0039  …       │
│        ϵ_t[33]  -0.0369  0.0681  0.0077   79.8381  124.1532  1.0105  …       │
│        ϵ_t[34]   0.5212  0.0783  0.0140   31.6137   49.5578  1.0135  …       │
│        ϵ_t[35]  -0.1254  0.0733  0.0141   32.4372   47.5360  1.0347  …       │
│        ϵ_t[36]   0.0078  0.0682  0.0078   79.3643   97.0484  1.0183  …       │
│        ϵ_t[37]  -0.4463  0.0831  0.0167   26.0215   53.1560  1.0043  …       │
│        ϵ_t[38]  -0.8763  0.1077  0.0178   36.7674  110.0285  1.0195  …       │
│        ϵ_t[39]  -0.0174  0.1015  0.0256   16.9070   88.2224  1.1023  …       │
│        ϵ_t[40]   0.5450  0.1246  0.0277   20.8237  100.1520  1.0098  …       │
│        ϵ_t[41]  -0.0911  0.1274  0.0164   63.4532   70.1062  1.0254  …       │
│        ϵ_t[42]  -0.1293  0.1776  0.0225   63.6180  124.5383  1.0020  …       │
│        ϵ_t[43]  -0.0738  0.4111  0.0452   83.1611  205.1149  1.0009  …       │
│        ϵ_t[44]  -0.1845  1.0468  0.1056  106.7008   79.2638  1.0018  …       │
│   init_incide…   2.9568  0.2548  0.0405   40.3687   75.0883  1.0075  …       │
│   cluster_fac…   0.0102  0.0061  0.0005  125.1904   79.0610  0.9996  …       │
╰──────────────────────────────────────────────────────────────────────────────╯

Prior versus posterior

Sampling the reporting-triangle model with Prior and overlaying it on the posterior with PairPlots.jl confirms the joint correction still identifies the shared parameters from the thinned triangle.

using PairPlots

prior_chain = sample(tri_post, Prior(), 1000; progress = false)
pp_keys = [@varname(damp), @varname(std),
    @varname(cluster_factor), @varname(init_incidence)]
pairplot(
    PairPlots.Series(tri_chain[pp_keys]; label = "posterior"),
    PairPlots.Series(prior_chain[pp_keys]; label = "prior"))

References

2. 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).