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 case study 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 case study, 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_priors = [truncated(Normal(0.8, 0.05), 0, 1),
truncated(Normal(0.1, 0.05), 0, 1)],
init_priors = [Normal(0.0, 0.2), Normal(0.0, 0.2)],
ϵ_t = HierarchicalNormal(std_prior = HalfNormal(0.1)))
renewal = Renewal(gen_distribution = Gamma(6.5, 0.62);
rt = latent, initialisation_prior = Normal(log(100.0), 0.1))
cases = LatentDelay(
Ascertainment(NegativeBinomialError(cluster_factor_prior = 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_prior = 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 case studies, 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.9; adtype = AutoMooncake(; config = nothing)),
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.025
┌ Warning: There were 14 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
┌ 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 17 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.0299 0.1914 0.0051 1379.4717 717.2011 0.9999 … │
│ ar_init[2] 0.3000 0.1488 0.0068 501.5587 529.9798 1.0028 … │
│ damp_AR[1] 0.7890 0.0451 0.0014 1098.2883 563.3529 1.0025 … │
│ damp_AR[2] 0.0877 0.0415 0.0018 506.9129 358.2660 1.0017 … │
│ std 0.0670 0.0291 0.0016 356.6384 566.1703 1.0019 … │
│ ϵ_t[1] 0.4400 0.9208 0.0223 1713.3410 719.3022 1.0019 … │
│ ϵ_t[2] 0.3713 0.9752 0.0261 1396.8584 720.3311 0.9995 … │
│ ϵ_t[3] 0.3589 0.9728 0.0249 1545.6705 817.9909 1.0012 … │
│ ϵ_t[4] 0.3261 0.9883 0.0268 1334.9201 656.6052 1.0022 … │
│ ϵ_t[5] 0.3029 0.9494 0.0230 1694.5900 766.6262 1.0105 … │
│ ϵ_t[6] 0.2683 0.9589 0.0223 1846.9847 766.5140 0.9991 … │
│ ϵ_t[7] 0.1790 1.0927 0.0291 1420.1815 627.9748 0.9990 … │
│ ϵ_t[8] 0.0891 1.0169 0.0288 1235.3927 599.3529 1.0013 … │
│ ϵ_t[9] 0.0822 1.0208 0.0334 975.6333 390.9054 0.9990 … │
│ ϵ_t[10] 0.0429 0.9289 0.0255 1335.7093 647.1558 0.9998 … │
│ ϵ_t[11] -0.0497 0.9908 0.0244 1667.8398 788.7719 1.0003 … │
│ ϵ_t[12] -0.1298 0.9535 0.0249 1475.1384 763.2700 0.9992 … │
│ ϵ_t[13] -0.2068 0.9522 0.0309 953.9643 424.0639 0.9998 … │
│ ϵ_t[14] -0.1529 0.8976 0.0253 1266.2370 589.4294 1.0003 … │
│ ϵ_t[15] -0.1337 0.9607 0.0291 1077.1187 610.4288 1.0010 … │
│ ϵ_t[16] -0.1425 0.9840 0.0270 1333.3419 673.3429 1.0004 … │
│ ϵ_t[17] -0.1582 0.9543 0.0253 1343.5090 445.2859 0.9998 … │
│ ϵ_t[18] -0.1470 0.9890 0.0288 1189.3723 679.6304 1.0011 … │
│ ϵ_t[19] -0.1912 0.9822 0.0324 927.9244 652.9587 1.0006 … │
│ ϵ_t[20] -0.1418 0.9867 0.0242 1631.6577 597.9713 0.9997 … │
│ ϵ_t[21] -0.1476 0.9579 0.0387 615.1095 829.2486 1.0041 … │
│ ϵ_t[22] -0.2493 0.9631 0.0226 1825.0371 812.7352 0.9992 … │
│ ϵ_t[23] -0.3183 0.9887 0.0294 1124.7365 690.2574 1.0008 … │
│ ϵ_t[24] -0.3382 0.8866 0.0256 1229.6665 782.4686 0.9993 … │
│ ϵ_t[25] -0.2524 0.9810 0.0305 1019.1349 621.7398 1.0016 … │
│ ϵ_t[26] -0.2186 0.9725 0.0305 1022.6145 813.5746 1.0006 … │
│ ϵ_t[27] -0.2207 0.9801 0.0252 1520.2426 615.3933 0.9996 … │
│ ϵ_t[28] -0.3378 0.9328 0.0246 1464.5234 659.6991 1.0201 … │
│ ϵ_t[29] -0.3120 0.9471 0.0255 1365.1400 594.7144 1.0008 … │
│ ϵ_t[30] -0.3460 0.9586 0.0269 1253.0186 555.2098 1.0029 … │
│ ϵ_t[31] -0.3378 0.9257 0.0257 1376.9714 597.9713 0.9995 … │
│ ϵ_t[32] -0.3360 0.9603 0.0322 904.5532 547.6946 0.9993 … │
│ ϵ_t[33] -0.4409 0.9794 0.0297 1095.0442 590.3806 0.9994 … │
│ ϵ_t[34] -0.4367 0.9851 0.0263 1439.5999 761.5980 0.9998 … │
│ ϵ_t[35] -0.4678 0.9482 0.0310 982.7334 595.7440 0.9998 … │
│ ϵ_t[36] -0.3568 0.9314 0.0267 1199.3129 840.5052 1.0058 … │
│ ϵ_t[37] -0.2224 0.9994 0.0290 1174.3821 564.1218 1.0018 … │
│ ϵ_t[38] -0.1496 0.9396 0.0282 1089.7641 397.5262 1.0000 … │
│ ϵ_t[39] 0.0025 0.9493 0.0279 1156.0119 683.4720 1.0008 … │
│ ϵ_t[40] 0.0809 0.9894 0.0267 1375.1235 785.8980 1.0026 … │
│ ϵ_t[41] 0.1439 0.9775 0.0252 1479.8596 696.7916 1.0003 … │
│ ϵ_t[42] 0.2897 0.9189 0.0259 1265.1412 302.4636 1.0015 … │
│ ϵ_t[43] 0.4260 0.9923 0.0262 1427.4933 626.0803 1.0000 … │
│ ϵ_t[44] 0.4613 0.9430 0.0265 1271.2586 741.3550 1.0014 … │
│ ϵ_t[45] 0.4261 1.0262 0.0323 996.1556 629.9894 1.0000 … │
│ ϵ_t[46] 0.2985 0.9771 0.0240 1635.5206 848.7930 0.9992 … │
│ ϵ_t[47] 0.2674 0.9142 0.0250 1340.5378 652.8045 0.9994 … │
│ ϵ_t[48] 0.2594 0.9940 0.0293 1160.5614 715.5603 0.9991 … │
│ ϵ_t[49] 0.3086 0.9896 0.0272 1318.9681 659.4834 1.0020 … │
│ ϵ_t[50] 0.3615 0.8960 0.0238 1406.7856 650.2854 1.0004 … │
│ ϵ_t[51] 0.3989 0.9827 0.0272 1263.2953 638.0945 0.9993 … │
│ ϵ_t[52] 0.2956 0.9597 0.0266 1299.7739 518.4951 0.9995 … │
│ ϵ_t[53] 0.1301 0.8715 0.0227 1453.4729 740.0597 1.0011 … │
│ ϵ_t[54] 0.0273 0.9302 0.0284 1084.7308 672.2245 1.0007 … │
│ ϵ_t[55] -0.0203 0.9665 0.0259 1467.3667 713.9182 0.9994 … │
│ ϵ_t[56] -0.0510 0.9723 0.0278 1211.7340 596.8552 1.0012 … │
│ ϵ_t[57] -0.1623 0.9685 0.0303 1029.9201 717.4876 1.0003 … │
│ ϵ_t[58] -0.0448 0.9812 0.0287 1158.1947 692.5252 1.0017 … │
│ ϵ_t[59] 0.0709 0.9893 0.0250 1553.0066 698.4482 1.0046 … │
│ ϵ_t[60] 0.2735 0.9514 0.0282 1152.5013 777.4144 1.0007 … │
│ ϵ_t[61] 0.2273 0.9918 0.0249 1615.5806 703.1755 0.9996 … │
│ ϵ_t[62] 0.1962 0.9492 0.0263 1316.0725 663.3716 1.0001 … │
│ ϵ_t[63] 0.0972 0.9184 0.0251 1347.0599 766.6262 0.9998 … │
│ ϵ_t[64] 0.0062 0.9461 0.0280 1114.5236 622.5472 1.0006 … │
│ ϵ_t[65] 0.0238 1.0002 0.0283 1261.3252 639.4779 1.0005 … │
│ ϵ_t[66] 0.0069 1.0140 0.0264 1476.3565 715.2958 1.0031 … │
│ ϵ_t[67] -0.0030 1.0170 0.0306 1107.9496 565.5722 1.0018 … │
│ ϵ_t[68] -0.0277 0.9518 0.0271 1255.2761 683.6879 0.9996 … │
│ init_incide… 4.6630 0.0938 0.0026 1250.0059 771.6962 0.9990 … │
│ cases.clust… 0.0943 0.0201 0.0006 1057.0063 413.2729 1.0010 … │
│ deaths.Asce… -4.1920 0.1081 0.0025 1800.7663 573.8593 1.0009 … │
│ deaths.clus… 0.1060 0.0749 0.0026 524.8919 269.5494 1.0047 … │
╰──────────────────────────────────────────────────────────────────────────────╯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_prior = HalfNormal(0.1)),
deaths = LatentDelay( # case→death delay
Ascertainment(NegativeBinomialError(cluster_factor_prior = 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_prior = HalfNormal(0.1)),
FixedIntercept(log(0.7))), LogNormal(1.5, 0.4)),
old = LatentDelay(
Ascertainment(NegativeBinomialError(cluster_factor_prior = 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. Swapping the single infection stratum for an infection-strata × time matrix (one row per genuinely distinct infection process) and the weights for estimated ones is the seam a partially-observed or cross-classified reporting structure grows from; modelling those separate infection strata jointly is future work (#45).
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).