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, and both keep the model's as the eventual total, so the nowcast is just 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.

naive_model = IDModel(renewal, error)
naive_post = as_turing_model(naive_model, observed_so_far, n)
naive_chain = sample(
    naive_post, NUTS(1000, 0.95; adtype = adt), MCMCThreads(), 250, 2;
    progress = false)

rt_obs = RightTruncate(error, ReportingCDF(reporting_delay; D = 10))
rt_model = IDModel(renewal, rt_obs)
rt_post = as_turing_model(rt_model, observed_so_far, n)
rt_chain = sample(
    rt_post, NUTS(1000, 0.95; adtype = adt), MCMCThreads(), 250, 2;
    progress = false)

tri_obs = ReportTriangle(error, delay_pmf)
tri_data = define_y_t(tri_obs, reported_triangle, eventual)
tri_model = IDModel(renewal, tri_obs)
tri_post = as_turing_model(tri_model, tri_data, n)
tri_chain = sample(
    tri_post, NUTS(1000, 0.95; adtype = adt), MCMCThreads(), 250, 2;
    progress = false)

┌ Info: Found initial step size
└   ϵ = 0.00625
┌ Info: Found initial step size
└   ϵ = 0.0125
┌ Warning: There were 49 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 49 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
┌ Info: Found initial step size
└   ϵ = 0.025
┌ 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/WteH7/src/mcmc/hmc.jl:483
┌ 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/WteH7/src/mcmc/hmc.jl:483
┌ Info: Found initial step size
└   ϵ = 0.00625
┌ Info: Found initial step size
└   ϵ = 0.003125

As a reference we also fit the plain model to the complete (untruncated) series, what the analyst would eventually see.

complete_post = as_turing_model(naive_model, eventual, n)
complete_chain = sample(
    complete_post, NUTS(1000, 0.95; adtype = adt), MCMCThreads(), 250, 2;
    progress = false)

┌ Info: Found initial step size
└   ϵ = 0.0125
┌ Info: Found initial step size
└   ϵ = 0.0125

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 is the robust, repeatable signal.

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

