Clinical transitions

Each case follows a disease timeline after infection: symptom onset, reporting, maybe admission to hospital, and recovery or death. Each event happens a delay after an earlier one. You describe the timeline as a list of events called progression and give it to a ModelSpec alongside the transmission process. Control measures (isolation, contact tracing, vaccination) go in the model's interventions; the course of disease goes in its progression.

Each event takes its timing as either a delay or a rate. A delay is a fixed number of days or a distribution of days. A rate r gives an exponentially distributed waiting time with mean 1/r days, as in a compartmental model; only the general Transition event below takes a rate.

Built-in events

The package has four ready-made events, each timed from symptom onset by default:

• Reporting: a symptomatic case is reported after a delay, with a probability that can be below 1.

• Hospitalisation: a proportion of cases are admitted after a delay.

• Death: a case dies after a delay, with a given probability.

• Recovery: every symptomatic case gets a recovery time.

Death and recovery are final outcomes: each case ends with one of them (see Competing outcomes below). Each event has its own delay distribution and probability:

using EpiBranch
using Distributions
using StableRNGs

# incubation period, infection to onset: mean about 5 days
clinical = clinical_presentation(incubation_period = LogNormal(1.5, 0.5))

progression = [
    Reporting(delay = LogNormal(1.0, 0.3)),                         # onset to report
    Hospitalisation(delay = LogNormal(2.0, 0.5), probability = 0.2), # onset to admission
    Death(delay = LogNormal(2.5, 0.4), probability = 0.05),          # onset to death
    Recovery(delay = LogNormal(2.0, 0.4)),                           # onset to recovery
]

model = ModelSpec(
    BranchingProcess(Poisson(2.0), Exponential(5.0));  # R = 2, generation time mean 5 days
    progression = progression, attributes = clinical)

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

cases = linelist(state)
first(cases[:, [:id, :date_infection, :date_onset, :date_reporting,
    :date_admission, :outcome, :date_outcome]], 5)

5×7 DataFrame

Row	id	date_infection	date_onset	date_reporting	date_admission	outcome	date_outcome
	Int64	Date	Date?	Date?	Date?	String?	Date?
1	1	2020-01-01	2020-01-06	2020-01-10	missing	recovered	2020-01-11
2	2	2020-01-09	2020-01-11	2020-01-15	missing	recovered	2020-01-17
3	3	2020-01-18	2020-01-24	2020-01-26	missing	recovered	2020-01-29
4	4	2020-01-21	2020-01-26	2020-01-29	missing	recovered	2020-02-01
5	5	2020-01-26	2020-01-28	2020-02-01	missing	recovered	2020-02-07

LogNormal(μ, σ) takes the mean and standard deviation of the logarithm of the delay in days, not of the delay itself: LogNormal(1.0, 0.3) has a mean of about 2.8 days. Exponential(θ) has mean θ days. Each row of the line list is a case, with the date of each event it reached and its final outcome.

Reporting and Hospitalisation draw their own time independently of Death and Recovery, so a draw can fall after the case's outcome, although a report or an admission cannot happen to a person who has already died or recovered. The engine resets such an event to "did not occur" (:reported/:admitted back to false, their time back to Inf): see censor_after_outcome!.

Where admission is meant to change the outcome rather than sit alongside it, chain Death/Recovery onto Hospitalisation with from = :admission_time instead, so the outcome delay is measured from admission:

progression_post_admission = [
    Hospitalisation(delay = LogNormal(2.0, 0.5), probability = 0.2),
    Death(delay = LogNormal(1.5, 0.4), probability = 0.1, from = :admission_time),
    Recovery(delay = LogNormal(2.0, 0.4)),
]
model = ModelSpec(
    BranchingProcess(Poisson(2.0), Exponential(5.0));
    progression = progression_post_admission, attributes = clinical
)
state = simulate(model; max_cases = 200, rng = StableRNG(42))

SimulationState(cases=200, individuals=200, gen=15, active)

A case not admitted still recovers from onset as before: Death is chained from :admission_time, which that case never reaches, so it never enters that race.

Asymptomatic cases

Only cases with a symptom onset go through events timed from onset. Cases without an onset, because they are asymptomatic or never observed, are skipped. clinical_presentation makes a proportion of cases asymptomatic with prob_asymptomatic, for example clinical_presentation(incubation_period = LogNormal(1.5, 0.5), prob_asymptomatic = 0.3). The default of 0 gives every case an onset.

