Household models

HouseholdProcess spreads infection within households. The population is partitioned into households; within a household every infectious member can infect every susceptible household-mate, with the timing of infectious contact drawn from a contact-interval kernel (Kenah 2011). It is a structure-driven model like NetworkProcess, and like it is simulated by the Sellke construction in continuous time (the exact generative model of its pairwise likelihood) rather than by the generation-based engine; a household is a small, depleting clique rather than a fixed graph.

It lives in the companion EpiHouseholds package.

Defining a household model

The contact-interval kernel is the one required input to HouseholdProcess. sizes gives the size of each household. The disease's natural history is a progression attached with a ModelSpec.

using EpiBranch
using EpiHouseholds
using Distributions
using StableRNGs

# 300 households of four, a Weibull contact interval, a six-day infectious period
model = ModelSpec(HouseholdProcess(fill(4, 300), Weibull(1.5, 3.0));
    progression = [Transition(:recovered; from = :infection, delay = 6.0, terminal = true)])

ModelSpec(HouseholdProcess(300 households, 1200 individuals, kernel=Weibull); 0 interventions, 1 transitions)

Simulating

simulate returns a SimulationState, and linelist renders the one-row-per-case table. Each household is seeded with one index, and the outbreak spreads within it.

state = simulate(model; rng = StableRNG(1))
df = linelist(state)
(cases = size(df, 1), indexes = count(df.index))

(cases = 1200, indexes = 300)

A flexible natural history

The infectious timeline is a progression of Transitions on the ModelSpec, exactly as for BranchingProcess. A latent period is a Transition(:infectious; from = :infection, …); an infectious period is a terminal removal transition timed from the state before it. The progression's states become line-list columns, so symptom onset, testing and recovery come straight out of the simulation. The kernel times each infectious contact from the infectious window's start; with a latent period present, that from state is derived as :infectious, otherwise :infection.

clinical = ModelSpec(HouseholdProcess(fill(4, 300), Weibull(1.5, 3.0));
    progression = [
        Transition(:infectious; from = :infection, delay = LogNormal(1.2, 0.4)),  # infection → infectiousness
        Transition(:recovered; from = :infectious, delay = Gamma(6, 1),           # infectiousness → recovery
            terminal = true)])
sort(propertynames(linelist(simulate(clinical; rng = StableRNG(2)))))

13-element Vector{Symbol}:
 :chain_id
 :date_infection
 :date_infectious
 :date_outcome
 :date_recovered
 :generation
 :household
 :id
 :index
 :infectious
 :outcome
 :parent_id
 :recovered

date_infectious and date_recovered appear because the progression writes :infectious_time and :recovered_time onto each case.

The pairwise likelihood

Infections are latent: the model generates them, and the progression maps each to its observable outcomes. pairwise_surv_loglik is the contact-process density of that infection layer, which household_infections reads out of a simulation. Because the Sellke construction is the likelihood's generative model, simulate → loglikelihood is an exact round trip, so the simulated outbreak recovers the kernel.

truth = ModelSpec(HouseholdProcess(fill(4, 500), Exponential(4.0));
    progression = [Transition(:recovered; from = :infection, delay = 6.0, terminal = true)])
data = household_infections(simulate(truth; rng = StableRNG(3)), truth)

ll(scale) = pairwise_surv_loglik(Exponential(scale), data)
grid = 2.0:0.5:6.0
grid[argmax([ll(s) for s in grid])]   # ≈ the true scale, 4.0

4.0

The kernel can also be a callable (infector, susceptible) -> Distribution that takes host ids, which allows covariate models such as adults transmitting faster than children. The simulator and every pairwise_surv_loglik form that takes household data call it with the ids in that order. One callable therefore works for both simulation and fitting, and fitting the simulated outbreak recovers each of its parameters.

Fitting with Turing

When the infection layer is observed (here it comes directly from the simulation), the likelihood slots into a Turing @model. Put a prior on the log contact rate and add the pairwise log-density to the target. The household structure is fixed across draws, so compile_household_pairs captures the pair layout once and each evaluation reuses it — no per-sample rebuild:

using Turing

layout = compile_household_pairs(data)   # the fixed pair structure, compiled once

@model function household_fit(data, layout)
    logβ ~ Normal(-1, 1)                # log within-household contact rate
    Turing.@addlogprob! pairwise_surv_loglik(Exponential(1 / exp(logβ)), data, layout)
end

chain = sample(StableRNG(4), household_fit(data, layout), NUTS(), 300; progress = false)
exp(-mean(chain[:logβ]))                # posterior mean contact-interval scale, ≈ 4.0

3.8347763521889964

The plain pairwise_surv_loglik(kernel, data; external_hazard) form re-derives the pair structure (a bucketed pass over households, one susceptible-grouped row list) on every call. HouseholdPairsLayout hoists that structural work out of the gradient loop: compile_household_pairs enumerates the ordered (susceptible, infector) rows once — everything that doesn't depend on the sampled parameters — and the three-argument pairwise_surv_loglik(kernel, data, layout) then evaluates the density in two allocation-free passes, reading the (possibly augmented) times on the fly. The two forms agree up to row order.

In real data the infection times are unobserved. A household @model then augments them and conditions the observed onsets and tests through the progression's delays, with pairwise_surv_loglik supplying the contact-process density of the augmented configuration. The layout stays valid across draws as long as the household structure and the set of ever-infected hosts are fixed — only the latent times move — so it is compiled once, outside the model, and reused.