Line lists and contacts

What would the surveillance data from a simulated outbreak look like? linelist turns a simulation into a line list: a table (a DataFrame, the Julia counterpart of an R data.frame) with one row per case, in the format of the simulist R package. contacts gives the matching table of contacts, infected or not.

Every line list has the columns id, parent_id (the infector's id), generation, chain_id and date_infection. Anything else the simulation records about a case, such as symptom onset, age or outcome, becomes a column as well. A recorded time whose name ends in _time becomes a date column, so onset_time appears as date_onset. A model with isolation or quarantine adds isolated and date_isolation, and when they last a finite time (a finite duration) also date_isolation_release, the date the person was released. To get a new column, record that value on each case during the simulation, for example as a population characteristic (see the clinical transitions tutorial).

Line list

The model below has a negative binomial offspring distribution with R = 1.5 and dispersion k = 0.5 (NegBin(R, k); smaller k means more superspreading) and a log-normal generation time. All delays are in days. LogNormal(μ, σ) takes the mean and standard deviation of the log: LogNormal(1.5, 0.5) is an incubation period with a median of about 4.5 days. Exponential(θ) has mean θ.

After symptom onset each case is reported (mean 3 days after onset), admitted to hospital with probability 0.2, and is due to die with probability 0.05. Every case also gets a recovery time. Death and recovery are competing outcomes: whichever comes first ends the case's course.

using EpiBranch
using Distributions
using DataFrames
using Dates
using StableRNGs

attrs = clinical_presentation(incubation_period = LogNormal(1.5, 0.5))

progression = [
    Reporting(delay = Exponential(3.0)),
    Hospitalisation(delay = Exponential(5.0), probability = 0.2),
    Death(delay = Exponential(14.0), probability = 0.05),
    Recovery(delay = Exponential(14.0)),
]

# Secondary cases: NegBin(R = 1.5, k = 0.5); generation time: LogNormal (days)
model = ModelSpec(BranchingProcess(NegBin(1.5, 0.5), LogNormal(1.6, 0.5));
    progression = progression, attributes = attrs)

rng = StableRNG(42)  # a fixed seed makes the results reproducible
state = simulate(model; condition = 50:200, max_cases = 200, rng = rng)

ll = linelist(state; reference_date = Date(2024, 1, 1))
first(ll, 5)

5×17 DataFrame

Row	id	parent_id	generation	chain_id	date_infection	date_admission	date_death_candidate	date_onset	date_outcome	date_recovery_candidate	date_reporting	admitted	asymptomatic	incubation_period	index	outcome	reported
	Int64	Int64	Int64	Int64	Date	Date?	Date?	Date?	Date?	Date?	Date?	Bool	Bool	Float64	Bool	String?	Bool
1	1	0	0	1	2024-01-01	missing	missing	2024-01-09	2024-01-14	2024-01-14	2024-01-11	false	false	8.83045	true	recovered	true
2	2	1	1	1	2024-01-08	2024-01-13	missing	2024-01-12	2024-01-24	2024-01-24	2024-01-13	true	false	4.00131	false	recovered	true
3	3	2	2	1	2024-01-14	missing	missing	2024-01-16	2024-01-21	2024-01-21	2024-01-17	false	false	2.06682	false	recovered	true
4	4	2	2	1	2024-01-15	missing	missing	2024-01-21	2024-02-25	2024-02-25	2024-01-25	false	false	6.56388	false	recovered	true
5	5	3	3	1	2024-01-17	missing	missing	2024-01-24	2024-01-26	2024-01-26	missing	false	false	6.65169	false	recovered	false

condition = 50:200 repeats the simulation until it produces an outbreak of between 50 and 200 cases (both ends included), and max_cases stops an outbreak that keeps growing. Conditioning on size mimics an observed outbreak of known size, but the outbreaks you get are a selected subset: their R and timing are not representative of all outbreaks the model produces. reference_date is the calendar date of time 0, from which all dates are counted.

A column appears only if the simulation recorded that information. Drop the Hospitalisation step and date_admission disappears. Drop clinical_presentation and there is no symptom onset. The reporting, admission and outcome delays are measured from onset, so simulate then stops with an error asking for clinical_presentation.

Demographics

demographics assigns each case an age (in years) and a sex when it is created. They appear in the line list as age and sex columns:

attrs_demo = [
    clinical_presentation(incubation_period = LogNormal(1.5, 0.5)),
    demographics(age_distribution = Normal(40, 15), prob_female = 0.55),
]

model = ModelSpec(BranchingProcess(NegBin(1.5, 0.5), LogNormal(1.6, 0.5));
    progression = progression, attributes = attrs_demo)

rng = StableRNG(42)
state = simulate(model; condition = 50:200, max_cases = 200, rng = rng)

ll = linelist(state; reference_date = Date(2024, 1, 1))
println("Age range: $(minimum(ll.age)) - $(maximum(ll.age))")
println("Female: $(round(count(==("female"), ll.sex) / nrow(ll) * 100, digits=1))%")

Age range: 0 - 80
Female: 58.0%

The proportion female should be close to the 55% asked for.

Age-specific case fatality risk

To make the case fatality risk depend on age, give Death a function that returns each case's probability of death. The function below uses risks of 0.1% under 15, 1% from 15 to 64 and 15% at 65 and over:

attrs_demo = [
    clinical_presentation(incubation_period = LogNormal(1.5, 0.5)),
    demographics(age_distribution = Uniform(0, 90)),
]

function cfr_by_age(rng, ind)
    age = ind.state[:age]
    if age < 15
        return 0.001
    elseif age < 65
        return 0.01
    else
        return 0.15
    end
end

age_stratified = [
    Death(delay = Exponential(14.0), probability = cfr_by_age),
    Recovery(delay = Exponential(14.0)),
]

model = ModelSpec(BranchingProcess(NegBin(1.5, 0.5), LogNormal(1.6, 0.5));
    progression = age_stratified, attributes = attrs_demo)

rng = StableRNG(42)
state = simulate(model; condition = 100:500, max_cases = 500, rng = rng)

ll = linelist(state; reference_date = Date(2024, 1, 1))

# Cases and deaths in each age band
for (lo, hi) in [(0, 14), (15, 64), (65, 90)]
    group = filter(row -> lo <= row.age <= hi, ll)
    n_died = count(==("died"), group.outcome)
    pct = nrow(group) > 0 ? round(n_died / nrow(group) * 100, digits=1) : 0.0
    println("Age $lo-$hi: $(nrow(group)) cases, $n_died deaths ($pct%)")
end

Age 0-14: 99 cases, 0 deaths (0.0%)
Age 15-64: 258 cases, 2 deaths (0.8%)
Age 65-90: 143 cases, 10 deaths (7.0%)

EpiBranch calls the function with two arguments: the random number generator rng and the case ind. cfr_by_age does not use rng, but it must accept it. ind.state[:age] is the age that demographics gave the case.

The proportion dying rises with age, but in the oldest band it falls well short of 15%.

The probability of death is not the case fatality risk

probability on Death is the chance that a case is due to die. Recovery still competes with it: a case due to die whose recovery time comes first recovers. Here death and recovery have the same delay distribution. About half the cases due to die therefore recover, and the proportion who die is about half of probability. For a fixed case fatality risk, use exclusive_probabilities to make death and recovery mutually exclusive, as shown in the clinical transitions tutorial.

The same approach works for any characteristic you assign to cases, such as a risk group or a comorbidity. See the clinical transitions tutorial for the other steps a case can go through.

The whole population

linelist lists cases only by default. Pass infected_only = false to list everyone in the population, as needed for a test-negative design, an attack rate by exposure, or an exposed/unexposed comparison. This fits models that simulate a fixed population from the start: random mixing (HomogeneousProcess), networks (NetworkProcess) and households (HouseholdProcess).

The model below has 200 people mixing at random. transmission_rate is the rate per day at which each infectious person makes infectious contacts, spread evenly over the population, and each infected person recovers after an exponentially distributed time with a mean of 5 days, which ends their infectiousness:

pool = ModelSpec(HomogeneousProcess(; transmission_rate = 0.6, population_size = 200);
    progression = [Transition(:recovered; from = :infection, delay = Exponential(5.0),
        terminal = true)],
    attributes = attrs)

pool_state = simulate(pool; n_initial = 2, rng = StableRNG(1))

pop = linelist(pool_state; reference_date = Date(2024, 1, 1), infected_only = false)
println("Population: $(nrow(pop)), infected: $(count(pop.infected))")
first(pop, 5)

5×14 DataFrame

Row	id	parent_id	generation	chain_id	infected	date_infection	date_onset	date_outcome	date_recovered	asymptomatic	incubation_period	index	outcome	recovered
	Int64	Int64	Int64	Int64	Bool	Date?	Date?	Date?	Date?	Bool	Float64	Bool?	String?	Bool?
1	1	153	10	16	true	2024-01-15	2024-01-20	2024-01-21	2024-01-21	false	4.70786	false	recovered	true
2	2	198	8	16	true	2024-01-15	2024-01-18	2024-01-25	2024-01-25	false	2.93948	false	recovered	true
3	3	68	8	16	true	2024-01-19	2024-01-23	2024-01-19	2024-01-19	false	3.9161	false	recovered	true
4	4	79	7	16	true	2024-01-16	2024-01-22	2024-01-27	2024-01-27	false	6.28259	false	recovered	true
5	5	197	5	16	true	2024-01-11	2024-01-16	2024-01-22	2024-01-22	false	5.07062	false	recovered	true

In a BranchingProcess model the rows are the cases plus the contacts they exposed who were not infected.

An uninfected row has missing for date_infection and for every date that follows from an infection: date_onset, the reporting, admission and outcome dates, and any date from your own _time values. Dates of events that can happen to a person whether or not they are infected are kept:

• date_trace, when the contact was traced;

• date_vaccination and date_immunity;

• date_isolation, when the contact was quarantined on being traced, and date_isolation_release, when that quarantine ended. If a traced contact was later due to be isolated at what would have been their symptom onset, that isolation is not shown (they were never infected, so never had an onset): both columns give the earlier quarantine if there was one, and are missing otherwise.

Columns that are not dates, such as asymptomatic, traced or vaccinated, are shown unchanged.

Contacts table

contacts returns one row per contact: the case (from), the contact (to), whether the contact was infected (infected), and the contact's generation and infection date:

ct = contacts(state; reference_date = Date(2024, 1, 1))
println("Total: $(nrow(ct)), Infected: $(count(ct.infected)), Not infected: $(nrow(ct) - count(ct.infected))")
first(ct, 5)

5×6 DataFrame

Row	from	to	infected	generation	infection_time	date_infection
	Int64	Int64	Bool	Int64	Float64	Date
1	1	2	true	1	2.48059	2024-01-03
2	1	3	true	1	3.02382	2024-01-04
3	1	6	true	1	6.60729	2024-01-07
4	1	7	true	1	9.02852	2024-01-10
5	2	4	true	2	5.40132	2024-01-06

Contacts who were not infected appear when something stopped transmission, such as isolation, quarantine or vaccination. This model has no interventions: every contact was infected.

Fixed-population models

For models that simulate a fixed population (HomogeneousProcess, HouseholdProcess, MixingProcess, NetworkProcess, RoutedNetwork), contacts lists every infector-infectee pair the race settled, but has no row for a household-mate or network neighbour who was exposed without being infected: these models keep no record of that.