Probabilities and delays that depend on the person

probability and delay also accept a function of the random number generator and the person, (rng, ind) -> value. ind.state[:age] is the person's age, set by demographics. This covers age-specific case fatality, delays that differ by risk group, and reporting that differs between groups. c ? a : b is Julia for "if c then a, otherwise b", like ifelse in R.

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

# Probability of death among symptomatic cases: 30% for 80+, 2% otherwise.
death_age = Death(
    delay = LogNormal(2.5, 0.4),
    probability = (rng, ind) -> ind.state[:age] >= 80 ? 0.3 : 0.02,
)

# Onset to admission: 1 day for under-30s, 5 days otherwise.
hosp_age = Hospitalisation(
    delay = (rng, ind) -> ind.state[:age] < 30 ? 1.0 : 5.0,
    probability = 0.2,
)

progression = [hosp_age, death_age, Recovery(delay = LogNormal(2.0, 0.4))]
model = ModelSpec(BranchingProcess(Poisson(2.0), Exponential(5.0));
    progression = progression, attributes = attrs)

rng = StableRNG(42)
state = simulate(model; max_cases = 2000, rng = rng)

n_died = count(ind -> ind.state[:outcome] == :died, state.individuals)
n_died_80plus = count(state.individuals) do ind
    ind.state[:outcome] == :died && ind.state[:age] >= 80
end
println("Deaths: $n_died, of which 80+: $n_died_80plus")

Deaths: 24, of which 80+: 17

People aged 80 and over are about one in nine of the population here, but their probability of death is 15 times higher, so they make up most of the deaths: 17 of 24 here. With so few deaths, the printed split varies between runs. Age is one example: any characteristic you give people through attributes (risk group, comorbidity, vaccination status) can be used the same way.

Events that depend on an earlier event

Some events happen only after another: only reported cases are admitted, only tested cases are treated. Make the probability of the later event 0 for everyone who has not had the earlier one. List the events in the order they happen, because an event can only depend on events listed before it.

# Only reported cases can be admitted; 20% of them are.
reported_hosp = Hospitalisation(
    delay = LogNormal(2.0, 0.5),
    probability = (rng, ind) -> get(ind.state, :reported, false) ? 0.2 : 0.0,
)

progression = [
    Reporting(delay = LogNormal(1.0, 0.3), probability = 0.5),
    reported_hosp,
]
model = ModelSpec(BranchingProcess(Poisson(2.0), Exponential(5.0));
    progression = progression, attributes = clinical)

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

n_admitted = count(ind -> ind.state[:admitted], state.individuals)
n_admitted_unreported = count(ind -> ind.state[:admitted] && !ind.state[:reported],
    state.individuals)
println("Admitted: $n_admitted, of which not reported: $n_admitted_unreported")

Admitted: 24, of which not reported: 0

No unreported case was admitted. get(ind.state, :reported, false) reads whether the case was reported, treating it as not reported if the value is missing. Conditions can be combined, for example admitting only cases that were reported and are not vaccinated.

Measuring delays from an earlier event

By default each built-in event is timed from symptom onset. The from argument times it from another event instead. Suppose cases are reported after a positive test, with the reporting delay measured from the test and not from onset. The general Transition event describes the test: Transition(:tested; from = :onset, ...) means that a case is tested a delay after onset, here with probability 0.9 (the proportion of symptomatic cases who get a positive test). It records whether the case was tested as :tested, and when as :tested_time. Reporting is then timed from the test with from = :tested_time:

testing = Transition(:tested; from = :onset, delay = LogNormal(0.5, 0.3),
    probability = 0.9)
reporting_post_test = Reporting(delay = LogNormal(0.0, 0.2), from = :tested_time)

progression = [testing, reporting_post_test]
model = ModelSpec(BranchingProcess(Poisson(2.0), Exponential(5.0));
    progression = progression, attributes = clinical)

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

first(linelist(state)[:, [:id, :date_onset, :tested, :date_tested, :date_reporting]], 5)

5×5 DataFrame

Row	id	date_onset	tested	date_tested	date_reporting
	Int64	Date?	Bool	Date?	Date?
