ComposableTuringIDModels.jl

Appearance

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.

julia
using ComposableTuringIDModels, Distributions, Random, Turing, Mooncake
using ADTypes: AutoMooncake
Random.seed!(1234)

data = IDData(gen_distribution = Gamma(6.5, 0.62))
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(data; 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: Intercept

The composed model assembles the renewal infection process and the two-stream observation model exactly like a single-stream study.

julia
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: Intercept

Passing 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.

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

julia
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 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/4hMHm/src/mcmc/hmc.jl:483
┌ Info: Found initial step size
└   ϵ = 0.00625
┌ 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/4hMHm/src/mcmc/hmc.jl:483
┌ Warning: There were 8 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/4hMHm/src/mcmc/hmc.jl:483

The 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 process; the sparse death stream is observed jointly rather than fit in isolation.

julia
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.0235  0.1947  0.0048  1746.5875  837.9152  1.0107  …      │
│     ar_init[2]   0.3121  0.1381  0.0051   763.6380  754.9844  1.0005  …      │
│     damp_AR[1]   0.7906  0.0450  0.0017   731.3182  413.0659  1.0088  …      │
│     damp_AR[2]   0.0885  0.0409  0.0018   444.4279  188.3290  0.9993  …      │
│            std   0.0662  0.0281  0.0012   532.2000  718.9306  1.0021  …      │
│         ϵ_t[1]   0.3905  0.8918  0.0246  1217.2689  610.5076  0.9991  …      │
│         ϵ_t[2]   0.4173  0.9659  0.0243  1592.9857  541.8490  0.9999  …      │
│         ϵ_t[3]   0.3590  0.9909  0.0251  1559.9846  833.3406  0.9997  …      │
│         ϵ_t[4]   0.3199  0.9727  0.0254  1462.4820  642.2763  1.0016  …      │
│         ϵ_t[5]   0.2722  0.9884  0.0253  1510.8946  875.6889  1.0064  …      │
│         ϵ_t[6]   0.1895  0.9808  0.0237  1643.8448  649.3610  0.9997  …      │
│         ϵ_t[7]   0.1506  0.9838  0.0228  1870.6447  682.1348  1.0000  …      │
│         ϵ_t[8]   0.0859  0.9993  0.0242  1702.6562  733.1849  1.0011  …      │
│         ϵ_t[9]   0.0299  0.9615  0.0277  1208.4816  671.7461  0.9993  …      │
│        ϵ_t[10]  -0.0031  1.0111  0.0222  2051.4299  621.6825  0.9994  …      │
│        ϵ_t[11]  -0.1115  0.9895  0.0243  1656.5500  627.4927  0.9995  …      │
│        ϵ_t[12]  -0.1526  0.9599  0.0257  1409.1509  765.2927  1.0001  …      │
│        ϵ_t[13]  -0.2014  0.9318  0.0240  1518.2055  873.8228  1.0044  …      │
│        ϵ_t[14]  -0.1275  0.9312  0.0182  2591.4822  773.1527  1.0015  …      │
│        ϵ_t[15]  -0.1653  0.9765  0.0262  1423.0402  722.0794  0.9999  …      │
│        ϵ_t[16]  -0.1556  0.9318  0.0220  1756.1740  792.1886  1.0009  …      │
│        ϵ_t[17]  -0.0954  1.0004  0.0271  1385.7044  585.4064  0.9992  …      │
│        ϵ_t[18]  -0.1458  0.9811  0.0273  1296.0606  739.4534  0.9991  …      │
│        ϵ_t[19]  -0.1586  0.9277  0.0169  2971.7775  676.5287  1.0064  …      │
│        ϵ_t[20]  -0.1897  1.0445  0.0274  1447.9894  665.7081  1.0099  …      │
│        ϵ_t[21]  -0.1472  0.9612  0.0250  1494.2785  806.8491  1.0032  …      │
│        ϵ_t[22]  -0.2292  1.0145  0.0320  1014.6518  515.1903  1.0007  …      │
│        ϵ_t[23]  -0.3206  0.9974  0.0224  1965.5654  709.3554  1.0039  …      │
│        ϵ_t[24]  -0.3399  0.9275  0.0204  2091.5900  697.9295  1.0003  …      │
│        ϵ_t[25]  -0.3007  1.0029  0.0256  1496.8387  593.8112  1.0112  …      │
│        ϵ_t[26]  -0.2136  0.9448  0.0227  1728.2583  693.7350  0.9995  …      │
│        ϵ_t[27]  -0.2169  0.9381  0.0272  1227.5362  613.4144  1.0014  …      │
│        ϵ_t[28]  -0.3057  0.9240  0.0249  1345.7831  644.1115  1.0000  …      │
│        ϵ_t[29]  -0.3235  0.9736  0.0263  1362.0906  750.5713  1.0112  …      │
│        ϵ_t[30]  -0.3959  0.9552  0.0227  1714.1395  669.1924  0.9995  …      │
│        ϵ_t[31]  -0.2930  0.9493  0.0212  2017.8069  717.2317  1.0080  …      │
│        ϵ_t[32]  -0.3195  0.9488  0.0214  1980.7259  621.0526  0.9991  …      │
│        ϵ_t[33]  -0.4417  0.8711  0.0232  1402.9607  678.9912  1.0011  …      │
│        ϵ_t[34]  -0.4519  0.8862  0.0217  1620.5415  790.3256  1.0063  …      │
│        ϵ_t[35]  -0.4250  0.9800  0.0309   987.3425  613.4144  0.9991  …      │
│        ϵ_t[36]  -0.3536  0.9840  0.0265  1310.3991  653.4660  1.0111  …      │
│        ϵ_t[37]  -0.2507  0.9999  0.0216  2104.6761  677.4645  1.0019  …      │
│        ϵ_t[38]  -0.1297  1.0391  0.0330   989.7077  737.8926  1.0013  …      │
│        ϵ_t[39]  -0.0169  0.9729  0.0241  1597.7865  737.4905  0.9992  …      │
│        ϵ_t[40]   0.0905  0.9626  0.0282  1164.3546  593.9812  0.9995  …      │
│        ϵ_t[41]   0.1770  0.9678  0.0273  1278.0951  670.7093  1.0006  …      │
│        ϵ_t[42]   0.2943  0.9644  0.0216  1973.4185  650.7920  0.9990  …      │
│        ϵ_t[43]   0.4344  0.9629  0.0247  1526.3535  783.1348  0.9990  …      │
│        ϵ_t[44]   0.4813  0.9421  0.0215  1910.8169  674.1631  0.9991  …      │
│        ϵ_t[45]   0.3608  0.9642  0.0215  2006.0688  677.9405  1.0013  …      │
│        ϵ_t[46]   0.3076  0.9357  0.0244  1485.0233  812.7352  1.0130  …      │
│        ϵ_t[47]   0.2463  0.9380  0.0238  1532.1462  880.2411  1.0004  …      │
│        ϵ_t[48]   0.3181  0.9503  0.0228  1737.6970  721.7165  1.0000  …      │
│        ϵ_t[49]   0.3298  0.9224  0.0261  1246.8015  844.5899  1.0021  …      │
│        ϵ_t[50]   0.3167  0.9358  0.0241  1475.8535  871.8009  0.9995  …      │
│        ϵ_t[51]   0.3626  0.9535  0.0245  1510.6200  581.1769  0.9991  …      │
│        ϵ_t[52]   0.2744  0.9105  0.0216  1761.5905  739.5609  0.9995  …      │
│        ϵ_t[53]   0.1922  0.9697  0.0247  1491.5509  697.1544  1.0059  …      │
│        ϵ_t[54]   0.0183  0.9834  0.0221  1979.1351  676.7590  0.9991  …      │
│        ϵ_t[55]  -0.0285  0.9281  0.0241  1498.5407  808.3068  0.9998  …      │
│        ϵ_t[56]  -0.0972  0.9325  0.0242  1501.8610  718.9306  0.9996  …      │
│        ϵ_t[57]  -0.0324  0.9387  0.0203  2137.4940  875.5169  0.9994  …      │
│        ϵ_t[58]   0.0231  0.9475  0.0279  1155.3321  769.2823  1.0005  …      │
│        ϵ_t[59]   0.1270  0.9485  0.0243  1523.0915  617.4658  1.0007  …      │
│        ϵ_t[60]   0.2076  0.9560  0.0185  2653.0337  622.2702  1.0006  …      │
│        ϵ_t[61]   0.2131  0.9699  0.0268  1318.4993  838.5173  0.9995  …      │
│        ϵ_t[62]   0.1961  0.9446  0.0200  2224.2826  719.7487  0.9990  …      │
│        ϵ_t[63]   0.1192  0.9311  0.0248  1421.7528  687.9418  1.0004  …      │
│        ϵ_t[64]  -0.0111  1.0040  0.0237  1804.2494  562.6525  1.0040  …      │
│        ϵ_t[65]   0.0258  0.9417  0.0233  1630.5426  764.5035  0.9992  …      │
│        ϵ_t[66]   0.0121  0.9509  0.0192  2445.4691  827.5536  0.9994  …      │
│        ϵ_t[67]   0.0006  1.0075  0.0303  1120.3813  613.7102  0.9998  …      │
│        ϵ_t[68]   0.0400  0.9743  0.0271  1288.9252  594.0827  0.9993  …      │
│   init_incide…   4.6574  0.0937  0.0026  1251.6224  780.6171  0.9992  …      │
│   cases.clust…   0.0953  0.0204  0.0007   805.1697  443.4929  1.0074  …      │
│   deaths.Asce…  -4.1928  0.1080  0.0029  1439.8290  632.4842  1.0004  …      │
│   deaths.clus…   0.1066  0.0819  0.0023   870.7337  548.7630  0.9999  …      │
╰──────────────────────────────────────────────────────────────────────────────╯

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   , and it needs no new construct and no mode flag — it is the same 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.

julia
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.48168753785353, 278.2756274314308, 290.9627856707276, 301.28182061015445, 312.3983423411687, 324.84589430964934, 335.20267973117086, 339.2005876922202  …  31.186316265782217, 28.694989824511666, 26.44369524061574, 24.485114690631903, 22.782521197898706, 21.169716185551678, 19.462217764470903, 17.66258532039226, 15.938218082601603, 14.341444599295773], deaths = [6.438898577278827, 6.52699593084055, 6.567662962367721, 6.564000144659306, 6.517047629226355, 6.424169592255845, 6.280901953934635, 6.084810397911528, 5.838057962715189, 5.547591893051596  …  1.1520649516076356, 1.0772106846356602, 1.0106802434475657, 0.9503586772567186, 0.8943779493996293, 0.8411996384131999, 0.7896930783313525, 0.739292731336043, 0.6900632656439981, 0.6424963952669397]), I_t = [124.82023240024644, 118.78877067901828, 125.1487735449749, 134.7185804215374, 158.58527604659736, 156.77778466087597, 176.86934435498142, 184.4853794947595, 201.74499799634336, 237.88281397858518  …  20.165685555311335, 18.842248244661604, 17.723028988752468, 15.751613393207471, 12.275154235021272, 11.73102792823999, 10.622872267713365, 9.386131507924796, 7.465721130334022, 5.331369642723152], Z_t = [0.5312446389953557, 0.3423504104670212, 0.2619472009124497, 0.2309590587560234, 0.32222357859830564, 0.25145249562109956, 0.3032385726661695, 0.2704627134438929, 0.28597582344222056, 0.3785932527760605  …  -0.30701422716001264, -0.2801358465279927, -0.25738877419206196, -0.3042302765638897, -0.4888920355466818, -0.46241069165440163, -0.46741840122631073, -0.4769601397674182, -0.5856294205914865, -0.8032070473681384])

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.

julia
(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.

julia
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:

julia
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 ​

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