(complete = recent_Rt(complete_post, complete_chain),
    naive = recent_Rt(naive_post, naive_chain),
    right_truncate = recent_Rt(rt_post, rt_chain),
    report_triangle = recent_Rt(tri_post, tri_chain))

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

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_AR[1]), the innovation scale (std), the observation overdispersion (cluster_factor) and the initial infections (init_incidence). sample returns a FlexiChains chain that summarystats summarises directly, giving point estimates and their uncertainty alongside the effective sample size and .

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  ar_init[1], damp_AR, ρ, 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  …       │
│     ar_init[1]   0.1084  0.2932  0.0399   54.1156  121.4212  1.0078  …       │
│        damp_AR   0.7846  0.0447  0.0120   14.6666   24.3737  1.1018  …       │
│              ρ   0.7846  0.0447  0.0120   14.6666   24.3737  1.1018  …       │
│            std   0.3152  0.0296  0.0080   13.9895   62.4321  1.0315  …       │
│         ϵ_t[1]   3.0074  0.7609  0.1202   39.6618   73.5218  1.0005  …       │
│         ϵ_t[2]  -0.0837  0.6554  0.0969   45.1590   91.3262  1.0382  …       │
│         ϵ_t[3]   1.2524  0.4957  0.0578   72.8974  129.9935  1.0060  …       │
│         ϵ_t[4]  -0.3367  0.4296  0.0649   43.9172   80.8746  1.0152  …       │
│         ϵ_t[5]  -0.4068  0.4345  0.0527   66.0618  179.1219  1.0078  …       │
│         ϵ_t[6]   2.9749  0.4331  0.0704   38.1948  146.3276  1.0277  …       │
│         ϵ_t[7]  -0.5139  0.3971  0.0519   54.5130   84.7886  1.0094  …       │
│         ϵ_t[8]  -0.2319  0.3064  0.0331   86.1893   94.5736  1.0233  …       │
│         ϵ_t[9]   2.3243  0.3170  0.0758   17.6728  146.6883  1.0682  …       │
│        ϵ_t[10]  -1.8610  0.3178  0.0546   35.1686   54.7904  1.0034  …       │
│        ϵ_t[11]   0.6792  0.2130  0.0263   67.5400  142.7571  1.0428  …       │
│        ϵ_t[12]   0.5795  0.2241  0.0303   56.4840  117.9778  1.0023  …       │
│        ϵ_t[13]   0.7046  0.2298  0.0355   48.5727   45.8115  1.0335  …       │
│        ϵ_t[14]  -0.1079  0.1876  0.0286   45.3207  114.7050  1.0179  …       │
│        ϵ_t[15]   1.3145  0.1867  0.0358   27.9670   93.7628  1.0558  …       │
│        ϵ_t[16]   0.3402  0.1873  0.0414   20.3260   63.9414  1.0849  …       │
│        ϵ_t[17]   0.3385  0.1608  0.0216   56.4605   93.8624  1.0341  …       │
│        ϵ_t[18]  -2.1120  0.2542  0.0640   16.0316   67.0725  1.0222  …       │
│        ϵ_t[19]   2.5780  0.2596  0.0644   16.5652  109.5566  1.0186  …       │
│        ϵ_t[20]   0.2843  0.1461  0.0285   26.4102   57.9988  1.0622  …       │
│        ϵ_t[21]  -0.3206  0.1437  0.0304   24.7140   36.0683  1.0288  …       │
│        ϵ_t[22]   0.8771  0.1216  0.0456    7.2978  133.7093  1.1188  …       │
│        ϵ_t[23]  -0.8285  0.1466  0.0318   21.8355   50.6233  1.0292  …       │
│        ϵ_t[24]   1.0053  0.1258  0.0496    6.4828   50.5214  1.1335  …       │
│        ϵ_t[25]  -0.5072  0.1208  0.0253   24.5958   32.9753  1.0067  …       │
│        ϵ_t[26]   0.4894  0.0862  0.0150   34.1629  150.9689  1.0292  …       │
│        ϵ_t[27]   0.6818  0.1022  0.0308   11.6543   74.7083  1.0948  …       │
│        ϵ_t[28]   0.2740  0.0891  0.0183   24.8455  131.6505  1.0631  …       │
│        ϵ_t[29]   0.1734  0.0946  0.0175   29.0117   64.9601  1.0446  …       │
│        ϵ_t[30]  -0.6264  0.1054  0.0218   25.2698   53.4113  1.0111  …       │
│        ϵ_t[31]  -0.5603  0.0917  0.0186   23.6200  137.8237  1.0010  …       │
│        ϵ_t[32]   0.2603  0.0705  0.0108   43.0246  101.8865  1.0059  …       │
│        ϵ_t[33]  -0.0364  0.0682  0.0079   75.8553  104.6227  1.0111  …       │
│        ϵ_t[34]   0.5200  0.0787  0.0141   31.5827   49.9605  1.0154  …       │
│        ϵ_t[35]  -0.1253  0.0732  0.0139   33.2577   48.9546  1.0329  …       │
│        ϵ_t[36]   0.0083  0.0679  0.0078   79.3147  100.6130  1.0128  …       │
│        ϵ_t[37]  -0.4454  0.0825  0.0170   24.7522   45.9898  1.0076  …       │
│        ϵ_t[38]  -0.8741  0.1081  0.0175   38.1229  108.3424  1.0151  …       │
│        ϵ_t[39]  -0.0173  0.1012  0.0216   26.4797   94.1350  1.0894  …       │
│        ϵ_t[40]   0.5435  0.1241  0.0278   20.2096  107.3457  1.0108  …       │
│        ϵ_t[41]  -0.0901  0.1267  0.0162   64.1568   56.2637  1.0246  …       │
│        ϵ_t[42]  -0.1299  0.1780  0.0229   61.3783  124.5383  1.0016  …       │
│        ϵ_t[43]  -0.0731  0.4112  0.0457   81.1615  226.8501  1.0007  …       │
│        ϵ_t[44]  -0.1853  1.0447  0.1059  106.7740   79.2638  1.0008  …       │
│   init_incide…   2.9553  0.2546  0.0404   40.3484   75.0883  1.0058  …       │
│   cluster_fac…   0.0101  0.0061  0.0005  115.8496   99.9814  1.0005  …       │
╰──────────────────────────────────────────────────────────────────────────────╯

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_AR), @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).