1	1	2020-01-06	false	missing	missing
2	2	2020-01-11	true	2020-01-13	2020-01-14
3	3	2020-01-12	true	2020-01-13	2020-01-14
4	4	2020-01-11	true	2020-01-13	2020-01-14
5	5	2020-01-21	true	2020-01-22	2020-01-23

Cases without a test have no report. In each row, the test date follows onset and the report date follows the test.

from accepts either the event or its recorded time

Transition, Reporting, Hospitalisation, Death and Recovery resolve from the same way: from = :onset and from = :onset_time both read :onset_time; from = :infection reads the case's infection time.

The general Transition covers any event in the timeline: it takes the name of the event, from (the earlier event it is timed from, by default infection), a delay or rate, a probability, and terminal = true for a final outcome.

Delays from infection

If you do not model symptom onset, time events from infection instead, and no clinical_presentation is needed. The built-in events take from as a function of the person:

Reporting(delay = LogNormal(2.0, 0.3), from = ind -> ind.infection_time)

Transition(:reported; from = :infection, ...) also times an event from infection (its default), but it records :reported_time. The line list then has a date_reported column in place of date_reporting, and weekly_incidence(state; by = :reporting) counts no cases. Use Reporting when you want the standard reporting dates. Each event has its own from: a timeline can admit cases a delay after onset and time their outcome from admission.

Combining with multi-type models and population characteristics

The disease timeline is separate from the transmission model. A multi-type model with age or risk groups can have outcomes that differ by type. A rule can look up the case's type, ind.state[:type] (1 for the first type, 2 for the second, in the order of the next-generation matrix), as well as its age:

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

# Probability of death depends on both type and age.
death_type_age = Death(
    delay = LogNormal(1.5, 0.4),  # onset to death: mean about 5 days
    probability = (rng, ind) -> begin
        base = ind.state[:age] >= 65 ? 0.15 : 0.01
        ind.state[:type] == 1 ? 0.5 * base : base  # children: half the risk
    end,
)

# Onset to admission depends on type.
hosp_type = Hospitalisation(
    delay = (rng, ind) -> ind.state[:type] == 1 ? 1.0 : 3.0,
    probability = 0.2,
)

progression = [
    hosp_type,
    death_type_age,
    Recovery(delay = LogNormal(2.0, 0.4)),
]

# Two types, children and adults. Columns are infectors, rows infectees:
# a child causes 2.0 child and 0.8 adult cases, an adult 0.5 child and
# 1.5 adult cases. Secondary cases are negative binomial with dispersion k = 0.5.
multitype = ModelSpec(
    BranchingProcess(
        [2.0 0.5; 0.8 1.5],
        R -> NegBin(R, 0.5),
        LogNormal(1.6, 0.5),
        type_labels = ["children", "adults"]);
    progression = progression, attributes = attrs_age)

rng = StableRNG(3)
state = simulate(multitype; max_cases = 500, rng = rng)

n_died_kids = count(state.individuals) do ind
    ind.state[:outcome] == :died && ind.state[:type] == 1
end
n_died_adults = count(state.individuals) do ind
    ind.state[:outcome] == :died && ind.state[:type] == 2
end
println("Deaths. Children: $n_died_kids. Adults: $n_died_adults.")

Deaths. Children: 5. Adults: 9.

Children and adults make up similar numbers of cases in this outbreak. The halved risk among children shows as fewer deaths, although the counts are small and noisy. Death here comes sooner after onset than in the first example, which leaves fewer of those who would die to recover first (see Competing outcomes). Relative susceptibility and infectiousness can also vary between people, through transmission_traits, and outcome rules can use any of these characteristics.

Competing outcomes

Death and recovery compete: each case ends with whichever final outcome happens first, the situation survival analysis calls competing risks. The outcome is recorded as :outcome, and its time as :outcome_time.

This affects what the probability of Death means. It is the probability that a case would die if it did not recover first. When recovery tends to come sooner than death, as in the first example on this page, some of the cases that would have died recover first. The proportion of cases dying is then below the probability given.

Further final outcomes join the competition through the general Transition with terminal = true. Here 10% of cases are lost to follow-up, a mean of about 5 days after onset, unless they die or recover first:

