Multiple observation streams: cases, deaths, and strata
Real-time surveillance rarely watches an epidemic through a single lens. The same infections surface as reported cases, hospital admissions, deaths, and often each of these split by age, region, or variant. These streams share one underlying infection process but differ in their reporting delay, ascertainment, and noise [4]. Fitting them jointly — one infection trajectory, several observation streams — propagates uncertainty correctly and lets a sparse stream (deaths) borrow strength from a dense one (cases).
This tutorial uses one construct, Split, for every multi-stream shape. Split observes the expected series arriving at the point where it sits in the pipeline through several named streams, so where you place it chooses the composition:
parallel — placed high, on infections: every stream observes the same
(cases and deaths each a delayed, ascertained fraction of ); cascade — placed low, after a shared layer: a later stream is observed downstream of an earlier one (deaths as a delayed fraction of the expected reported cases);
strata — one stream per data-defined group (an age band).
How Split threads streams
Every observation model in the package returns the uniform pair (; y_t, expected): the sampled observations y_t and the pre-error expected series the error was scored against. Exposing expected is what lets Split do all three shapes with one mechanism. Split feeds each stream the expected series reaching it, and — because Split is itself an observation model — a shared modifier can run before it. Split((cases = …, deaths = …)) on its own splits infections (parallel), while LatentDelay(Split((cases = …, deaths = …)), pmf) applies a common delay first and then splits, so a stream nested inside another stream's pipeline sits downstream of it (cascade).
The threaded quantity is the expected, not the realised, series
A downstream stream reads its upstream stream's expected (pre-error) series, never its realised, sampled counts. So a cascade threads the mean reported cases into deaths, not a noisy draw. The case where an observation depends on another stream's realised (error-corrupted) observation — feeding sampled cases, not expected cases, into deaths — is not covered here and is out of scope for now.
Split also prefixes each stream's sampled variables with the stream name automatically, so the streams stay distinct without any manual prefix layer.
Parallel: cases and deaths from shared infections
We drive the streams with a renewal infection process, exactly as in the renewal tutorial, and observe it through two pipelines. Cases are a short-delay, high-ascertainment negative-binomial stream. Deaths are a long-delay stream whose ascertainment — the infection-fatality ratio — is itself estimated: each stream is a full observation model, so its ascertainment can be a fixed fraction or, as here, a latent Intercept model with a prior.
using ComposableTuringIDModels, Distributions, Random, Turing, Mooncake
using ADTypes: AutoMooncake
Random.seed!(1234)
latent = AR(
damp = [truncated(Normal(0.8, 0.05), 0, 1),
truncated(Normal(0.1, 0.05), 0, 1)],
init = [Normal(0.0, 0.2), Normal(0.0, 0.2)],
ϵ_t = HierarchicalNormal(std = HalfNormal(0.1)))
renewal = Renewal(; generation_time = Gamma(6.5, 0.62),
rt = latent, initialisation = Normal(log(100.0), 0.1))
cases = LatentDelay(
Ascertainment(NegativeBinomialError(cluster_factor = HalfNormal(0.1)),
FixedIntercept(log(0.6))), # ~60% case ascertainment
LogNormal(1.6, 0.5)) # short infection→report delay
deaths = LatentDelay(
Ascertainment(NegativeBinomialError(cluster_factor = HalfNormal(0.1)),
Intercept(Normal(log(0.015), 0.25))), # estimated ~1.5% IFR
LogNormal(2.8, 0.4)) # long infection→death delay
parallel = Split((cases = cases, deaths = deaths))Split
├─ cases: LatentDelay
│ └─ model: Ascertainment
│ ├─ model: NegativeBinomialError
│ └─ latent: PrefixLatentModel
│ └─ model: FixedIntercept
└─ deaths: LatentDelay
└─ model: Ascertainment
├─ model: NegativeBinomialError
└─ latent: PrefixLatentModel
└─ model: InterceptThe composed model assembles the renewal infection process and the two-stream observation model exactly like a single-stream study.
model = IDModel(renewal, parallel)IDModel
├─ infection: Renewal
│ └─ rt: AR
│ └─ ϵ_t: HierarchicalNormal
└─ observation: Split
├─ cases: LatentDelay
│ └─ model: Ascertainment
│ ├─ model: NegativeBinomialError
│ └─ latent: PrefixLatentModel
│ └─ model: FixedIntercept
└─ deaths: LatentDelay
└─ model: Ascertainment
├─ model: NegativeBinomialError
└─ latent: PrefixLatentModel
└─ model: InterceptPassing missing data simulates a synthetic outbreak. The per-stream data contract is a NamedTuple keyed by stream name, and the returned generated_y_t is a NamedTuple of the two simulated series.
n = 70
sim = as_turing_model(model, (cases = missing, deaths = missing), n)()
y = sim.generated_y_t
(total_cases = sum(skipmissing(y.cases)), total_deaths = sum(skipmissing(y.deaths)))(total_cases = 6576, total_deaths = 82)Fitting conditions on both streams at once. We draw a full chain with NUTS, matching the other tutorials, and differentiate with Mooncake, the recommended backend for this package (see Automatic differentiation backend).
ydata = (cases = y.cases, deaths = y.deaths)
posterior = as_turing_model(model, ydata, n)
chain = sample(
posterior, NUTS(0.95; adtype = AutoMooncake(; config = nothing)),
MCMCThreads(), 250, 2; progress = false)┌ Info: Found initial step size
└ ϵ = 0.025
┌ Info: Found initial step size
└ ϵ = 0.00625
┌ 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/WteH7/src/mcmc/hmc.jl:483
┌ Warning: There were 7 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 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:483The two streams keep their own overdispersion parameters — Split prefixes them cases.cluster_factor and deaths.cluster_factor — while sharing the one infection trajectory, and the deaths stream's estimated IFR intercept (deaths.Ascertainment.intercept) is recovered alongside them. The dense case stream pins the shared
using MCMCChains
summarystats(chain)╭─FlexiSummary (9 statistics) ─────────────────────────────────────────────────╮
│ iter collapsed │
│ chain collapsed │
│ ↓ stat = [mean, std, mcse, ess_bulk, ess_tail, rhat, q5, q50, q95] │
│ │
│ Parameters (77) ── AbstractPPL.VarName │
│ Float64 ar_init[1], ar_init[2], damp_AR[1], damp_AR[2], 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], ϵ_t[50], ϵ_t[51], │
│ ϵ_t[52], ϵ_t[53], ϵ_t[54], ϵ_t[55], ϵ_t[56], ϵ_t[57], ϵ_t[58], │
│ ϵ_t[59], ϵ_t[60], ϵ_t[61], ϵ_t[62], ϵ_t[63], ϵ_t[64], ϵ_t[65], │
│ ϵ_t[66], ϵ_t[67], ϵ_t[68], init_incidence, cases.cluster_factor, │
│ deaths.Ascertainment.intercept, deaths.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.0208 0.2048 0.0115 315.6986 292.0383 1.0092 … │
│ ar_init[2] 0.2999 0.1489 0.0055 732.7980 425.5089 1.0024 … │
│ damp_AR[1] 0.7888 0.0464 0.0021 467.3546 457.1127 1.0045 … │
│ damp_AR[2] 0.0913 0.0420 0.0020 381.9998 207.7303 0.9991 … │
│ std 0.0705 0.0326 0.0016 399.2614 341.8850 1.0003 … │
│ ϵ_t[1] 0.4029 0.9405 0.0386 589.8113 293.6967 1.0008 … │
│ ϵ_t[2] 0.3833 1.0264 0.0361 837.0544 347.5764 0.9983 … │
│ ϵ_t[3] 0.3096 1.0927 0.0357 925.8997 313.0064 1.0170 … │
│ ϵ_t[4] 0.3156 0.9653 0.0360 707.6225 360.1761 0.9994 … │
│ ϵ_t[5] 0.2831 1.1361 0.0467 597.3467 237.7621 1.0054 … │
│ ϵ_t[6] 0.2684 0.9213 0.0409 525.9125 197.5997 0.9983 … │
│ ϵ_t[7] 0.1342 1.0461 0.0369 832.3686 328.4725 0.9985 … │
│ ϵ_t[8] 0.1365 0.9614 0.0355 733.1619 165.6432 1.0114 … │
│ ϵ_t[9] 0.0629 1.0326 0.0407 648.7640 229.0981 0.9995 … │
│ ϵ_t[10] -0.0494 0.9850 0.0390 628.6982 360.2849 1.0000 … │
│ ϵ_t[11] -0.0508 0.9987 0.0371 743.2081 323.9116 0.9981 … │
│ ϵ_t[12] -0.1491 1.0132 0.0397 648.3733 236.1593 0.9980 … │
│ ϵ_t[13] -0.1924 0.8694 0.0351 617.7515 353.1196 0.9989 … │
│ ϵ_t[14] -0.1502 0.8829 0.0302 877.1226 205.8649 1.0063 … │
│ ϵ_t[15] -0.1968 0.9086 0.0395 524.0056 372.7810 0.9994 … │
│ ϵ_t[16] -0.1427 0.9559 0.0321 877.2210 333.2991 1.0023 … │
│ ϵ_t[17] -0.1577 0.9935 0.0375 687.2348 293.4143 0.9999 … │
│ ϵ_t[18] -0.1469 0.9962 0.0350 832.0828 373.1149 1.0006 … │
│ ϵ_t[19] -0.1378 0.9429 0.0398 541.0798 295.8846 0.9981 … │
│ ϵ_t[20] -0.1427 0.9437 0.0353 741.9526 332.4821 0.9980 … │
│ ϵ_t[21] -0.1878 0.9506 0.0358 702.1864 283.0639 0.9990 … │
│ ϵ_t[22] -0.2330 0.9004 0.0389 528.7714 350.3450 0.9997 … │
│ ϵ_t[23] -0.3164 1.0132 0.0397 677.3278 400.7875 1.0079 … │
│ ϵ_t[24] -0.3501 0.9529 0.0320 888.4442 315.3736 1.0020 … │
│ ϵ_t[25] -0.2400 0.9349 0.0361 668.1323 181.7616 1.0213 … │
│ ϵ_t[26] -0.2112 0.8979 0.0329 760.5667 438.4716 1.0130 … │
│ ϵ_t[27] -0.1772 1.0127 0.0430 554.9067 375.4466 1.0070 … │
│ ϵ_t[28] -0.2809 0.9967 0.0376 716.6970 321.5575 1.0044 … │
│ ϵ_t[29] -0.3338 0.9813 0.0412 548.8337 330.8633 0.9992 … │
│ ϵ_t[30] -0.3603 0.9635 0.0341 798.1880 338.1784 1.0067 … │
│ ϵ_t[31] -0.3027 0.9838 0.0352 776.6440 355.1647 0.9995 … │
│ ϵ_t[32] -0.3087 0.9216 0.0354 675.9847 339.7730 1.0134 … │
│ ϵ_t[33] -0.4439 0.9234 0.0440 431.0583 302.3305 1.0069 … │
│ ϵ_t[34] -0.4351 1.0125 0.0320 1006.4291 178.0587 1.0029 … │
│ ϵ_t[35] -0.4590 1.0596 0.0446 570.3955 307.2606 1.0087 … │
│ ϵ_t[36] -0.3986 0.9481 0.0412 526.1370 237.3993 1.0014 … │
│ ϵ_t[37] -0.2621 0.9743 0.0393 620.3654 293.2963 1.0029 … │
│ ϵ_t[38] -0.0314 1.0121 0.0370 738.3547 294.0855 1.0055 … │
│ ϵ_t[39] -0.0233 0.9967 0.0371 716.6462 372.2645 0.9982 … │
│ ϵ_t[40] 0.0984 0.9230 0.0342 788.8379 306.5749 1.0001 … │
│ ϵ_t[41] 0.1884 1.0356 0.0407 630.9454 372.2645 1.0009 … │
│ ϵ_t[42] 0.2958 0.9549 0.0385 581.1140 237.7131 1.0010 … │
│ ϵ_t[43] 0.4177 0.8771 0.0363 584.6554 283.3097 0.9984 … │
│ ϵ_t[44] 0.4980 1.0005 0.0355 764.9530 375.4466 1.0037 … │
│ ϵ_t[45] 0.3737 1.0298 0.0459 496.0364 342.9173 0.9980 … │
│ ϵ_t[46] 0.2857 0.9767 0.0394 622.9815 311.3205 1.0005 … │
│ ϵ_t[47] 0.2300 0.9717 0.0380 672.4036 209.7994 0.9981 … │
│ ϵ_t[48] 0.2903 1.0349 0.0380 758.5418 259.6868 1.0089 … │
│ ϵ_t[49] 0.3694 0.9295 0.0365 642.7308 299.0191 0.9984 … │
│ ϵ_t[50] 0.2973 0.9543 0.0322 905.3370 343.0167 0.9985 … │
│ ϵ_t[51] 0.3912 0.8956 0.0362 622.5015 337.8505 1.0179 … │
│ ϵ_t[52] 0.2793 0.9766 0.0356 759.5582 450.5307 0.9994 … │
│ ϵ_t[53] 0.1183 0.9380 0.0403 536.6579 317.1508 1.0082 … │
│ ϵ_t[54] 0.0477 0.9927 0.0363 736.2650 264.5985 0.9986 … │
│ ϵ_t[55] -0.1000 0.9469 0.0347 710.6326 327.4149 1.0159 … │
│ ϵ_t[56] -0.0975 1.0011 0.0420 580.8342 349.5218 1.0112 … │
│ ϵ_t[57] -0.0216 0.9942 0.0392 669.2090 366.8390 1.0009 … │
│ ϵ_t[58] -0.0176 1.0233 0.0338 937.5874 473.4051 1.0033 … │
│ ϵ_t[59] 0.0986 0.8959 0.0423 465.0118 282.4825 1.0051 … │
│ ϵ_t[60] 0.2520 1.0046 0.0425 553.9923 240.0420 1.0232 … │
│ ϵ_t[61] 0.1623 0.9890 0.0399 636.0954 226.5309 1.0123 … │
│ ϵ_t[62] 0.1916 0.9878 0.0380 744.9858 382.0824 1.0021 … │
│ ϵ_t[63] 0.1007 1.0493 0.0418 608.1062 277.7236 1.0083 … │
│ ϵ_t[64] 0.0609 0.9113 0.0368 625.5258 353.1196 1.0025 … │
│ ϵ_t[65] 0.0273 0.9526 0.0330 819.0542 279.7720 0.9982 … │
│ ϵ_t[66] 0.0404 1.0565 0.0366 826.3426 397.1099 1.0019 … │
│ ϵ_t[67] 0.0493 0.9167 0.0389 562.5097 337.3775 1.0077 … │
│ ϵ_t[68] -0.0023 0.9467 0.0371 664.1802 462.0680 1.0065 … │
│ init_incide… 4.6547 0.0913 0.0040 530.0366 470.2391 1.0009 … │
│ cases.clust… 0.0942 0.0197 0.0009 438.5740 309.9413 0.9991 … │
│ deaths.Asce… -4.1884 0.1052 0.0052 416.1446 240.6091 0.9981 … │
│ deaths.clus… 0.1089 0.0754 0.0026 558.3697 228.1906 0.9997 … │
╰──────────────────────────────────────────────────────────────────────────────╯Cascade: deaths downstream of reported cases
In the parallel model, cases and deaths both branch off infections, so a reporting artefact in the case series (a weekend dip, an ascertainment change) does not touch deaths. Sometimes we want the opposite: deaths modelled as a delayed fraction of the reported cases, so whatever is reflected in cases propagates into deaths. That is a cascade Split placed lower in the stack. Share the infection→case-report delay, then split: the cases stream applies its error to the delayed expectation, and the deaths stream sits downstream, delayed again by the case-report→death interval and scaled by the fatality fraction.
cascade = LatentDelay( # infection→case delay
Split((
cases = NegativeBinomialError(cluster_factor = HalfNormal(0.1)),
deaths = LatentDelay( # case→death delay
Ascertainment(NegativeBinomialError(cluster_factor = HalfNormal(0.1)),
FixedIntercept(log(0.02))),
LogNormal(2.2, 0.3)))),
LogNormal(1.6, 0.5))
cascade_model = IDModel(renewal, cascade)
cas = as_turing_model(cascade_model, (cases = missing, deaths = missing), n)()(generated_y_t = (cases = Union{Missing, Int64}[270, 218, 297, 216, 337, 350, 349, 317, 345, 235 … 18, 24, 42, 36, 27, 32, 26, 28, 15, 16], deaths = Union{Missing, Int64}[missing, missing, missing, missing, missing, missing, missing, missing, missing, missing … 2, 0, 1, 1, 2, 1, 1, 0, 0, 0]), expected_y_t = (cases = [227.09297655006515, 243.34766824296463, 261.4816875378536, 278.27562743143085, 290.96278567072767, 301.28182061015445, 312.3983423411687, 324.84589430964934, 335.202679731171, 339.20058769222027 … 31.186316265782253, 28.694989824511705, 26.443695240615767, 24.485114690631935, 22.78252119789874, 21.169716185551703, 19.462217764470928, 17.66258532039228, 15.938218082601624, 14.34144459929579], deaths = [6.438898577278827, 6.5269959308405525, 6.567662962367721, 6.564000144659306, 6.517047629226356, 6.424169592255847, 6.280901953934637, 6.084810397911528, 5.83805796271519, 5.547591893051596 … 1.1520649516076364, 1.0772106846356608, 1.0106802434475666, 0.9503586772567195, 0.8943779493996301, 0.8411996384132008, 0.7896930783313534, 0.7392927313360437, 0.6900632656439989, 0.6424963952669404]), I_t = [124.82023240024644, 118.78877067901828, 125.14877354497487, 134.7185804215374, 158.58527604659736, 156.77778466087594, 176.86934435498142, 184.48537949475954, 201.74499799634336, 237.88281397858523 … 20.165685555311363, 18.84224824466163, 17.72302898875249, 15.751613393207492, 12.275154235021288, 11.731027928240005, 10.62287226771338, 9.386131507924807, 7.465721130334034, 5.331369642723158], Z_t = [0.5312446389953557, 0.3423504104670212, 0.26194720091244966, 0.23095905875602335, 0.32222357859830564, 0.25145249562109956, 0.3032385726661695, 0.2704627134438929, 0.28597582344222056, 0.3785932527760605 … -0.30701422716001253, -0.28013584652799267, -0.2573887741920619, -0.3042302765638896, -0.4888920355466817, -0.4624106916544015, -0.4674184012263106, -0.47696013976741813, -0.5856294205914864, -0.8032070473681383])The Split sits after the shared case delay and before the error leaves, so the deaths stream's expected input is the delayed-and-ascertained expected cases, not the raw infections: it is both scaled by the fatality fraction and shortened by the case delay.
(cases_expected_length = length(cas.expected_y_t.cases),
deaths_expected_length = length(cas.expected_y_t.deaths),
deaths_are_a_fraction_of_cases =
sum(cas.expected_y_t.deaths) < sum(cas.expected_y_t.cases))(cases_expected_length = 55, deaths_expected_length = 38, deaths_are_a_fraction_of_cases = true)Strata: one stream per age band
A stratified stream — one observation series per age band, region, or variant — is again the same construct, here composed with the renewal infection process and observed through one named stream per band. Each band is a full observation model, so its delay and ascertainment can differ, and its parameters are namespaced by the band name.
strata_obs = Split((
young = LatentDelay(
Ascertainment(NegativeBinomialError(cluster_factor = HalfNormal(0.1)),
FixedIntercept(log(0.7))), LogNormal(1.5, 0.4)),
old = LatentDelay(
Ascertainment(NegativeBinomialError(cluster_factor = HalfNormal(0.1)),
FixedIntercept(log(0.4))), LogNormal(1.8, 0.4))))
strata_model = IDModel(renewal, strata_obs)
strata_sim = as_turing_model(
strata_model, (young = missing, old = missing), n)().generated_y_t
map(s -> sum(skipmissing(s)), strata_sim) # totals per band(young = 2733, old = 1504)The streams above each observe the same infections. When the streams instead draw on a weighted mix of infections — one band, another band, and a summed total — the same Split carries an observation-strata × infection-strata weight matrix, and a single template model is replicated once per data stream. Split(template, W) projects the infection series reaching it through W, so it composes inside an IDModel like any other observation model: the infections come from the modelled process, not a hand-built series. One weight matrix covers the one-to-one (an identity map), many-to-one (an aggregation row summing infection strata into one stream), and many-to-many (a general matrix) infection → observation cases.
Here the renewal process supplies one infection stratum, and W maps it onto a young band, an old band, and their total:
W = reshape([0.7, 0.3, 1.0], 3, 1) # young, old, and their total
weighted = Split(LatentDelay(PoissonError(), LogNormal(1.6, 0.5)), W)
weighted_model = IDModel(renewal, weighted)
age = as_turing_model(
weighted_model, (young = missing, old = missing, total = missing), n)()
map(s -> sum(skipmissing(s)), age.generated_y_t) # simulated total per band(young = 3718, old = 1619, total = 5585)The aggregate total stream sees the summed expected infections of both bands — its expected series is exactly young .+ old.
Here the single renewal process supplied one infection stratum, broadcast through W. When the strata are genuinely separate infection processes — several distinct regions, say, each with its own latent — swap the single infection model for CombineInfections: it draws each process independently and stacks the results into the same infection-strata × time matrix Split/StrataMap already expect, so IDModel(CombineInfections([...]), Split(template, W)) maps several distinct infection processes onto streams end-to-end. For one process carried across a strata axis instead, with partially pooled per-stratum deviations, see Stratify and Partial pooling across groups. It also composes with Renewal's mixing slot, see Coupled patch models.
References
- K. Sherratt, S. Abbott, S. R. Meakin, J. Hellewell, J. D. Munday, N. Bosse, M. Jit and S. Funk. Exploring surveillance data biases when estimating the reproduction number: with insights into subpopulation transmission of COVID-19 in England. Philosophical Transactions of the Royal Society B 376, 20200283 (2021).