progression = [
    Death(delay = LogNormal(2.5, 0.4), probability = 0.05),
    Recovery(delay = LogNormal(2.0, 0.4)),
    Transition(:lost; from = :onset, delay = LogNormal(1.5, 0.5),
        probability = 0.1, terminal = true),
]
model = ModelSpec(BranchingProcess(Poisson(2.0), Exponential(5.0));
    progression = progression, attributes = clinical)

rng = StableRNG(42)
state = simulate(model; max_cases = 2000, rng = rng)

outcomes = [ind.state[:outcome] for ind in state.individuals if haskey(ind.state, :outcome)]
println("Outcome counts: ",
    (died = count(==(:died), outcomes),
     recovered = count(==(:recovered), outcomes),
     lost = count(==(:lost), outcomes)))

Outcome counts: (died = 16, recovered = 1827, lost = 157)

Most cases recover. The deaths and losses to follow-up, 16 and 157 of 2,000 cases, are fewer than 5% and 10% of cases, because recovery often comes first. The same approach gives admission to intensive care after admission to hospital, or outcomes that depend on treatment.

An exact case fatality ratio

To fix the proportion of symptomatic cases who die, the outcomes must be mutually exclusive. If death and recovery each happened independently, with probabilities 0.05 and 0.95, some cases would get both and some neither. The proportion dying would then differ from 5%. exclusive_probabilities makes the outcomes mutually exclusive: each case has exactly one of them, with the probabilities given. Recovery always happens and has no probability to set. The recovery outcome is therefore a general Transition with its own probability:

death_p, recovered_p = exclusive_probabilities([0.05, 0.95])
exact_cfr = [
    Death(delay = LogNormal(2.5, 0.4), probability = death_p),
    Transition(:recovered, from = :onset, delay = LogNormal(2.0, 0.4),
        probability = recovered_p, terminal = true),
]
model = ModelSpec(BranchingProcess(Poisson(2.0), Exponential(5.0));
    progression = exact_cfr, attributes = clinical)

rng = StableRNG(42)
state = simulate(model; max_cases = 300, rng = rng)
symptomatic = [ind for ind in state.individuals if !isnan(onset_time(ind))]
n_missing = count(ind -> !haskey(ind.state, :outcome), symptomatic)
n_died = count(ind -> get(ind.state, :outcome, nothing) == :died, symptomatic)
println("Symptomatic cases: ", length(symptomatic), ", missing an outcome: ", n_missing)
println("Died: ", n_died, " of ", length(symptomatic))

Symptomatic cases: 300, missing an outcome: 0
Died: 12 of 300

Every symptomatic case has an outcome. The proportion who died is close to 5%, with the difference due to chance.

ModelSpec warns when every final outcome has a fixed probability below 1, because some cases could then end without an outcome. Probabilities from exclusive_probabilities do not trigger the warning.

Likelihood of the clinical timeline

To fit reporting or outcome delays to observed cases, you need the likelihood of each case's timeline. progression_loglik gives the log-likelihood of the simulated cases' timelines under the progression: for each event, the probability that it did or did not happen and, when it happened, the probability density of its delay. It takes the simulated outbreak that simulate returned:

progression = [
    Reporting(delay = LogNormal(1.0, 0.3), probability = 0.7),
    Recovery(delay = LogNormal(2.0, 0.4)),
]
model = ModelSpec(BranchingProcess(Poisson(2.0), Exponential(5.0));
    progression = progression, attributes = clinical)

rng = StableRNG(42)
state = simulate(model; max_cases = 200, rng = rng)
progression_loglik(model, state)

-769.1002759016546

The result is a log-likelihood: comparing it across parameter values shows which fit the timelines better, with higher values fitting better. pairwise_surv_loglik gives the likelihood of who was infected and when, summed over everyone who could have infected each case. The two added together are the log-likelihood of the whole outbreak, infection times and clinical timelines, including the parts that are not directly observed. Cases that never reached the event an event is timed from (an asymptomatic case for an event timed from onset) contribute nothing.

Delays that depend on the person cannot be evaluated

A delay given as a function (rng, ind) -> ... can be simulated but has no density, and progression_loglik throws an error on it. To use the likelihood, give each delay as a distribution or a fixed number of days.

To write your own kind of event, see Custom clinical transitions, and The likelihood of a custom transition for its likelihood.