API reference

This page documents the functions and types of EpiBranch, grouped by what they are used for. Most analyses need only a few of them:

describe transmission with BranchingProcess and combine it with the natural history, interventions and population characteristics in a ModelSpec;
add control measures such as Isolation, ContactTracing and RingVaccination;
run outbreaks with simulate and summarise them with containment_probability, linelist or chain_statistics;
calculate exact results such as extinction_probability and chain_size_distribution, and fit models to data with loglikelihood.

Delays and times are in whatever unit the delay distributions use; the documentation uses days throughout. Names written as EpiBranch.name are not made available by using EpiBranch and have to be written in full. The network and household names, such as NetworkProcess and HouseholdProcess, come from the satellite packages and need using EpiNetwork or using EpiHouseholds. A function whose name ends in ! changes one of its arguments in place. For the hooks the engine dispatches on, that argument is the simulation state or the person being resolved, and the first argument is the intervention, step or action the method is written for. The last section, For extension authors, lists what you need only to write a new model, intervention or observation model.

Transmission models

A transmission model says who infects whom and when. BranchingProcess is the offspring-distribution model used by epichains and ringbp, optionally with several types of case. HomogeneousProcess is a stochastic SIR/SEIR epidemic in a closed population. MixingProcess is the same in a population that mixes by group. NetworkProcess, RoutedNetwork and HouseholdProcess spread infection over a fixed set of contacts. ModelSpec adds the rest of the scenario.

EpiBranch.BranchingProcess Type

julia

struct BranchingProcess{W<:Tuple, P, L} <: TransmissionModel

julia

BranchingProcess(offspring, generation_time; population_size = NoPopulation())
BranchingProcess(offspring)
BranchingProcess(windows::Infectiousness...)

A stochastic branching process: each case independently infects a random number of others, drawn from the offspring distribution offspring, and the outbreak grows as a tree from the index cases. A negative binomial, NegBin(R, k) with mean R and dispersion k, is the usual choice; smaller k means more superspreading.

generation_time is the distribution of the time in days from a case's infection to the infection of each of its secondary cases. Leave it out to study only chain sizes and lengths, without timing. population_size sets a finite population: each potential infection succeeds with probability equal to the share of people not yet infected, so transmission slows as cases accumulate (unlimited by default).

The transmission model describes only who infects whom and when. Add the natural history, interventions, population characteristics and reporting with a ModelSpec.

Examples

julia

# R = 2.5, k = 0.16; generation time log-normal with log-mean 1.6 and
# log-sd 0.5 (mean about 5.6 days)
BranchingProcess(NegBin(2.5, 0.16), LogNormal(1.6, 0.5))

# no timing: chain sizes and lengths only
BranchingProcess(NegBin(0.8, 0.5))

# two types (e.g. children and adults): M[i, j] is the mean number of
# type-i cases infected by one type-j case
M = [1.2 0.4; 0.3 0.9]
BranchingProcess(M, R -> NegBin(R, 0.16), LogNormal(1.6, 0.5))

# add isolation with a ModelSpec
ModelSpec(BranchingProcess(NegBin(2.5, 0.16), LogNormal(1.6, 0.5));
    attributes = clinical_presentation(incubation_period = LogNormal(1.6, 0.5)),
    interventions = [Isolation(onset_to_isolation_delay = Exponential(2.0), duration = 7.0)])

For transmission by several routes, or starting and stopping at points in the natural history, build the process from Infectiousness windows. The forms above make a single one that starts at infection.

Fields

infectiousness::Tuple
population_size::Any
n_types::Int64
type_labels::Any

EpiBranch.Infectiousness Type

julia

struct Infectiousness{O, F, U, K}

julia

Infectiousness(offspring; from = :infection, until = (), kernel = NoGenerationTime())

When, and how much, a case transmits by one route. offspring is the distribution of the number of people infected by this route; transmission starts when the case reaches state from and stops at the earliest of the until states. Give a BranchingProcess several of these for different routes, such as community and funeral transmission of Ebola.

offspring: a distribution of secondary cases, such as NegBin(R, k) (mean R, dispersion k), or a function (rng, ind) -> n returning the number infected by a given case (one count per type in a multi-type model). The function may also take the outbreak state as a third argument (see draw_offspring).
from: the state at which transmission starts, :infection by default. Any other name refers to a state set by a step of the natural history, such as :infectious or :died; transmission starts only once the case reaches it.
until: the states at which transmission stops, for example (:recovered, :died) for community transmission or (:buried,) for a funeral. The earliest one reached ends transmission, and infections that would have happened later do not. Empty by default. Isolation is set with the Isolation intervention, not here.
kernel: the time in days from from to each infection: a distribution, a function of the individual returning a distribution, (ind) -> ..., or NoGenerationTime() to place every infection at the from time. With the default from = :infection and no until, this is the generation time. With until it is the contact interval (the time to a contact that would infect if nothing stopped it), and the generation time follows from which comes first, the contact or the end of transmission, so do not also shorten it by hand. A kernel measured from a later state, such as onset, can put a contact before the case's own infection. No exposure happens then, so that contact is dropped, and a case makes fewer contacts on average than its offspring distribution draws, by the share of the kernel that falls before infection.

Examples

julia

using EpiBranch, Distributions

progression = [
    Transition(:infectious, from = :infection, delay = Gamma(4.0, 2.0)),
    Transition(:died, from = :infectious, delay = Gamma(4.0, 2.0),
        probability = 0.6, terminal = true),
    Transition(:recovered, from = :infectious, delay = Gamma(5.0, 2.0),
        terminal = true),
    Transition(:buried, from = :died, delay = 2.0),
]
community = Infectiousness(NegBin(1.2, 0.5);
    from = :infectious, until = (:recovered, :died), kernel = Exponential(4.0))
funeral = Infectiousness(Poisson(0.5);
    from = :died, until = (:buried,), kernel = Uniform(0.0, 2.0))
model = ModelSpec(BranchingProcess(community, funeral); progression)

Fields

offspring::Any
from::Any
until::Any
kernel::Any

EpiBranch.MultiTypeOffspring Type

julia

struct MultiTypeOffspring{M<:(AbstractMatrix{<:Real}), F, R<:(AbstractVector{<:Real}), A<:(AbstractMatrix{<:Real})}

julia

MultiTypeOffspring(offspring_matrix, dist_fn)

The offspring distribution of a multi-type branching process (for example children and adults, or health-care workers and the community), built from a next-generation matrix M: M[i, j] is the mean number of type-i cases infected by one type-j case.

A type-j case draws how many people it infects in total from dist_fn(R_j), where R_j is the sum of column j, then splits them at random between types in proportion to that column. The mean number infected is the mean of dist_fn(R_j), which equals R_j only if the distribution has mean R_j. A function such as R -> Poisson(θ * R) scales the whole matrix by θ.

Warning

dist_fn must take a mean if the column sums are meant to be the reproduction numbers. Use R -> NegBin(R, k) (mean R, dispersion k). R -> NegativeBinomial(R, 0.3) from Distributions.jl reads R as a number of successes, not a mean: with it, a matrix whose dominant eigenvalue is 1.045 gives a process with reproduction number 2.44. reproduction_number reports the reproduction number the model actually has.

BranchingProcess(M, dist_fn, generation_time) builds one of these; reproduction_number and extinction_probability use the same matrix and distribution.

Fields

offspring_matrix::AbstractMatrix{<:Real}
dist_fn::Any
R_by_type::AbstractVector{<:Real}
alloc_probs::AbstractMatrix{<:Real}

EpiBranch.HomogeneousProcess Type

julia

struct HomogeneousProcess{T<:Real} <: TransmissionModel

julia

HomogeneousProcess(; transmission_rate, population_size,
                   from = nothing, until = (:recovered, :died, :isolated))

A closed population of population_size individuals that mix homogeneously: every infectious person exerts the same force of infection on every susceptible. It is simulated in continuous time as mass action, every pair meeting at the same rate, which reproduces the exact stochastic SIR final-size law. Transmission is the per-infective rate transmission_rate (β, so β/N to each susceptible).

The natural history comes from the progression of a ModelSpec, as for a BranchingProcess. With only a recovery step, Transition(:recovered; rate = γ, terminal = true), this is an SIR model; adding a latent period (a step to :infectious) makes it SEIR. Further steps (onset, hospitalisation, death) appear in the line list.

from: the state at which a case becomes infectious. Left as nothing, it is :infectious when the progression has a latent period, otherwise :infection.
until: the states that end infectiousness, by default (:recovered, :died, :isolated).

Interventions attach through two seams. An intervention that removes a case from transmission, such as Isolation, shortens its infectious window and curtails spread; one whose effect is a per-contact block, such as a leaky isolation or a vaccine's efficacy, is resolved against each contact a pair makes, and a blocked contact leaves the pair to meet again, so blocking a fraction of the contacts thins the force of infection by the same fraction. Per-individual susceptibility and infectiousness reach that thinning directly: they scale the rate at which each pair meets. Four interventions are reported with a warning rather than applied, because mass action settles one case at a time and never produces the batches of new contacts their dosing and tracing act on: MassVaccination's rollout, GroupVaccination, RingVaccination and ContactTracing, which has no contacts to trace under random mixing anyway. Doses a person already holds when the outbreak starts still protect them; only new dosing during the run is dropped. Control expressed as a removal Transition in the progression always applies.

The population is simulated until extinction or max_time, whichever comes first; with max_time, individuals whose infection would fall later stay uninfected. The other simulate termination controls (max_cases, max_generations, stopping_rules other than MaxTime) do not apply, and simulate warns if you set one. This holds for the structure-driven models generally.

A run costs about one draw per contact made inside an infectious window. A window that never closes, as in an SI model with no removal transition, makes each case contact every member, so such a run costs O(N) per case unless max_time bounds it.

Example

julia

using EpiBranch, Distributions
model = ModelSpec(
    HomogeneousProcess(; transmission_rate = 2.0, population_size = 3000);
    progression = [Transition(:recovered; from = :infection, rate = 1.0, terminal = true)])
state = simulate(model; n_initial = 5)

Fields

population_size::Int64
transmission_rate::Real
from::Union{Nothing, Symbol}
until::Tuple

EpiBranch.MixingProcess Type

julia

struct MixingProcess{F} <: TransmissionModel

julia

MixingProcess(; population_size, mixing_by = (), rate,
             from = nothing, until = (:recovered, :died, :isolated))

A closed population of population_size individuals, split into mixing groups by mixing_by, a tuple of attribute keys already set on each individual (:age_band, :ses, :patch; set by the composed attributes, not by this process). A member's group is the tuple of those attribute values; mixing_by = () puts everyone in the single group (), recovering HomogeneousProcess.

rate(infector_group, target_group) is the rate at which one infective in the first group contacts one member of the second. Homogeneous mixing is β/N for every pair; a contact matrix M[b, a] giving the rate at which one infective in group a contacts the n[b] members of group b between them gives M[b, a]/n[b]. It is simulated as mass action on the same continuous-time race as HomogeneousProcess.

The process describes the transmission alone. The natural history is a progression of Transitions attached with a ModelSpec, exactly as for HomogeneousProcess: from is the state the infectious window opens at, derived from the progression when left nothing; until names the removal states that close it.

The population is simulated until extinction or max_time, whichever comes first; the other simulate termination controls do not apply, and simulate warns if you set one.

Every contact is a contact between two named people, so competing risks, per-individual susceptibility and infectiousness, and risks that read the infector all apply as they do on any other route of the race.

Example

julia

using EpiBranch, Distributions

# Two age bands of equal size, tagged on creation; band 1 mixes more.
band = (rng, ind) -> (ind.state[:age_band] = rand(rng, 1:2))
M = [3.0 0.5; 0.5 0.5]
rate = (from, to) -> M[to[1], from[1]] / 1500

model = ModelSpec(
    MixingProcess(; population_size = 3000, mixing_by = (:age_band,), rate);
    progression = [Transition(:recovered; from = :infection, rate = 1.0, terminal = true)],
    attributes = band)
state = simulate(model; n_initial = 5)

Fields

population_size::Int64
mixing_by::Tuple
rate::Any
from::Union{Nothing, Symbol}
until::Tuple

EpiNetwork.NetworkProcess Type

julia

NetworkProcess(adjacency, kernel; from = nothing,
               until = (:recovered, :died, :isolated), external_hazard = 0.0,
               obs_end = Inf)

Transmission over a fixed contact network: who can infect whom is set by the network, and each infectious person infects a neighbour after a random delay, provided they are still infectious by then. Ending someone's infectious period early (recovery, isolation) therefore cuts their onward transmission.

The network is the whole population. adjacency[i] lists the people person i is in contact with (numbered from 1).

kernel is the contact interval, the delay in days from a person becoming infectious to their infecting contact with a neighbour. It can be:

one continuous distribution shared by every contact, such as Exponential(2.0) (mean 2 days);
a function (infector, susceptible) -> Distribution of the two people's numbers, for covariates that change transmission between particular pairs;
a PairKernel, which can also use the infector's infection time (for calendar-time effects) and each person's attributes and history;
a vector of vectors parallel to adjacency, with kernel[i][k] the distribution for person i's k-th listed neighbour.

The natural history is a progression of Transitions attached with a ModelSpec, as for HouseholdProcess. A latent period is Transition(:infectious; from = :infection, delay = …) and the infectious period ends with a terminal removal transition; onset, testing and the rest are further transitions that the line list reports. from is the state from which the contact interval is measured. Left as nothing, it is :infectious when a latent period produces that state and :infection otherwise. until names the states that end the infectious period, so adding a Transition(:isolated, …) cuts onward transmission.

external_hazard adds infections from outside the network: a constant rate per person per day, or a continuous distribution of the time (in days) at which each person would be infected from outside. To let it differ between people, give a vector with one rate or distribution per node, or a function j -> rate or distribution of node j's number. These introductions happen only in the first obs_end days, which must be finite when external_hazard is used.

Example

julia

using EpiNetwork, EpiBranch, Distributions
adjacency = [[2, 5], [1, 3], [2, 4], [3, 5], [4, 1]]   # ring of 5 people
model = ModelSpec(NetworkProcess(adjacency, Exponential(2.0));
    progression = [Transition(:recovered; from = :infection, delay = 6.0, terminal = true)])

EpiNetwork.RoutedNetwork Type

julia

RoutedNetwork(windows; from = nothing, external_hazard = 0.0, obs_end = Inf)

Transmission over a contact network with several routes, such as household and community, each with its own contacts, timing and response to control measures. Isolation can then end a case's community transmission while they go on infecting the people they live with.

Example

Households as fully connected groups, community contact as a sparser network over the same people. Self-isolation ends community transmission and leaves the household route running. In this example household_adjacency and community_adjacency are contact lists, one vector of people per person:

julia

using EpiNetwork, EpiBranch, Distributions

household = RouteWindow(:household; until = (:recovered,),
    kernel = Weibull(1.5, 3.0), reach = household_adjacency)
community = RouteWindow(:community;
    until = (:recovered, EpiBranch.INTERVENTION_REMOVAL),
    kernel = Exponential(12.0), reach = community_adjacency)

model = ModelSpec(RoutedNetwork([household, community]);
    progression = [Transition(:recovered; from = :infection, delay = 7.0,
        terminal = true)],
    interventions = [Isolation(onset_to_isolation_delay = Exponential(2.0), duration = 7.0)],
    attributes = clinical_presentation(incubation_period = LogNormal(1.0, 0.3)))

Because the household route does not list INTERVENTION_REMOVAL, a case who isolates keeps infecting their household to the end of their infectious period, as self-isolation at home does. Each route's kernel describes transmission as if the case were never removed; the realised number of secondary cases falls according to which routes isolation cuts.

Routes

windows is a collection of RouteWindows, each with:

reach: the route's contact list (reach[i] lists the people i can infect on this route), so different routes can connect different pairs of the same people;
kernel: the contact interval along those contacts, in days, in any of the forms NetworkProcess accepts (one distribution, a function (infector, susceptible) -> Distribution, a PairKernel, or a per-contact vector parallel to reach). A PairKernel that reads people's records lists them in its watches, and only the routes reading a record that changes have their contacts redrawn;
until: the states that end transmission on this route, so one route can be cut and another left alone. Include EpiBranch.INTERVENTION_REMOVAL for a route that isolation or quarantine should end;
from: when the route's infectiousness starts. The default, nothing, takes the model's from or, if that is also nothing, the start the progression implies (:infectious when a latent period produces it). Name :infection for a route infectious from the moment of infection, or a later state for one that starts later, such as a funeral route from :died;
contacts_from: when the route's contacts become the case's contacts for tracing. Leave it at :infection for standing relationships; give a funeral route contacts_from = :died so its contacts are traced only once the funeral has happened;
traceable: the probability that a case can name a contact on this route. Contact tracing reaches only named contacts. Keep the default 1.0 for people a case can always name, such as their household, and lower it for casual contact.

All routes cover the same people, so every reach must have the same length. external_hazard and obs_end are as for NetworkProcess.

Whether a case names a given contact is decided once, when the case's contacts are passed to tracing. A contact reachable on several routes is named with the highest of those routes' traceable probabilities, and can be traced from the earliest contacts_from time among the routes that name them; for example, with traceable = 0.2 on the community route, a case can name one community contact in five but everyone they live with, including household members they also meet in the community. The intervention then traces a named contact with its own probability, so the two multiply.

A case infected on one of these routes records the route's name in its :infection_route, and a community introduction records :external. Both appear in linelist, so cases can be broken down by setting.

EpiHouseholds.HouseholdProcess Type

julia

HouseholdProcess(sizes, kernel; from = nothing, until = (:recovered, :died, :isolated),
                 external_hazard = 0.0, obs_end = Inf)

Transmission within households: each infectious person infects their household members after a random delay, provided they are still infectious by then, and infections from outside the household arrive at a community rate. sizes gives the size of each household, so there are sum(sizes) people in length(sizes) households.

kernel is the within-household contact interval, the delay in days from a person becoming infectious to their infecting contact with a household member. It can be any continuous distribution on the positive reals, a function (infector, susceptible) -> Distribution of the two people's numbers for covariates, or a PairKernel, which can also use the infector's infection time and each person's attributes and history.

The natural history is a progression of Transitions attached with a ModelSpec. A latent period is Transition(:infectious; from = :infection, delay = …) and the infectious period ends with a terminal removal transition; onset, testing and the rest are further transitions that the line list reports. from is the state from which the contact interval is measured. Left as nothing, it is :infectious when a latent period produces that state and :infection otherwise. until names the states that end the infectious period.

external_hazard is infection from the community: a constant rate per person per day, or a distribution of the time (in days) at which each person would be infected from outside, for a rate that changes over time. To let it differ between people, households or trial arms, give a vector with one rate or distribution per person, or a function j -> rate or distribution of person j's number, as the kernel is a function of two people's numbers. A person with a zero rate is never infected from outside. These introductions happen only in the first obs_end days, which must be finite when external_hazard is used. Without an external hazard, each household starts with one index case at day 0, unless initial_cases is given, in which case only those people are infected at day 0.

A progression listing Transition(:susceptible_again; from = :recovered, delay = …) lets a household member be infected again once their immunity has waned: a still-infectious household-mate, or a later community introduction, can reach them from susceptible_again_time onward, and their earlier infection is kept on Individual.episodes rather than overwritten (see close_episode!). Without such a step every member is infected at most once, as before.

What isolation means here

The only transmission this process represents is within a household; community infection enters as external_hazard introductions, and there is no contact between households. An intervention that ends a case's infectious period here therefore stops them infecting their own household, which physically means removing them from the household: hospitalisation, or transfer to an isolation facility. It does not model self-isolation at home, which would leave household transmission going and, in most settings, increase it. Nor is there a community route for a case to be isolated from while they stay infectious to the people they live with, which is what self-isolation does. Representing that needs transmission split into a household route and a community route, so a control measure can cut one and leave the other; see the route windows in the design notes. Read isolation and quarantine on this process as removal from the household, and choose the parameters accordingly.

Interventions on this process:

Isolation removes a case from their isolation time for its duration, cutting their secondary cases; with duration = Inf this ends their infectious period. Leaky isolation multiplies the rate at which each household contact is infected by post_isolation_transmission while the case is isolated, which lowers their chance of infection proportionally less than it lowers the rate.
ContactTracing treats a case's household members as their contacts; quarantining a traced contact stops that contact's transmission for the quarantine's duration. Tracing can only reach household members infected after the case's own course of infection is known, which in a fast-spreading household means the ones infected later.
A leaky vaccine's efficacy, and per-person susceptibility and infectiousness, reduce the rate at which each household contact is infected. An all-or-nothing vaccine instead fully protects its share of the vaccinated from their immunity date.
RingVaccination and GroupVaccination work, including with Scheduled and CapacityConstrained, and a ring dose's post_exposure_efficacy can stop a household member's own infection before onset. A capacity limit is one budget shared by all households, renewed every period days of the outbreak. RingVaccination's eligibility_window admits a household member's dose before their own infection time is known, and withdraws it once that infection settles and the window turns out not to cover it, freeing any capacity it used.
MassVaccination doses every household's members before its own race starts, needing no household membership at all.
A removal Transition in the progression always applies.

Limits:

Existing protection can instead be given through transmission_traits, a PairKernel, or an intervention with its own competing_risk.

Example

julia

using EpiHouseholds, EpiBranch, Distributions
model = ModelSpec(HouseholdProcess([3, 4, 2], Weibull(1.5, 3.0));
    progression = [Transition(:recovered; from = :infection, delay = 6.0, terminal = true)])

EpiHouseholds.household_sizes Function

julia

household_sizes(model) -> Vector{Int}

The size of each household in model.

EpiHouseholds.household_members Function

julia

household_members(model) -> Vector{Vector{Int}}

The individual ids belonging to each household in model, so household_members(model)[h] lists the members of household h.

EpiHouseholds.external_hazard Function

julia

external_hazard(model)

The community force of infection in model: a constant rate, a distribution of introduction times, or a per-person vector or function of those, as given to HouseholdProcess. Zero means no community introductions.

EpiHouseholds.observation_end Function

julia

observation_end(model) -> Float64

The end, in days, of the window in which model's community introductions can happen. Inf when model has no external_hazard.

EpiBranch.ModelSpec Type

julia

ModelSpec(process; progression, interventions, attributes, observation, recorder)

Combine a transmission model (who infects whom, and when) with the rest of an outbreak scenario. Pass the result to simulate, or to loglikelihood(data, model) to compare it with observed data (the data are not stored in the model).

progression: the natural history of a case, as a vector of Transitions and other steps such as Recovery and Death: latent period, symptom onset, hospitalisation, outcome.
interventions: control measures, such as Isolation, ContactTracing and RingVaccination.
attributes: population characteristics drawn for each person, such as clinical_presentation (incubation period, asymptomatic fraction) and demographics (age, sex). Give one, or several in a vector.
observation: which outbreaks are recorded, such as MinimumSize (only chains of at least a given size). Report individual cases with a Reporting step in progression.
recorder: rarely needed; see ContactRecorder.

Every keyword left out keeps the value the transmission model already has, so ModelSpec(process) behaves exactly like process.

Examples

julia

using EpiBranch, Distributions

model = ModelSpec(
    BranchingProcess(NegBin(2.5, 0.16), Gamma(2.5, 2.0));
    attributes = clinical_presentation(incubation_period = LogNormal(1.6, 0.5)),
    interventions = [Isolation(onset_to_isolation_delay = Exponential(2.0),
        duration = 14.0)],
)
state = simulate(model; max_cases = 500)

EpiBranch.single_type_offspring Function

julia

single_type_offspring(model::TransmissionModel) -> Any

julia

single_type_offspring(model::TransmissionModel)

The offspring distribution of a single-type model: usually a distribution such as NegBin(R, k), or anything else chain_size_distribution accepts, such as ClusterMixed. Raises an error for a multi-type model.

The closed-form results read the offspring distribution through this function, so a new transmission model that defines it can use all of them: reproduction_number, extinction_probability, epidemic_probability, probability_contain, proportion_transmission, proportion_cases_individual, proportion_cases_offspring, proportion_cluster_size and chain_size_distribution.

EpiBranch.has_single_type_law Function

julia

has_single_type_law(model::TransmissionModel) -> Bool

julia

has_single_type_law(model::TransmissionModel)

Whether single_type_offspring(model) returns an offspring specification with a closed-form chain_size_distribution method, rather than single_type_offspring throwing or the offspring having no such method. chain_size_distribution and loglikelihood(::ChainSizes, model) check this before taking the analytical fast path, rather than each catching the errors single_type_offspring and chain_size_distribution throw for a multi-type/multi-window model or an offspring with no closed form.

The default calls single_type_offspring and conservatively reports false (fall back to simulation) for anything it throws an ArgumentError or MethodError over, and likewise checks whether chain_size_distribution has a method for the result, so a custom TransmissionModel need not specialise this to get correct, if not always fastest, behaviour.

EpiBranch.offspring_law_unchanged Function

julia

offspring_law_unchanged(_) -> Any

julia

offspring_law_unchanged(x) -> Bool

Whether x leaves a model's single-type offspring law exactly as the offspring distribution draws it, with nothing reshaping how many of the drawn contacts end up as cases. loglikelihood(::ChainSizes, model), chain_size_distribution and the closed-form incidence likelihood take their analytical path, scoring the bare offspring distribution, only once the process, every entry of attributes and every progression transition each answer true here; any one of them answering false sends the evaluation to simulation instead.

For an attribute or a progression transition, the default is true: most of what either can hold — the natural history, population characteristics such as age or a group label, reporting — has no bearing on the offspring law at all. transmission_traits overrides it to false once susceptibility or infectiousness departs from its neutral value of 1, since either one thins how many of the drawn contacts actually become cases. A Tuple or AbstractVector of attributes answers true only once every entry does.

A bare Function or closure answers the conservative false instead, whatever it does: an attributes entry can set ind.susceptibility or ind.infectiousness directly rather than going through transmission_traits (its docstring recommends exactly this for rules it doesn't cover), and this trait has no way to tell that apart from one that leaves them alone. demographics and groups return a dedicated callable type rather than a closure so they keep the true default; give a custom attributes function one too if it should.

A custom AbstractClinicalTransition answers the same conservative false, for the same reason: initialise_individual! and resolve_individual! receive the case directly and can set its susceptibility or infectiousness just as a bare attributes function can. The five built-in transitions (Transition, Reporting, Hospitalisation, Recovery, Death) touch neither field; each overrides it back to true, and one of your own that doesn't touch them either should do the same.

For a TransmissionModel the default is the conservative false instead, as for EpiBranch.infection_likelihood_compatible: a new model opts in once it has checked that nothing of its own — a finite population, a window that can close early, or anything else — reshapes the offspring law away from the bare distribution. BranchingProcess opts in except for a finite population_size (fewer susceptible people than an unlimited offspring law assumes) or an infectiousness window with a non-empty until (transmission can stop before every drawn contact happens).

Population characteristics

Functions that give each case its characteristics when it is created: symptom onset and asymptomatic status, age and sex, susceptibility and infectiousness, group membership, and willingness to be vaccinated. Pass them to ModelSpec as attributes.

EpiBranch.clinical_presentation Function

julia

clinical_presentation(
;
    incubation_period,
    prob_asymptomatic
)

julia

clinical_presentation(; incubation_period, prob_asymptomatic = 0.0)

Give each case a symptom onset time and decide whether it is asymptomatic. Pass the result as attributes to ModelSpec.

incubation_period: distribution of the time from infection to symptom onset, in days. Each symptomatic case's onset time is its infection time plus a draw from it.
prob_asymptomatic: probability that a case never develops symptoms (default 0). Asymptomatic cases have no onset time (NaN). Give a number, a distribution (one probability drawn per case) or a function of the random number generator and the individual, (rng, ind) -> probability, for example to make it depend on age.

Isolation needs these onset times, and linelist reports them as date_onset.

Examples

Symptomatic-only with a log-normal incubation period:

julia

attributes = clinical_presentation(incubation_period = LogNormal(1.6, 0.5))

With 30% asymptomatic:

julia

attributes = clinical_presentation(
    incubation_period = LogNormal(1.6, 0.5),
    prob_asymptomatic = 0.3,
)

Per-individual asymptomatic probability drawn from a Beta:

julia

attributes = clinical_presentation(
    incubation_period = LogNormal(1.6, 0.5),
    prob_asymptomatic = Beta(2, 8),
)

Age-conditional (children much more likely to be asymptomatic; list after demographics so :age is set first):

julia

attributes = [
    demographics(age_distribution = Uniform(0, 90)),
    clinical_presentation(
        incubation_period = LogNormal(1.6, 0.5),
        prob_asymptomatic = (rng, ind) -> ind.state[:age] < 18 ? 0.6 : 0.2,
    ),
]

Clinical progression (natural history)

The steps of a case's natural history after infection, each with a probability and a delay in days: becoming infectious after a latent period, symptom onset, reporting, hospitalisation, recovery or death. Pass them to ModelSpec as progression.

EpiBranch.AbstractClinicalTransition Type

Parent type of the steps in a case's natural history: Transition, Reporting, Hospitalisation, Recovery and Death. List several in a model's progression; each case passes through them in order as it is settled, after its population characteristics are set and before the interventions act on it.

Steps that end a case, such as death or recovery, are called terminal. When a case could reach several terminal steps, the earliest one becomes its outcome (outcome and date_outcome in the line list), so death and recovery compete and whichever comes first wins.

Interventions are the policy applied to a case (AbstractIntervention); progression steps are the biology that happens to it. Writing a new kind of step, such as loss to follow-up, needs the methods below.

Once a case has reached its outcome, an event one of the other steps recorded for a later time is reset to "did not happen", since nothing more happens to a case that has died or recovered.

For extension authors

A new step defines initialise_individual! (set default state on a new case) and resolve_individual! (draw whether and when the step happens, and write the result to individual.state under keys it owns). A terminal step also returns true from is_terminal, implements terminal_event, and should implement terminal_target so the check that every infectious window closes can see the new end state.

While a step resolves, it sees the same simulation state the interventions do: state.max_infection_time at the case's infection time, state.current_generation at its generation, and state.cumulative_cases counting the cases settled up to and including it.

EpiBranch.Transition Type

julia

struct Transition{D, P, F} <: AbstractClinicalTransition

julia

Transition(state; from = :infection, delay = …, rate = …, probability = 1.0, terminal = false)

One step in a case's natural history, such as becoming infectious after a latent period or recovering after an infectious period: the case reaches state a delay (in days) after it reached from, with probability probability. If the step does not happen, because it was not drawn or the case never reached from, the case simply never reaches state.

from: the earlier state the delay is measured from. :infection (the default) is the time of infection; from = :onset measures from symptom onset, and in general from = :s measures from the time recorded for state s (from = :onset_time works too). A function of the individual returning a time is also accepted. If the case never reached from (including asymptomatic cases, whose onset time is NaN), the step is skipped. If no attribute or other step records that state at all, the simulation stops with an error when it starts.
delay or rate (give exactly one): delay is a fixed number of days, a distribution, or a function of the random number generator and the individual, (rng, ind) -> ..., drawn for each case. rate = r is the compartmental-model alternative, an exponentially distributed delay with mean 1 / r days (delay = Exponential(1 / r)).
probability: the chance the step happens, a number or a function (rng, ind) -> ....
terminal = true marks a step that ends the case. When a case can reach several terminal steps, the earliest becomes its outcome (:outcome, :outcome_time).

Each case records state as reached (true/false) and the time it was reached (Inf if never), under the name Symbol(state, :_time), so Transition(:infectious, ...) gives :infectious_time. Reporting, Hospitalisation, Recovery and Death are ready-made versions of this step with fixed names.

A non-terminal transition (terminal = false, the default) that resolves after the case's outcome is reset to "did not occur": see censor_after_outcome!.

Examples

julia

# latent period: infection → onset of infectiousness
Transition(:infectious, from = :infection, delay = LogNormal(1.0, 0.4))

# infectious period as a recovery rate (exponential, mean 1/γ)
Transition(:recovered, from = :infectious, rate = 1 / 6, terminal = true)

# severity branch, then death from the severe state
Transition(:severe, from = :onset,  delay = Gamma(2, 2), probability = 0.3)
Transition(:died,   from = :severe, delay = Gamma(2, 3), probability = 0.6, terminal = true)

Fields

state::Symbol
time_key::Symbol
delay::Any
probability::Any
from::Any
terminal::Bool

EpiBranch.Reporting Type

julia

struct Reporting{D, P, F} <: AbstractClinicalTransition

julia

Reporting(; delay, probability = 1.0, from = :onset_time)

Case reporting: each case is reported with probability probability, delay days after symptom onset (or after from, if given). Each case records :reported (true/false) and :reporting_time (Inf if never reported).

delay: a distribution, a fixed number of days, or a function of the random number generator and the individual, (rng, ind) -> ....
probability: the chance a case is reported, a number or a function (rng, ind) -> ... (for example a different detection probability per risk group).
from: the time the delay is measured from. Symptom onset by default; another recorded time such as :test_time (:test works too) or :admission_time; or a function of the individual, such as ind -> ind.infection_time.

Cases that never reached from are not reported. With the default from, this includes asymptomatic cases from clinical_presentation, whose onset time is NaN. The time measured from must be one that the model's attributes set or that a step in its progression records; the simulation checks this when it starts, and stops with an error if neither does. A report drawn for a time after the case has died or recovered is reset to "not reported".

simulate, progression_loglik, the chain-size likelihood and end_of_outbreak_probability all read this step's reporting probability. It replaces the deprecated PerCaseObservation.

Examples

julia

# report 80% of symptomatic cases, on average 3 days after onset
Reporting(delay = Gamma(2.0, 1.5), probability = 0.8)

Fields

delay::Any
probability::Any
from::Any

EpiBranch.Hospitalisation Type

julia

struct Hospitalisation{D, P, F} <: AbstractClinicalTransition

julia

Hospitalisation(; delay, probability = 0.2, from = :onset_time)

Hospital admission: each case is admitted with probability probability (default 0.2), delay days after symptom onset (or after from). Each case records :admitted (true/false) and :admission_time (Inf if never admitted). delay, probability and from take the same forms as in Reporting; cases that never reached from are not admitted.

Examples

julia

# 20% admitted, a median of about 7 days after onset
Hospitalisation(delay = LogNormal(2.0, 0.5), probability = 0.2)

To admit only cases that meet some other condition (reported, traced, vaccinated, ...), return 0.0 from a probability function when the condition does not hold. Here only reported cases can be admitted; list Reporting before Hospitalisation in the progression so reporting is drawn first:

julia

Hospitalisation(
    delay = LogNormal(2.0, 0.5),
    probability = (rng, ind) -> get(ind.state, :reported, false) ? 0.2 : 0.0
)

An admission drawn for a time after the case has died or recovered is reset to "not admitted". To have admission change the outcome instead of being dropped, give the terminal step from = :admission_time.

Fields

delay::Any
probability::Any
from::Any

EpiBranch.Death Type

julia

struct Death{D, P, F} <: AbstractClinicalTransition

julia

Death(; delay, probability, from = :onset_time)

Death, a step that ends the case: a case dies with probability probability, delay days after symptom onset (or after from), unless an earlier terminal step ends the case first.

The realised case fatality ratio can be lower than probability

Death competes with other terminal steps, and the earliest wins. With a competing Recovery, a case drawn to die still recovers if its recovery time comes first, so fewer than probability of cases die (with equal delays, about half as many). For an exact case fatality ratio, replace Recovery with a terminal Transition and split the two outcomes with exclusive_probabilities, so each case gets exactly one:

julia

death_p, recovered_p = exclusive_probabilities([0.64, 0.36])
progression = [
    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),
]

Giving death and recovery independent probabilities CFR and 1 - CFR does not work: about CFR * (1 - CFR) of cases would qualify for both and another CFR * (1 - CFR) for neither, leaving them with no outcome. The alternative is to make the delay to death reliably shorter than any competing delay.

probability has no default, because the case fatality ratio is too disease-specific for one; pass it explicitly, even if it is 0.0. It is a number or a function of the random number generator and the individual, (rng, ind) -> ..., for example to make it depend on age:

julia

Death(delay = LogNormal(2.5, 0.4),
      probability = (rng, ind) -> ind.state[:age] >= 80 ? 0.3 : 0.02)

delay is a distribution, a fixed number of days or such a function. from takes the same forms as in Reporting.

A vaccine that lowers mortality (a RingVaccination or MassVaccination with severity_efficacy) acts through the same function, using severity_efficacy and immunity_time:

julia

Death(delay = LogNormal(2.5, 0.4),
      probability = (rng, ind) ->
          immunity_time(ind) <= onset_time(ind) ?
              0.7 * (1 - severity_efficacy(ind)) : 0.7)

Comparing immunity_time(ind) with onset_time(ind) means a dose whose protection had not developed by symptom onset gives none; testing is_vaccinated(ind) alone would count it as protective. Combine it with an age-dependent case fatality ratio by multiplying that ratio by 1 - severity_efficacy(ind) once the case is immune.

Fields

delay::Any
probability::Any
from::Any

EpiBranch.Recovery Type

julia

struct Recovery{D, F} <: AbstractClinicalTransition

julia

Recovery(; delay, from = :onset_time)

Recovery, a step that ends the case: every case that reaches from (symptom onset by default) recovers delay days later, unless an earlier terminal step such as Death ends the case first. The earliest terminal step becomes the case's :outcome and :outcome_time.

delay is a distribution, a fixed number of days or a function of the random number generator and the individual, (rng, ind) -> ... (for example a recovery delay that depends on age). from takes the same forms as in Reporting; cases that never reached it get no recovery time. Recovery has no probability: it always happens once from is reached.

Examples

julia

# recover a mean of 10 days after onset, unless death comes first
Recovery(delay = Gamma(4.0, 2.5))

Fields

delay::Any
from::Any

EpiBranch.exclusive_probabilities Function

julia

exclusive_probabilities(ps::AbstractVector{<:Real}) -> Any

julia

exclusive_probabilities(ps::AbstractVector{<:Real}) -> Vector

Split cases between mutually exclusive outcomes in fixed proportions, for example exactly 64% die and 36% recover. Returns one value per outcome, to pass as each outcome's probability.

Giving each outcome a plain probability instead draws them independently: with death at p and recovery at 1 - p, about p(1 - p) of cases would qualify for both (and take whichever comes first) and another p(1 - p) for neither, so the realised case fatality ratio drifts from p. With exclusive_probabilities, each case falls into exactly one outcome.

ps must be non-negative and sum to at most 1. Any shortfall below 1 is the probability that none of the outcomes happens; pair it with an unconditional terminal step, or expect some cases to have no outcome.

Pass each returned value straight to probability. Wrapping it in a function of your own makes it an ordinary 0-or-1 probability, which loses the shared split in the likelihood; fold any per-case modifier into ps instead. Keep every outcome of the group in the progression, since progression_loglik evaluates the group as one event.

Examples

julia

death_p, recovered_p = exclusive_probabilities([0.64, 0.36])
progression = [
    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),
]

EpiBranch.is_terminal Function

julia

is_terminal(_::AbstractClinicalTransition) -> Bool

Whether this step ends the case (as death or recovery do). false by default.

EpiBranch.terminal_event Function

julia

terminal_event(
    _::AbstractClinicalTransition,
    individual
) -> Union{Nothing, Tuple{Any, Symbol}}

julia

terminal_event(transition, individual) -> Union{Nothing, Tuple{Float64, Symbol}}

For a terminal step, the time (days since the start of the simulation) and end state at which it would end this case, such as (11.3, :died), or nothing if it does not happen to this case. The earliest across all terminal steps becomes the case's outcome, recorded as :outcome and :outcome_time. Only called for terminal steps, after every step has been drawn.

EpiBranch.terminal_target Function

julia

terminal_target(
    _::AbstractClinicalTransition
) -> Union{Nothing, Symbol}

julia

terminal_target(transition::AbstractClinicalTransition) -> Union{Nothing, Symbol}

The end state a terminal step records, such as :died, known before any case is simulated. Model checks use it to warn when an infectious window's until misses an end state a case can reach. A terminal step that leaves this at the default nothing is skipped by that warning. Non-terminal steps never need it.

EpiBranch.progression_loglik Function

julia

progression_loglik(
    spec::ModelSpec,
    individuals::AbstractVector{<:Individual}
) -> Any

julia

progression_loglik(spec::ModelSpec, individuals) -> Float64

The log-likelihood of the cases' clinical timelines under the steps of the model's progression (reporting, hospitalisation, death or recovery, each measured from the time it starts from). Use it to estimate the delays in those steps, such as a reporting delay, from a line list with event times. An incubation period drawn by clinical_presentation is not included unless onset is itself a step of the progression.

individuals is a vector of Individual or a SimulationState from simulate, or a hand-built equivalent holding the times each progression step records. People who were exposed but never infected add nothing. The result is the sum over cases and steps of transition_loglik.

It covers natural history only. Who infected whom and when is the transmission part, evaluated by pairwise_surv_loglik; add the two for the full log-likelihood of an outbreak with known infection times. Neither includes an observation model.

Examples

julia

using EpiBranch, Distributions, StableRNGs

progression = [
    Reporting(delay = LogNormal(1.0, 0.3)),
    Recovery(delay = LogNormal(2.0, 0.4)),
]
spec = ModelSpec(
    BranchingProcess(Poisson(2.0), Exponential(5.0));
    progression = progression,
    attributes = clinical_presentation(incubation_period = LogNormal(1.5, 0.5))
)
state = simulate(spec; max_cases = 100, rng = StableRNG(1))
progression_loglik(spec, state)

Interventions

Control measures applied during an outbreak, passed to ModelSpec as interventions. Isolation applies to cases; quarantine applies to traced contacts who are not (yet) known cases.

EpiBranch.AbstractIntervention Type

julia

abstract type AbstractIntervention

Parent type of all interventions: Isolation, ContactTracing, RingVaccination, MassVaccination and GroupVaccination. Wrap one in Scheduled to start or stop it at a given time or case count, or in CapacityConstrained to limit how many people it can reach, per period or in total. Pass a vector of them as interventions to ModelSpec.

Measures that limit the number of secondary cases directly (a cap per case, a limit on gathering size) are written as an offspring distribution that depends on the simulation state, passed to BranchingProcess.

For extension authors: a new intervention is a subtype that defines one or more of initialise_individual!, resolve_individual!, trace_contacts!, intervention_actions and competing_risk. To work with Scheduled it also defines intervention_time and reset!, whose defaults (-Inf and doing nothing) suit an intervention that always acts immediately. See the Extending EpiBranch guide.

Fields

Isolation

EpiBranch.Isolation Type

julia

struct Isolation{E<:IsolationEligibility, D, S, U} <: AbstractIntervention

julia

Isolation(; onset_to_isolation_delay, duration,
          eligibility = SymptomaticOnly(), test_sensitivity = 1.0,
          post_isolation_transmission = 0.0)

Isolate cases a delay after symptom onset, so they stop (or reduce) onward transmission. Isolation applies to cases; traced contacts who are not yet known cases are quarantined through ContactTracing.

Arguments

onset_to_isolation_delay: days from symptom onset to isolation. A number, a distribution such as Exponential(2.0) (mean 2 days), or a function of the random number generator and the individual, (rng, ind) -> ..., drawn for each case. The function can read information another intervention has recorded on the case, for example to shorten the delay once a case in the same household has been detected.
duration (required): days a case stays isolated before release. A number, a distribution or a function (rng, ind) -> .... Inf isolates for good; with zero the case is still recorded as isolated, so tracing can start from it, but no transmission is prevented. There is no default because indefinite isolation is a modelling choice to make explicitly.
eligibility: who can be isolated, SymptomaticOnly (default) or AllCases (for example with mass testing).
test_sensitivity: probability that an eligible case is detected and isolated (default 1). A number, a distribution or a function (rng, ind) -> ..., drawn once per case.
post_isolation_transmission: fraction of normal transmission that continues while isolated, between 0 and 1 (default 0, perfect isolation; 0.25 would be leaky isolation).

After a finite duration the case transmits again at its normal rate from the release time (isolation_release_time) until its infectious period ends. A finite duration therefore gives the same result as Inf only when the case stops being infectious before it is released.

An isolation at or after the case's own outcome (recovery, death or any other terminal Transition) still stops transmission, so a route that opens at the outcome, such as a funeral, is cut. It does not count as a detection: is_isolated stays false, so OnIsolation tracing and the line list do not count it. records_isolation makes that choice and can be changed through a custom eligibility rule.

Examples

julia

iso = Isolation(onset_to_isolation_delay = Exponential(2.0), duration = Inf)
model = BranchingProcess(NegBin(2.5, 0.16), LogNormal(1.6, 0.5))
spec = ModelSpec(model; interventions = [iso], attributes = clinical_presentation(incubation_period = LogNormal(1.5, 0.5)))
state = simulate(spec; max_cases = 500)

# Leaky isolation lasting 14 days, open to all cases
Isolation(
    onset_to_isolation_delay = 1.0, duration = 14.0,
    eligibility = AllCases(), post_isolation_transmission = 0.25,
)

For extension authors: the case records :isolated, :isolation_time, :isolation_release_time and :test_positive in ind.state.

Fields

eligibility::IsolationEligibility
onset_to_isolation_delay::Any
test_sensitivity::Any
post_isolation_transmission::Float64
duration::Any

EpiBranch.IsolationEligibility Type

julia

abstract type IsolationEligibility

julia

IsolationEligibility

Rule for which cases can be isolated: SymptomaticOnly (the default) or AllCases. Whether an eligible case is actually isolated also depends on the intervention's test_sensitivity.

To write a new rule, define a subtype and a method of is_eligible_for_isolation.

Fields

EpiBranch.SymptomaticOnly Type

julia

struct SymptomaticOnly <: IsolationEligibility

Only symptomatic cases can be isolated (the default for Isolation).

Fields

EpiBranch.AllCases Type

julia

struct AllCases <: IsolationEligibility

Every case is eligible, including asymptomatic ones (mass-testing scenarios).

Fields

Contact tracing

Whose contacts are traced, what share of contacts is found, how long tracing takes and what happens to a traced contact.

EpiBranch.ContactTracing Type

julia

struct ContactTracing{E<:TraceEligibility, F<:TraceRate, D<:TraceDelay, A<:TraceAction} <: AbstractIntervention

julia

ContactTracing(eligibility, probability, isolation_to_trace_delay, action;
               depth = 1)
ContactTracing(; probability, isolation_to_trace_delay, action,
               eligibility = SymptomaticParent(), depth = 1)

Trace the contacts of cases, and quarantine (or just record) the contacts found. Isolation applies to cases (Isolation); quarantine applies to traced contacts who are not yet known cases.

Arguments

eligibility: which cases have their contacts traced and from when, for example OnSymptomOnset (clinical suspicion, from onset), OnLabConfirmation, OnIsolation, TraceEveryone or TraceNobody. Combine with &, | and !. The keyword form defaults to SymptomaticParent (symptomatic and isolated).
probability: share of an eligible case's contacts that tracing finds.
isolation_to_trace_delay: distribution of days from when tracing of the case starts (its isolation by default, onset for OnSymptomOnset) to each contact being reached.
action (required): what happens to a contact who is found, Quarantine for a given duration or FlagOnly to record the trace without quarantining. There is no default, because policies differ in whether and for how long traced contacts are quarantined.
depth: how far tracing reaches: 1 (default) traces contacts, 2 also contacts of contacts, and so on. Must be at least 1.

The keyword quarantine_on_trace is deprecated: quarantine_on_trace = true means action = Quarantine(duration = Inf) and false means action = FlagOnly(). It still works but gives a deprecation warning.

Examples

julia

# Trace on symptoms, without waiting for confirmation; quarantine for 7 days
ContactTracing(OnSymptomOnset(), 0.8, Exponential(1.0), Quarantine(duration = 7.0))

# Wait for lab confirmation; quarantine for 14 days
ContactTracing(OnLabConfirmation(), 0.6, Exponential(2.0), Quarantine(duration = 14.0))

# Trace suspected or confirmed cases
ContactTracing(
    OnSymptomOnset() | OnLabConfirmation(), 0.7, Exponential(1.5), Quarantine(duration = 7.0)
)

# Keyword form with the default eligibility (symptomatic and isolated)
ContactTracing(
    probability = 0.7, isolation_to_trace_delay = Exponential(1.0),
    action = Quarantine(duration = 7.0)
)

Ring depth

With depth = 2 the traced contacts of a case have their own contacts traced in turn: the contacts of contacts that level-2 ring vaccination targets. Each infected, eligible case starts a fresh ring of radius depth. Uninfected ring members are followed for one more generation so the ring can extend past them, but they do not infect anyone. A traced contact who later becomes an eligible case starts its own ring like any other case; add !PreviouslyTraced() to the eligibility to interview each case only once.

Pair with RingVaccination to vaccinate the whole ring:

julia

[ContactTracing(OnSymptomOnset(), 0.8, Exponential(1.0), Quarantine(duration = 7.0); depth = 2),
 RingVaccination(efficacy = 0.9)]

Requirements and output

Needs symptom onset and the asymptomatic flag from clinical_presentation, and, depending on the eligibility, isolation or test results from Isolation. Each contact records whether it was traced and quarantined, and trace_time, the day it was reached, from which interventions acting on traced contacts are timed.

Custom eligibility

A new eligibility rule is a subtype of TraceEligibility with a method of is_eligible. Case characteristics set by clinical_presentation or demographics are in infector.state:

julia

struct SymptomaticOver65 <: TraceEligibility end

function EpiBranch.is_eligible(::SymptomaticOver65, infector, contact, state)
    return !is_asymptomatic(infector) && get(infector.state, :age, 0) >= 65
end

Fields

eligibility::TraceEligibility
trace_rate::TraceRate
isolation_to_trace_delay::TraceDelay
action::TraceAction
depth::Int64

EpiBranch.SourceTracing Type

julia

struct SourceTracing{E<:TraceEligibility, F<:TraceRate, D<:TraceDelay, A<:TraceAction} <: AbstractIntervention

julia

SourceTracing(eligibility, probability, isolation_to_trace_delay, action)
SourceTracing(; probability, isolation_to_trace_delay, action,
              eligibility = SymptomaticParent())

Trace a case back to its own infector, once the case itself meets eligibility, and act on the infector found this way (quarantine or just record it). Where ContactTracing reaches forward, to the people a case infects, SourceTracing reaches backward, to the person who infected it, through the same TraceEligibility, TraceRate, TraceDelay and TraceAction protocol: eligibility and its timing are read off the case found eligible, as they are for ContactTracing, and the action is applied to its infector instead of to a forward contact.

Arguments

Same as ContactTracing: eligibility (which cases are traced back, default SymptomaticParent), probability (share of eligible cases whose infector is found), isolation_to_trace_delay and action (required).

Engine differences

By the time a case is itself confirmed and traced back, a branching process has usually already generated its infector's secondary cases in one step, together rather than over time, so there is nothing left for SourceTracing to shorten there; it still records the link. On a continuous-time model — the household and network processes, and the mass-action pools, where ContactTracing has no pairwise structure to act along at all — an infector's contacts already proposed but not yet resolved are checked against every intervention's competing_risk as each is reached, so quarantining the infector here blocks its remaining onward transmission without undoing any infection already settled.

Examples

julia

# Quarantine a case's own infector once the case is isolated
SourceTracing(
    probability = 0.7, isolation_to_trace_delay = Exponential(1.0),
    action = Quarantine(duration = 14.0)
)

Requirements and output

Same requirements as ContactTracing for the same eligibility. The infector found this way is marked traced (and quarantined, under Quarantine), with its trace_time and, under the key :source_traced_by, the id of the case whose tracing reached it — the reverse of ContactTracing's :traced_by, which always names an earlier- infected node. :source_traced_by is kept separate because it names a later one instead.

Fields

eligibility::TraceEligibility
trace_rate::TraceRate
isolation_to_trace_delay::TraceDelay
action::TraceAction

EpiBranch.TraceEligibility Type

julia

abstract type TraceEligibility

julia

TraceEligibility

Rule for which cases have their contacts traced, and from when. Built in: OnSymptomOnset, OnLabConfirmation, OnIsolation, SymptomaticParent (the default), PreviouslyTraced, TraceEveryone and TraceNobody. Combine them with & (and), | (or) and ! (not), for example OnSymptomOnset() | OnLabConfirmation() to trace suspected or confirmed cases.

To write a new rule, define a subtype and a method of is_eligible.

Fields

EpiBranch.OnSymptomOnset Type

julia

struct OnSymptomOnset <: TraceEligibility

Trace contacts of symptomatic cases (clinical suspicion), starting from symptom onset, so tracing can begin before lab confirmation or without it. See trigger_time.

Fields

EpiBranch.OnLabConfirmation Type

julia

struct OnLabConfirmation <: TraceEligibility

Trace contacts of cases that tested positive (lab confirmation), starting from the case's isolation.

Fields

EpiBranch.OnIsolation Type

julia

struct OnIsolation <: TraceEligibility

Trace contacts of cases that have been isolated, starting from isolation.

Fields

EpiBranch.TraceEveryone Type

julia

struct TraceEveryone <: TraceEligibility

Trace the contacts of every case, whatever its status, starting from the case's isolation.

Fields

EpiBranch.TraceNobody Type

julia

struct TraceNobody <: TraceEligibility

Trace no contacts.

Fields

EpiBranch.PreviouslyTraced Type

julia

struct PreviouslyTraced <: TraceEligibility

Trace contacts only of cases that were themselves traced as someone else's contact. Negate it to stop a traced contact who later becomes a case from being interviewed again: SymptomaticParent() & !PreviouslyTraced(). By default such a case is interviewed like any other, on every transmission model.

Fields

EpiBranch.AlwaysEligible Type

Another name for TraceEveryone: every contact is traced.

EpiBranch.SymptomaticParent Type

julia

struct SymptomaticParent <: TraceEligibility

Trace contacts of cases that are symptomatic and isolated, starting from isolation. This is the default eligibility of ContactTracing. OnSymptomOnset() & OnIsolation() picks out the same cases but starts from the later of onset and isolation, so for a case isolated before its own onset the two trace from different times (see trigger_time).

Fields

EpiBranch.NoTracing Type

Another name for TraceNobody: no contact is traced.

EpiBranch.AnyOf Type

julia

struct AnyOf{T<:Tuple} <: TraceEligibility

Trace if any of the listed conditions holds. Usually written a | b.

Fields

conditions::Tuple

EpiBranch.AllOf Type

julia

struct AllOf{T<:Tuple} <: TraceEligibility

Trace only if all of the listed conditions hold. Usually written a & b.

Fields

conditions::Tuple

EpiBranch.NoneOf Type

julia

struct NoneOf{T<:Tuple} <: TraceEligibility

Trace only if none of the listed conditions holds. !a is the shorthand for a single condition.

Fields

conditions::Tuple

EpiBranch.TraceRate Type

julia

abstract type TraceRate

julia

TraceRate

Rule for whether an eligible contact is actually found by tracing. Built in: ConstantRate, a fixed probability per contact. To write a new rule, define a subtype and a method of traces.

Fields

EpiBranch.ConstantRate Type

julia

struct ConstantRate <: TraceRate

Each eligible contact is traced with the same probability p.

Fields

p::Float64

EpiBranch.TraceDelay Type

julia

abstract type TraceDelay

julia

TraceDelay

Rule for the delay, in days, from the case's tracing start (its isolation by default; see trigger_time) to its contact being reached. Built in: ConstantDelay, one distribution for every contact. To write a new rule, define a subtype and a method of draw_trace_delay.

Fields

EpiBranch.ConstantDelay Type

julia

struct ConstantDelay{D<:Distributions.Distribution} <: TraceDelay

Every contact's tracing delay (days) is drawn from the same distribution dist.

Fields

dist::Distributions.Distribution

EpiBranch.TraceAction Type

julia

abstract type TraceAction

julia

TraceAction

What happens to a contact once traced: Quarantine or FlagOnly. To write a new action, define a subtype and a method of apply_trace!.

Fields

EpiBranch.Quarantine Type

julia

struct Quarantine{D} <: TraceAction

julia

Quarantine(; duration)

Quarantine a traced contact from the time they are reached (or from their own isolation, if that is earlier), so they cannot infect others while quarantined. A trace that never arrives (trace time Inf) marks the contact as traced and quarantined without changing any isolation already in place.

duration (required) is how long the quarantine lasts, in days: a number, a distribution or a function of the random number generator and the individual, (rng, ind) -> ..., drawn once per trace. There is no default, because policies differ on how long quarantine should last. Inf keeps the contact quarantined until the end of their infectious period. A finite duration matters for a contact who escapes the traced exposure and is infected later through another route: that person transmits normally once the quarantine ends.

ind.state[:trace_time] is already set by the time duration is called. A rule phrased as an end day counted from exposure, such as "until 14 days after exposure", can subtract the trace time from that day to get a duration. A rule that ends quarantine at symptom onset can compare it against onset_time(ind). A contact traced more than once (ring propagation keeping the widest budget any infector gave it, or more than one child triggering SourceTracing on the same parent) keeps its earliest trace in :trace_time. On a later trace, that stored time predates the one duration is being called for, so a rule that subtracts it from a fixed end day will run longer than intended.

Fields

duration::Any

EpiBranch.FlagOnly Type

julia

struct FlagOnly <: TraceAction

Record the contact as traced without quarantining them: they keep transmitting until they would be isolated as a case. If they develop symptoms, Isolation then isolates them at whichever comes first, their own self-reported isolation or the later of the trace time and their symptom onset. An asymptomatic contact is never isolated this way.

Fields

Vaccination

Ring vaccination of traced contacts, mass vaccination of the population and vaccination of a case's group, with leaky or all-or-nothing protection.

EpiBranch.AbstractVaccination Type

julia

abstract type AbstractVaccination <: AbstractIntervention

Parent type of the vaccination interventions: RingVaccination (traced contacts), GroupVaccination (everyone in a case's group) and MassVaccination (the population on a schedule). They differ in who is vaccinated and when; the dose's effects are the same for all three and are set by these keywords (collected in a VaccineEffect):

Effect	What it does	Keyword
Against infection	Probability that an exposure of a vaccinated person does not infect them, once immune	`efficacy`
Against severe disease	Probability that a vaccinated person's own illness is milder, once immune	`severity_efficacy`
Delay	Days from vaccination to immunity	`delay_to_immunity`
Waning	Fraction of protection left a given number of days after immunity	`waning`
Leaky or all-or-nothing	How `efficacy` acts	`mode`
Dose name	Tells doses of a multi-dose schedule apart	`dose_label`

For extension authors: a new vaccination subtype decides eligibility and nothing else. It holds a VaccineEffect and returns it from vaccine_effect, the only route shared code takes to those parameters, and it shares the competing_risk that blocks a vaccinated contact's infection once their immunity has developed. An effect that only some vaccinations have, such as RingVaccination's post_exposure_efficacy, stays on the type that has it and draws its per-dose value in the record_effect_draws! hook.

RingVaccination adds two effects on infections a contact already has: post_exposure_efficacy (stopping the infection) and onward_efficacy (reducing onward transmission).

efficacy, severity_efficacy and delay_to_immunity each take a number, a distribution or a function of the random number generator and the individual, (rng, ind) -> .... A distribution or function is drawn once per vaccinated person when the dose is given, so the same person keeps the same efficacy and immunity delay against every exposure.

Leaky and all-or-nothing vaccines

With mode = LeakyMode() (the default) every exposure of a vaccinated person is prevented with probability efficacy. With mode = AllOrNothingMode() a share efficacy of vaccinated people are fully protected against infection from their immunity date on and the rest get no protection against infection. The mode acts on efficacy alone: severity_efficacy, post_exposure_efficacy and onward_efficacy apply to everyone vaccinated. AllOrNothingMode cannot yet be combined with waning (this raises an error).

The two modes differ only with repeated exposure

In a branching process each contact is a separate person exposed once, so the two modes give the same infection probability per contact (though not the same random draws). They differ in the mass-action, household and network models, where a person can be exposed more than once: under all-or-nothing a protected person stays protected against every exposure, while a leaky vaccine is worn down by repeated exposure.

Severity

severity_efficacy does not stop transmission. It acts only where a clinical transition in progression (such as Death or Hospitalisation) reads it through the severity_efficacy and immunity_time functions in its probability; RingVaccination has a worked example.

Waning

waning is a function of the days dt since the person became immune, giving the fraction of protection left (dt = 0 at immunity), or nothing (default) for constant protection. It multiplies efficacy (and, for RingVaccination, onward_efficacy and post_exposure_efficacy) at each exposure, scaling each person's own drawn efficacy from their own immunity date. It matters most when vaccination comes long before exposure. Post-exposure protection acts at the moment of immunity and so uses waning(0). waning does not apply to severity_efficacy; RingVaccination shows how to make that wane too.

Multi-dose schedules

Give each dose its own dose_label (default :default), for example :prime and :boost. Each dose then has its own vaccination record and wanes from its own immunity date. The doses act independently, so an exposure gets through all of them with probability prod(1 - eff_i * w_i), where eff_i is dose i's efficacy and w_i the fraction it retains: a prime at 0.6 and a boost at 0.7, both at full strength, together prevent 0.88 of exposures. is_vaccinated reports the dose with the default label; pass its dose_label keyword to check another.

Fields

EpiBranch.VaccineEffect Type

julia

struct VaccineEffect{E, SV, D, W, M<:AbstractEffectMode}

julia

VaccineEffect(; efficacy, severity_efficacy = 0.0, delay_to_immunity = 0.0,
    waning = nothing, mode = LeakyMode(), dose_label = :default)

The effects of a vaccine dose shared by every vaccination intervention, the same whoever receives it and whenever it is given. The built-in vaccinations take these as keywords and build the VaccineEffect themselves; build one directly only when writing a new vaccination intervention (see A custom vaccination).

efficacy: for a leaky vaccine, the probability that each exposure is prevented once immune; for an all-or-nothing vaccine, the probability of being fully protected.
severity_efficacy: probability that the vaccinated person's own illness is milder once immune (default 0).
delay_to_immunity: days from vaccination to immunity (default 0).
waning: a function of days since immunity giving the fraction of protection left, or nothing (default) for constant protection. Not yet available with mode = AllOrNothingMode().
mode: LeakyMode (default) or AllOrNothingMode.
dose_label: name of the dose, so several doses can be told apart (default :default).

The first three each take a number, a distribution or a function (rng, ind) -> ..., drawn once per vaccinated person when the dose is given. See AbstractVaccination for how the effects combine.

Fields

efficacy::Any
severity_efficacy::Any
delay_to_immunity::Any
waning::Any
mode::AbstractEffectMode
dose_label::Symbol

EpiBranch.RingVaccination Type

julia

struct RingVaccination{V<:VaccineEffect, C, DD, W, PE, OE} <: AbstractVaccination

julia

RingVaccination(; efficacy, coverage = 1.0, delay_to_immunity = 0.0,
                post_exposure_efficacy = 0.0, onward_efficacy = 0.0,
                severity_efficacy = 0.0, eligibility_window = Inf,
                dose_delay = 0.0, requires_dose = nothing, waning = nothing,
                mode = LeakyMode(), dose_label = :default)

Vaccinate the contacts found by ContactTracing (ring vaccination), including post-exposure prophylaxis (PEP, as in pepbp). Needs ContactTracing among the interventions.

Examples

julia

# Ring vaccination of contacts and contacts of contacts
[ContactTracing(OnSymptomOnset(), 0.8, Exponential(1.0), Quarantine(duration = 7.0); depth = 2),
 RingVaccination(efficacy = 0.9, delay_to_immunity = 10.0, coverage = 0.8)]

# Post-exposure prophylaxis: a dose that takes effect before onset
# stops the infection with probability 0.6
RingVaccination(efficacy = 0.0, post_exposure_efficacy = 0.6)

Arguments

efficacy: probability that each later exposure of the vaccinated contact is prevented once immune (see AbstractVaccination for leaky and all-or-nothing vaccines).
coverage: probability that a traced contact is actually vaccinated (consent, absence, exclusion criteria, logistics). Default 1. A number, a distribution or a function (rng, contact) -> ..., for example age dependent, or reading vaccine_acceptance to cluster refusal by group.
delay_to_immunity: days from vaccination to immunity. Default 0, as for PEP; set it for a vaccine that takes time to protect.
post_exposure_efficacy: probability that the dose stops an infection the contact already has (default 0). Needs incubation_period from clinical_presentation.
onward_efficacy: probability that each onward transmission by the vaccinated contact is prevented once immune (default 0).
severity_efficacy: probability that the vaccinated contact's own illness is milder once immune (default 0). See "Severity" below.
eligibility_window: days after exposure beyond which a contact is no longer vaccinated (default Inf), as in filovirus protocols where vaccination more than about 21 days after exposure is pointless. A number or a function (rng, contact) -> ....
dose_delay, requires_dose: for second and later doses (below).
waning, mode, dose_label: see AbstractVaccination.

delay_to_immunity, post_exposure_efficacy, onward_efficacy and severity_efficacy take the same forms as efficacy (a number, a distribution or a function (rng, ind) -> ...), drawn once per contact when the dose is given. dose_delay is drawn once when the dose is scheduled.

The window is checked on the day of vaccination, so a contact vaccinated just inside it with a long delay_to_immunity is still vaccinated, and whether immunity comes in time is then decided at each exposure. On a network or household model, a contact not yet infected when the dose is due is vaccinated as though the window were open, since their own exposure is not yet known; once the race settles their infection, the dose is withdrawn if the window turns out not to cover it, counting against no capacity from then on.

The window is measured to this dose, dose_delay included

A second dose is given dose_delay days after the trace, and the window is checked against that day. Copying a first dose's eligibility_window onto a booster therefore rejects almost every booster, since dose_delay alone usually exceeds the window. Leave a later dose at the default Inf unless the protocol has a deadline on that dose itself.

Post-exposure and onward effects

post_exposure_efficacy can stop an infection when immunity (vaccination plus delay_to_immunity) arrives after the exposure but before the contact's symptom onset. A stopped infection transmits as usual until immunity arrives and not at all afterwards. It never has symptoms, so nothing triggered by onset happens (isolation, tracing, onset-timed transitions), and no transition in the progression takes effect at or after that time, so there is no later hospitalisation, death or other outcome; transitions before it stand. A stopped infection is still a case: it is counted by chain_statistics and listed by linelist with no onset date and a date_infection_aborted column. Where immunity is already in place at the exposure, the dose instead prevents the infection with the same probability. That is all it does for a contact who never has symptoms, so such a contact exposed before immunity is not protected. post_exposure_efficacy therefore already covers every exposure efficacy would; set one or the other, since setting both counts the protection twice.

onward_efficacy leaves the contact's infection and illness in place but prevents each of their transmissions after immunity with this probability, whenever their onset is. It therefore also reduces onward cases from contacts already infected when vaccinated. post_exposure_efficacy and onward_efficacy act independently, giving a vaccine that stops some infections and reduces the infectiousness of the rest. Setting efficacy and onward_efficacy to the same value gives a vaccine that acts equally on susceptibility and infectiousness.

efficacy changes nothing when contacts are traced after their infector is isolated

By default ContactTracing reaches a contact once its infector has been isolated, and perfect Isolation already prevents every later transmission to that contact. A dose acting only through efficacy then leaves results unchanged, with or without quarantine (action = FlagOnly()). efficacy has infections to prevent only when a contact can still be infected after being traced: with leaky isolation (post_isolation_transmission > 0), when tracing starts before isolation (for example eligibility = OnSymptomOnset()), or in a depth > 1 ring through people who keep transmitting after being traced. onward_efficacy acts on the contact's own later transmission, which quarantine already prevents, so it matters when tracing does not quarantine; post_exposure_efficacy acts on the infection the contact already has. Under quarantine, post_exposure_efficacy can lower containment, because a stopped infection never has symptoms and so no longer triggers tracing of the people it infected before the dose.

Severity

severity_efficacy does not change transmission. It acts only where a transition in the progression reads it, through the severity_efficacy and immunity_time functions, for example a lower case fatality ratio among vaccinated cases:

julia

Death(delay = LogNormal(2.5, 0.4),
      probability = (rng, ind) ->
          immunity_time(ind) <= onset_time(ind) ?
              0.7 * (1 - severity_efficacy(ind)) : 0.7)

The check immunity_time(ind) <= onset_time(ind) gives protection only to a person immune by their onset; checking is_vaccinated alone would count a dose whose immunity had not yet developed.

severity_efficacy does not wane: waning applies only to the protection against infection, so with waning set, protection against infection fades while protection against severe outcomes stays. To make it fade too, apply the decay in the transition, from immunity to onset:

julia

decay(dt) = exp(-dt / 180)
Death(delay = LogNormal(2.5, 0.4),
      probability = (rng, ind) -> begin
          dt = onset_time(ind) - immunity_time(ind)
          dt >= 0 ? 0.7 * (1 - severity_efficacy(ind) * decay(dt)) : 0.7
      end)

dt >= 0 is the same immunity check, and is false for a person with no onset (NaN), who keeps the unvaccinated probability. Passing decay as the dose's waning too makes both effects fade at the same rate.

Second and later doses

A dose is given at the trace, or dose_delay days after it. requires_dose names a dose the contact must have received by the time this dose is due; a contact without it does not get this dose. A prime-boost schedule is two RingVaccinations:

julia

[
    RingVaccination(efficacy = 0.6, delay_to_immunity = 21.0,
        coverage = 0.8, dose_label = :prime),
    RingVaccination(efficacy = 0.5, dose_delay = Uniform(28.0, 42.0),
        delay_to_immunity = 14.0, coverage = 0.9,
        requires_dose = :prime, dose_label = :boost),
]

The boost is given on average five weeks after the trace to 90% of those primed (the other 10% lost to follow-up) and protects 14 days later. Doses act independently, so reaching 80% protection in total from a 60% prime needs efficacy = 0.5 on the boost ((0.8 - 0.6) / (1 - 0.6)), the protection it adds among those the prime left unprotected.

List a dose after the dose it requires: interventions are applied in order, so a boost listed first would find no prime and never be given. The required dose may come from any vaccination, such as a MassVaccination prime; a contact whose prime falls after the boost is due is not boosted. Between two ring doses, a boost dose_delay shorter than the prime's is rejected when the ModelSpec is built, since the boost could never be given. A dose_delay drawn from a distribution is rejected when even its longest delay falls before the prime's shortest, and gives a warning when the two ranges overlap, because contacts whose draws come out in the wrong order go without the boost. A function, or a distribution without a known range, cannot be checked in advance; a dose whose draw falls before the required dose is then skipped for that contact.

A dose is recorded when it is due, whether or not the contact was infected in the meantime, so counts of later doses are counts of doses scheduled.

In network and household models, infectors settle in order of infection, not of isolation, so a contact can be traced late from one infector and later traced earlier still from a second infector whose own isolation came first. A contact not yet settled has their dose moved to the earliest trace found rather than staying at the first; a settled contact's dose, and a dose another vaccination already gave, keep their date.

waning scales efficacy, post_exposure_efficacy and onward_efficacy by the time since this contact became immune. Post-exposure protection acts the moment immunity arrives and uses waning(0), so a waning that builds up first, such as dt -> min(1, dt / 14), stops no infections.

For extension authors: each contact's dose is recorded in ind.state under :vaccinated, :vaccination_time, :vaccine_efficacy, :immunity_time and :severity_efficacy (with the dose_label as a suffix for a non-default label), plus :post_exposure_efficacy and :onward_efficacy when those vary between people.

Fields

effect::VaccineEffect
coverage::Any
dose_delay::Any
requires_dose::Union{Nothing, Symbol}
eligibility_window::Any
post_exposure_efficacy::Any
onward_efficacy::Any

EpiBranch.MassVaccination Type

julia

struct MassVaccination{V<:VaccineEffect, T, C} <: AbstractVaccination

julia

MassVaccination(; efficacy, eligibility_time, coverage = 1.0,
                severity_efficacy = 0.0, delay_to_immunity = 0.0,
                waning = nothing, mode = LeakyMode(), dose_label = :default)

Vaccinate the population on a schedule, independently of contact tracing.

Arguments

eligibility_time: the day each person is vaccinated; Inf means never. A number (everyone on that day), a distribution such as Exponential(60.0) for a slow random rollout (drawn per person), or a function (rng, ind) -> ... for an age-stratified or other rule. Reading vaccine_acceptance here clusters refusal by group: for example (rng, ind) -> ind.state[:vaccine_acceptance] > 0.5 ? 30.0 : Inf leaves out whole groups whose acceptance is at most 0.5, so uptake is P(acceptance > 0.5).
coverage: probability that a person eligible at eligibility_time is actually vaccinated (consent, absence, logistics). Default 1. A number, a distribution or a function (rng, ind) -> ..., as for RingVaccination's coverage.
efficacy, delay_to_immunity, severity_efficacy: take the same forms, drawn once per vaccinated person, for example an age-dependent efficacy or immunity taking one to three weeks. See AbstractVaccination.
waning: decline of protection against infection from the day of immunity (see AbstractVaccination); matters most when people are vaccinated long before exposure. It does not apply to severity_efficacy, which acts only through a clinical transition that reads it, as for RingVaccination.
mode, dose_label: see AbstractVaccination. Use two MassVaccinations with different dose_labels for a multi-dose rollout.

Each contact's vaccination day is drawn when a branching process creates the contact (index cases are not vaccinated), and they count as vaccinated if it is finite. They are protected against an exposure only if their vaccination day plus delay_to_immunity comes before it.

Closed-form offspring distribution

A campaign already complete before the simulation starts (eligibility_time a number with eligibility_time + delay_to_immunity <= 0) protects every case alike, so the effect on the offspring distribution has a closed form: each secondary case is prevented independently with probability coverage * efficacy, and a NegativeBinomial(R, k) offspring distribution becomes NegativeBinomial(R * (1 - coverage * efficacy), k) (a Poisson likewise, on its mean). extinction_probability, probability_contain and the other analytical functions use this when coverage, efficacy and delay_to_immunity are numbers and waning === nothing; a rollout still under way during the simulation changes each case's exposure to vaccinated contacts over time, which has no such form, and those functions raise an error instead.

Examples

Whole population eligible on day 30:

julia

MassVaccination(efficacy = 0.85, eligibility_time = 30.0)

Per-individual rollout draws from a distribution, for a vaccine whose immunity takes one to three weeks to develop:

julia

MassVaccination(efficacy = 0.85,
    eligibility_time = Exponential(60.0),
    delay_to_immunity = Uniform(7.0, 21.0))

Age-stratified rollout (65+ on day 30, younger on day 90), with age-dependent efficacy:

julia

MassVaccination(
    efficacy = (rng, ind) -> ind.state[:age] >= 65 ? 0.7 : 0.9,
    eligibility_time = (rng, ind) -> ind.state[:age] >= 65 ? 30.0 : 90.0,
    delay_to_immunity = 14.0,
)

Prime-and-boost schedule (two doses with different labels):

julia

[
    MassVaccination(efficacy = 0.6,  eligibility_time = 30.0,
        delay_to_immunity = 14.0, dose_label = :prime),
    MassVaccination(efficacy = 0.9,  eligibility_time = 60.0,
        delay_to_immunity = 14.0, dose_label = :boost),
]

Fields

effect::VaccineEffect
eligibility_time::Any
coverage::Any

EpiBranch.GroupVaccination Type

julia

struct GroupVaccination{V<:VaccineEffect, E<:TraceEligibility, C, DD} <: AbstractVaccination

julia

GroupVaccination(; efficacy, eligibility = OnLabConfirmation(),
                 coverage = 1.0, dose_delay = 0.0, group_key = :group,
                 severity_efficacy = 0.0, delay_to_immunity = 0.0,
                 waning = nothing, mode = LeakyMode(), dose_label = :default)

Vaccinate everyone in a case's group (a community, health area or household) once a case there is confirmed or suspected: the usual fallback when no ring of contacts can be traced. Every member is vaccinated, whether or not they have any traced link to the case. Groups come from the groups population characteristic (stored under group_key, default :group).

Examples

julia

# Vaccinate the whole community 2 days after a case there is lab-confirmed
GroupVaccination(efficacy = 0.7, eligibility = OnLabConfirmation(), dose_delay = 2.0)

Arguments

eligibility: which cases trigger vaccination of their group, using the same rules as ContactTracing but applied to the case: OnLabConfirmation() (default) once a case has tested positive, OnSymptomOnset() on suspicion alone, combined with &, | and !. The rule's trigger_time sets when the group counts as having a case.
dose_delay: days from that trigger to vaccination, for example the time a vaccination team needs to reach the group. A number, a distribution or a function, drawn once per member, so members can be reached on different days.
coverage, efficacy, severity_efficacy, delay_to_immunity, waning, mode, dose_label: as for RingVaccination. Pairing coverage with vaccine_acceptance on the same group_key clusters refusal by group. severity_efficacy acts only through a clinical transition that reads it and does not wane.

A group is vaccinated at the trigger time of its earliest eligible case; later cases in the same group do not move the date. Members who appear later (a case infected later in the same group) are vaccinated at that same date plus dose_delay, so they are protected only from exposures after that, as for members already present. Doses grow with group size, where RingVaccination doses grow with ring size.

Group vaccination as a fallback to ring vaccination

Listing a RingVaccination before a GroupVaccination with the same dose_label makes the group dose a fallback: a member the ring has already vaccinated is skipped, so the group dose reaches only those the ring did not:

julia

[ContactTracing(OnLabConfirmation(), 0.7, Exponential(1.0), Quarantine(duration = 7.0)),
 RingVaccination(efficacy = 0.9),
 GroupVaccination(efficacy = 0.6)]

Needs the group characteristic and whatever eligibility needs, such as the test result from Isolation for the default OnLabConfirmation().

Unlike RingVaccination and MassVaccination, which act only on new contacts as they are made, GroupVaccination also checks a case's group once its own infection is settled, on every engine: an index case, or one whose own trigger (such as lab confirmation) only becomes true after it last made a contact, is still reached.

In network and household models, the dose can also reach the settling case itself and members not yet infected, subject to scheduling and capacity. Cases there are settled in order of infection, not of eligibility, so a member can turn out eligible earlier than the person who infected them. A member whose infection is not yet settled has their dose moved earlier if an earlier trigger is found later; a member already settled keeps their dose date, as does a dose another vaccination already gave.

For a campaign with repeat visits, coverage is the final probability of being reached, and dose_delay can describe the first successful visit. A person not covered is not marked as having refused: record willingness as a population characteristic when that distinction matters. See Repeat campaign visits for an example using these inputs.

Fields

effect::VaccineEffect
eligibility::TraceEligibility
coverage::Any
dose_delay::Any
group_key::Symbol

EpiBranch.AbstractEffectMode Type

julia

abstract type AbstractEffectMode

How a vaccine's efficacy acts: LeakyMode (each exposure is prevented with probability efficacy) or AllOrNothingMode (a share efficacy of vaccinated people are fully protected, the rest not at all).

For extension authors: a new mode defines EpiBranch.realised_efficacy, and EpiBranch.realise_prior_dose! only if a dose recorded before the simulation should be handled differently; see A custom effect mode.

Fields

EpiBranch.LeakyMode Type

julia

struct LeakyMode <: AbstractEffectMode

Leaky vaccine (the default): each exposure of a vaccinated person is prevented independently with probability efficacy.

Fields

EpiBranch.AllOrNothingMode Type

julia

struct AllOrNothingMode <: AbstractEffectMode

All-or-nothing vaccine: each vaccinated person is, with probability efficacy, fully protected against infection at every exposure from their immunity date on, and otherwise not protected against infection at all. The dose's other effects, such as severity_efficacy, apply to everyone vaccinated.

Fields

Timing and limited resources

Start or stop an intervention at a given time or case count, or cap how many people it can reach in each period.

EpiBranch.Scheduled Type

julia

struct Scheduled{I<:AbstractIntervention, F, T<:Tuple} <: EpiBranch.InterventionWrapper

julia

Scheduled(intervention; start_time = nothing, end_time = nothing,
          start_after_cases = nothing, start_after = nothing)
Scheduled(intervention, condition)

Start (and optionally stop) an intervention at a given day or once the outbreak has reached a given number of cases, for example a policy introduced two weeks into an outbreak.

Examples

julia

iso = Isolation(onset_to_isolation_delay = Exponential(2.0), duration = 7.0)

# Isolation from day 14
Scheduled(iso; start_time = 14.0)

# Tracing once there have been 50 cases
Scheduled(
    ContactTracing(
        probability = 0.5, isolation_to_trace_delay = Exponential(1.0),
        action = Quarantine(duration = 7.0),
    );
    start_after_cases = 50,
)

# Isolation between days 10 and 30
Scheduled(iso; start_time = 10.0, end_time = 30.0)

# Isolation once 20 cases have been reported
Scheduled(iso; start_after = ReportedCases(20))

# Isolation from seven days after the fifth reported case
Scheduled(iso; start_after = After(ReportedCases(5), 7.0))

# Any condition on the simulation, here from the third generation on
Scheduled(iso, state -> state.current_generation >= 3)

Arguments

start_time, end_time: days (as decimals, such as 14.0) between which the intervention acts, read against the day the action itself falls on (a case's isolation, a contact's trace) rather than the case's infection time, exactly as a scheduled vaccination's dose is checked against its own date.
start_after_cases: number of infections after which it acts, counting every infection as it happens, whether or not it is ever reported.
start_after: a trigger such as ReportedCases, which counts only cases reported by that day, or Infections. start_after_cases = n is the same as start_after = Infections(n). Wrap one in After to start a fixed delay after it opened, or write one that reads an admitted InterventionAction's own individual, through the three-argument is_triggered!, to gate it by group rather than by the whole population.
condition: a function of the simulation state returning true while the intervention should act.

Give at least one; several keywords must all hold. An isolation or trace dated before start_time, or after end_time, is undone with reset!: with start_time = 10.0, a case infected on day 8 whose symptoms begin on day 9, with a 2-day delay putting their isolation on day 11, is isolated after all, even though they were infected before the policy began.

Isolation and quarantine already in place end when the schedule ends (after end_time, or once condition returns false), whatever duration they were given, on every transmission model. Only protection the intervention keeps, such as vaccination, continues (see persistent_competing_risks). Everyone still records the intervention's starting information, such as "not isolated", before it starts.

A finite isolation or quarantine duration can end it before the schedule does, never later.

For extension authors: an intervention that proposes its actions through intervention_actions has each action checked against the schedule at the action's own time (seen by the condition as state.max_infection_time). Among the rest, one that defines its own intervention_time is run whether or not the schedule is active at the case's or contact's own clock, and its effect undone with reset! once intervention_time is checked against the condition the same way, with state.max_infection_time moved to that time for the check; a trigger that reads other population-wide state, such as ReportedCases, sees it as it stands now, not as it will once that time comes, the same approximation a scheduled vaccination's dose already accepts. An intervention relying on the default intervention_time (-Inf) is instead skipped outright while the schedule is inactive at that clock, as every hook always was: intervention_time comes paired with a reset! that undoes whatever it names, and without one there is nothing to correct the effect with afterwards.

Fields

intervention::AbstractIntervention
condition::Any
start_time::Float64
population_dependent::Bool
can_lapse::Bool
triggers::Tuple
time_aware::Bool

EpiBranch.is_active Function

julia

is_active(
    _::AbstractIntervention,
    _::SimulationState
) -> Union{Missing, Bool}

Whether an intervention is in effect at this point of the simulation (for example after a Scheduled start). Default: always.

EpiBranch.AbstractTrigger Type

julia

abstract type AbstractTrigger

julia

AbstractTrigger

Base type for a start_after trigger passed to Scheduled, as in Scheduled(iv; start_after = ReportedCases(20)). A subtype implements is_triggered! to say whether its condition holds now, and may override the trait trigger_can_lapse to say whether that can turn false again. Scheduled reads both to set its own population_dependent and can_lapse flags, so a trigger written outside the package works without editing Scheduled.

A trigger that also implements the three-argument is_triggered! can gate an admitted action by whom it is for rather than only by the population as a whole, and one that implements opened_at can be wrapped in After to open a fixed delay after its own condition first held.

Fields

EpiBranch.Infections Type

julia

struct Infections <: AbstractTrigger

julia

Infections(n)

A start_after trigger for Scheduled that opens once there have been n infections. Each infection counts as it happens, including one that is never reported or is reported later. Scheduled(iv; start_after = Infections(n)) is the same as Scheduled(iv; start_after_cases = n).

Fields

n::Int64

EpiBranch.ReportedCases Type

julia

struct ReportedCases <: AbstractTrigger

julia

ReportedCases(n)

A start_after trigger for Scheduled that opens once n cases have been reported by the current day of the simulation. A case reported after that day does not count yet, unlike in Infections. The model's progression needs a Reporting step to record the reports, except when n = 0, which opens from the start.

Fields

n::Int64

EpiBranch.opened_at Function

julia

opened_at(
    _::AbstractTrigger,
    _::SimulationState
) -> Union{Nothing, Real}

julia

opened_at(trigger, state::SimulationState) -> Union{Nothing, Real}

The day on which trigger turned true, if it can name one, or nothing. After reads this to start its delay from the moment a trigger opened, rather than from whenever the gate next happens to be checked, which can be well after the condition itself became true.

The default is nothing, the conservative answer for a trigger written outside the package: most triggers need not implement this to work with Scheduled directly, only to be wrapped in After. Infections and ReportedCases both answer exactly, from the n-th infection or report.

EpiBranch.After Type

julia

struct After{T<:AbstractTrigger} <: AbstractTrigger

julia

After(trigger, delay)

A start_after trigger for Scheduled that opens delay days after trigger itself opened, found through opened_at rather than from whenever the combined gate is next checked, which can land well after trigger's own condition first held.

Examples

julia

# Isolation from seven days after the fifth reported case
Scheduled(iso; start_after = After(ReportedCases(5), 7.0))

trigger must itself answer opened_at with something other than the conservative default nothing; Infections and ReportedCases both do. Building an After around a trigger that does not throws straight away, rather than leaving a gate that can silently never open.

Fields

trigger::AbstractTrigger
delay::Float64

EpiBranch.CapacityConstrained Type

julia

struct CapacityConstrained{I<:AbstractIntervention, F} <: EpiBranch.InterventionWrapper

julia

CapacityConstrained(intervention; budget_per_period, period=Inf,
                    carry_over=true, priority=default_capacity_priority)

Limit how many people an intervention can reach, such as a fixed number of vaccine doses per day: budget_per_period people every period days.

Examples

julia

# Ring vaccination with 5 doses a day
CapacityConstrained(RingVaccination(efficacy = 0.8); budget_per_period = 5.0, period = 1.0)

# Group vaccination with 200 doses in total
CapacityConstrained(GroupVaccination(efficacy = 0.8); budget_per_period = 200.0)

Arguments

budget_per_period: number of people that can be served per period.
period: length of a period in days, measured on the date each dose is actually given. Inf (default) gives a single budget for the whole outbreak.
carry_over: whether unused capacity is added to later periods (default true); with false each period has its own allowance.
priority: function (individual, state) -> number ordering who is served first, lowest first, with ties kept in order. The default, default_capacity_priority, serves traced contacts in the order they were traced.

Each case's actions are offered together, in the order the engine settles cases, and each is checked against the budget of the period its own proposed date falls in, whatever day the case itself was settled on. People turned away use no capacity and are not queued, though they may be offered the intervention again if reached again later. Only whole people are served, so a fraction of remaining capacity serves nobody. Interventions with the same dose label share one budget, and doses given by another intervention under that label count against it; use different labels for separate budgets. capacity_usage reports the capacity used and available.

Scheduled and CapacityConstrained can be combined in either order: both check each action's own proposed date.

Supported models

In network and household models, ring, group and mass vaccination can all be capacity-limited, with any eligibility_window: a dose a network or household race withdraws for falling outside the window (see RingVaccination) counts against no capacity from then on, and the budget it used becomes available to the next candidate. simulate warns that an intervention will have no effect when it relies on a hook these models do not call, such as one that reaches its targets only through the deprecated apply_post_transmission! batch hook; for a custom intervention, the warning is the thing to check. The mass-action models, HomogeneousProcess and MixingProcess, do not trace contacts, so ring vaccination cannot act in either. Contact tracing itself cannot be capacity-limited in any model.

For extension authors: a new intervention can be capacity-limited by defining intervention_actions and capacity_key. Group vaccination offers every eligible member of a triggered group as a separate person to serve. Interventions without intervention_actions, offered through the deprecated apply_post_transmission! batch hook, propose no dated action, so each is checked against the budget of the period the admitting call itself falls in rather than the day it is given.

Fields

intervention::AbstractIntervention
budget_per_period::Float64
period::Float64
carry_over::Bool
priority::Any

EpiBranch.capacity_usage Function

julia

capacity_usage(
    cc::CapacityConstrained,
    state::SimulationState
) -> NamedTuple{(:used, :available), <:Tuple{Int64, Any}}

julia

capacity_usage(cc::CapacityConstrained, state::SimulationState) -> (used, available)

The number of doses (or other units of capacity) dated by now, and the number that could have been given by now, at the end of a simulation or at the point state has reached. With carry_over = true (the default) both count the whole outbreak: available is budget_per_period for every period begun so far, and used counts every dose dated no later than it. With carry_over = false both refer to the current period only: available is one budget_per_period, and used counts the doses dated within it, including one given in that period by another intervention sharing the same capacity.

EpiBranch.default_capacity_priority Function

julia

default_capacity_priority(individual, state) -> Real

julia

default_capacity_priority(individual, state) -> Float64

Default order in which CapacityConstrained serves people: first come, first served by the day they were traced. Someone who was not traced comes last, after everyone traced at the same point of the simulation.

Transmission routes

A route of transmission, such as community, household or funeral contact, with its own start, end and timing during a case's course of infection.

EpiBranch.RouteWindow Type

julia

struct RouteWindow{K, R}

julia

RouteWindow(name; from = nothing, until, kernel, reach = name,
            contacts_from = :infection, traceable = 1.0)

A route of transmission, such as community, household or funeral, with its own start, end and timing over a case's natural history. Used with RoutedNetwork (in EpiNetwork), where each route reaches a different set of people in the network.

name labels the route (:community, :household, :funeral, …).
from is the state at which this route's infectiousness begins. :infection opens it at the infection time itself; any other state opens it at that state's time, following the <state>_time convention that Transition writes. The default, nothing, opens it at :infectious when the natural history has a step to it, otherwise at infection.
until is the tuple of states that end the route; it closes at the earliest of them. A state ends only the routes that list it, so a control measure can end one route and leave another open. List EpiBranch.INTERVENTION_REMOVAL for a route that the model's interventions (isolation, or quarantine after being traced) should end.
kernel is the contact interval on this route: the time in days from the route opening to a contact that would infect if the route were still open.
reach tells the model who the route reaches. RoutedNetwork needs the route's contact list here (reach[i] lists the people i can infect on this route); a model of your own can read another value, such as the default, name.
contacts_from is the state from which the people the route reaches are the case's contacts, which is what contact tracing acts on. The default, :infection, suits a route over standing relationships such as a household, whose members are contacts however early the case is isolated. A route whose contacts come about through an event names that event's state, such as :died for a funeral: its contacts exist only if the event happens before the route is cut, and cannot be traced before it.
traceable is the probability that a case can name a given contact made on this route, so that contact tracing can find it. People can name the people they live with but not the strangers they stood next to, so a household route might keep the default 1.0 and a community route of casual encounters take something much lower. true and false also work, as 1.0 and 0.0. Naming comes before the tracing policy: a contact is traced only if the case names it and ContactTracing then traces it, so the overall probability is traceable times the tracing probability. Use traceable for what the relationship allows and the tracing probability for how well the programme performs, and count each limit in only one of them. Each model applies the probability to its contacts in its own way (see EpiNetwork.RoutedNetwork).

Examples

In these examples community_contacts, household_contacts and funeral_contacts are contact lists, one vector of people per person.

A case that stops mixing in the community when it isolates but keeps infecting the people it lives with:

julia

community = RouteWindow(:community; from = :infectious,
    until = (:recovered, EpiBranch.INTERVENTION_REMOVAL), kernel = Exponential(4.0),
    reach = community_contacts)
household = RouteWindow(:household; from = :infectious,
    until = (:recovered,), kernel = Weibull(1.5, 3.0), reach = household_contacts)

Ebola transmission at funerals, a route that only opens once the case has died and closes at burial:

julia

funeral = RouteWindow(:funeral; from = :died, until = (:buried,),
    kernel = Exponential(1.0), contacts_from = :died, reach = funeral_contacts)

A case that recovers never opens the funeral route, and with contacts_from = :died contact tracing also treats a survivor as having no funeral contacts.

A community route whose contacts are mostly strangers, of whom a case can name one in five:

julia

community = RouteWindow(:community;
    until = (:recovered, EpiBranch.INTERVENTION_REMOVAL),
    kernel = Exponential(4.0), traceable = 0.2, reach = community_contacts)

Fields

name::Symbol
from::Union{Nothing, Symbol}
until::Tuple
kernel::Any
reach::Any
contacts_from::Symbol
traceable::Float64

EpiBranch.window_open Function

julia

window_open(ind::Individual, w::RouteWindow) -> Any

julia

window_open(individual, window)

Time (days) at which transmission by route window starts for a case, or Inf if the case never reaches its from state, in which case the route transmits nothing. A route with from = nothing opens when the case becomes infectious if the natural history has that step, otherwise at infection.

EpiBranch.window_close Function

julia

window_close(ind::Individual, w::RouteWindow) -> Any
window_close(
    ind::Individual,
    w::RouteWindow,
    interventions
) -> Any

julia

window_close(individual, window, interventions = ())

Time (days) at which transmission by route window stops for a case: the earliest of its until states, or Inf if none has been reached. Each route is ended separately, so the same event can end one route and leave another open.

A route listing EpiBranch.INTERVENTION_REMOVAL also closes when the interventions remove the case for good (perfect isolation, or quarantine, with no end); leaky or time-limited isolation acts on each contact instead. Pass the interventions the simulation used to get the same closing time as the simulation.

EpiBranch.INTERVENTION_REMOVAL Constant

julia

INTERVENTION_REMOVAL

List this in a RouteWindow's until to let control measures end transmission along that route: the route then stops when the model's interventions (isolation, quarantine after tracing) remove the case.

Other states in until come from the natural history, so a route ends when the case recovers, dies or is buried. Removal by an intervention is different because whether it removes the case at all depends on the intervention: perfect isolation takes the case out of transmission entirely, whereas leaky isolation only reduces it. Listing INTERVENTION_REMOVAL asks each intervention when it removes the case (through infectious_removal_time).

This is why routes are worth separating. A community route lists it, so isolating a case ends their community transmission; a household route does not, so the case goes on infecting the people they live with. The same applies to the reduction from leaky isolation, but not to a vaccine's protection, which applies on every route; see risk_applies.

Running simulations

simulate runs one outbreak or many. Stopping rules end a run early, for example after a maximum number of cases or days. A run returns a SimulationState holding one Individual per person in the simulation: the cases and the contacts they exposed for a BranchingProcess, everyone for models with a fixed population.

EpiBranch.simulate Function

julia

simulate(
    model::TransmissionModel;
    n_initial,
    initial_cases,
    max_cases,
    max_generations,
    max_time,
    stopping_rules,
    rng,
    condition,
    max_attempts
) -> Any

julia

simulate(model::TransmissionModel;
         max_cases=10_000, max_generations=100, max_time=nothing,
         n_initial=nothing, initial_cases=nothing, stopping_rules=nothing,
         rng=Random.default_rng(), condition=nothing, max_attempts=10_000)

Simulate one outbreak and return it as a SimulationState. Pass the result to linelist, contacts, chain_statistics or weekly_incidence to get tables.

model is a transmission model such as a BranchingProcess, or a ModelSpec that adds natural history, interventions, population characteristics and reporting to it. A model on its own runs with none of these. Times are in the units of the model's delays; the documentation uses days throughout.

Keywords

n_initial: number of index cases (default 1).
max_cases, max_generations, max_time (days): stop the run once any of these is reached. Set one to nothing to remove that limit. The run always stops when transmission dies out. The built-in models settle cases in order of infection time, so a run stops at exactly max_cases cases (or with the index cases, if there are more of them), leaves out any infection after max_time, and draws no contacts for a case at generation max_generations. A model written against generate_offspring or contacts_of runs one generation at a time and finishes the generation that crosses a limit, so it can end past max_cases or with infections after max_time. The last cases of a run the built-in models stop by max_cases or max_time have not made all their contacts, so per-case quantities such as offspring counts or realised generation intervals are right-truncated; max_generations alone leaves every case below the bound complete. On a branching process, a contact whose infection was still to be decided when the run stopped is left out, unless tracing or an intervention's actions had already reached it; it is then kept with no infection time.
stopping_rules: a vector of AbstractStoppingRules for finer control over when the run ends. When given, max_cases, max_generations and max_time are ignored.
rng: the random number generator; pass a seeded one, such as Xoshiro(42), for reproducible runs. Pass CommonRandomNumbers instead to compare scenarios on shared random numbers: two runs built from the same seed then share every draw until one of them changes what happens to a particular person or pair.
condition: a range of final outbreak sizes, such as 10:1000. The simulation is repeated until the total number of cases falls inside it (at most max_attempts tries), for example to leave out outbreaks that died out after a few cases. An error is raised if no run qualifies.
initial_cases: for NetworkProcess, RoutedNetwork and HouseholdProcess only, the IDs of the people infected at time zero (an empty vector is allowed). It replaces the random choice of index cases. Give either initial_cases or n_initial, not both; IDs outside the population raise an error. With an external_hazard, everyone else can still be infected from outside from time zero. Other models raise an error when given a vector here.

The probability and delay of each step in a case's natural history are set on the transitions in the progression of a ModelSpec. Each takes a constant or a function of the random number generator and the individual, (rng, ind) -> value, so rates and delays can depend on age or risk group.

Note

On a model over a fixed population (HomogeneousProcess, MixingProcess, NetworkProcess, RoutedNetwork and HouseholdProcess) the max_cases and max_generations keyword defaults have no effect, so a run over a large population is not capped unless you ask for it. A stopping_rules entry, or a non-default max_cases or max_generations, still ends such a run.

Examples

julia

using EpiBranch, Distributions, Random

model = BranchingProcess(NegBin(2.5, 0.16), Gamma(2.5, 2.0))
state = simulate(model; max_cases = 500, rng = Xoshiro(1))
linelist(state)

# keep only outbreaks with between 10 and 500 cases
state = simulate(model; max_cases = 500, condition = 10:500, rng = Xoshiro(1))

julia

simulate(
    model::TransmissionModel,
    n::Int64;
    n_initial,
    initial_cases,
    max_cases,
    max_generations,
    max_time,
    stopping_rules,
    rng,
    parallel
) -> Vector

julia

simulate(model, n::Int; parallel=false, kwargs...)

Simulate n independent outbreaks from the same model and return them as a vector of SimulationStates. It takes the same stopping keywords (max_cases, max_generations, max_time, stopping_rules), n_initial, initial_cases and rng as the single-outbreak method; condition is not available here. The usual next step is containment_probability, the share of these outbreaks that died out.

With parallel = true the runs are spread over the CPU threads Julia was started with (for example julia --threads 4). The results are reproducible for a given seeded rng. When Julia has more than one thread, they are not the same outbreaks as with parallel = false.

Examples

julia

using EpiBranch, Distributions, Random

model = BranchingProcess(NegBin(2.5, 0.16), Gamma(2.5, 2.0))
states = simulate(model, 1000; max_cases = 5000, rng = Xoshiro(1))
containment_probability(states; max_cases = 5000)

EpiBranch.SimOpts Type

julia

struct SimOpts

julia

SimOpts(; n_initial, initial_cases, max_cases, max_generations, max_time, stopping_rules)

Settings for how a simulation starts (n_initial index cases, or the initial_cases IDs) and when it stops. simulate takes the same keywords and builds this for you. Natural history, population characteristics and interventions are set with a ModelSpec.

The run stops at the first generation at which any of the stopping_rules (AbstractStoppingRule) says so, and always when transmission dies out (Extinction). max_cases, max_generations and max_time (days) are shortcuts that create the matching rules; they are ignored when stopping_rules is given.

julia

SimOpts(max_cases = 500)   # [Extinction(), MaxCases(500), MaxGenerations(100)]
SimOpts(max_generations = 20, max_time = 90.0)
SimOpts(stopping_rules = [MaxCases(1000), MyCustomRule()])

Fields

n_initial::Int64
initial_cases::Union{Nothing, Vector{Int64}}
stopping_rules::Vector{AbstractStoppingRule}

EpiBranch.AbstractStoppingRule Type

julia

abstract type AbstractStoppingRule

julia

AbstractStoppingRule

A rule for ending a simulation. The run stops as soon as any of its rules says so; pass them to simulate as stopping_rules, or use the max_cases, max_generations and max_time keywords, which create the first three of these:

MaxCases: stop once the outbreak reaches a number of cases.
MaxGenerations: stop after a number of generations.
MaxTime: stop once infections reach a time, in days.
Extinction: stop when transmission has died out. Always included, even when you give your own stopping_rules.

To write your own rule, define a new type and a method EpiBranch.should_stop(rule::MyRule, state::SimulationState) (the EpiBranch. prefix and the ::SimulationState are both needed):

julia

struct MaxChainLength <: AbstractStoppingRule
    n::Int
end
EpiBranch.should_stop(r::MaxChainLength, state::SimulationState) =
    maximum(ind.generation for ind in state.individuals; init = 0) >= r.n

The mass-action, network and household models run a continuous-time race rather than stepping by generation; they poll should_stop after each settled case, for any rule whose effect is not already captured by time_bound and generation_bound alone, as MaxTime and MaxGenerations are. A rule bounded this way also defines honoured_without_should_stop. The max_cases and max_generations keyword defaults have no effect on these models, so a run over a large fixed population is not capped unless you ask: give an explicit value, or a stopping_rules entry, to end such a run early. The Extending EpiBranch guide has a worked example in its section on stopping rules.

Fields

EpiBranch.Extinction Type

julia

struct Extinction <: AbstractStoppingRule

Stop when transmission has died out (no one is left who can still infect others).

Fields

EpiBranch.MaxCases Type

julia

struct MaxCases <: AbstractStoppingRule

Stop once the outbreak has at least n cases in total.

Fields

n::Int64

EpiBranch.MaxGenerations Type

julia

struct MaxGenerations <: AbstractStoppingRule

Stop once n generations have been simulated.

Fields

n::Int64

EpiBranch.MaxTime Type

julia

struct MaxTime <: AbstractStoppingRule

Stop once the latest infection is at or after time t (days since the start of the outbreak).

Fields

t::Float64

EpiBranch.SimulationState Type

julia

mutable struct SimulationState{T<:Real, R<:Random.AbstractRNG, P, A}

One simulated outbreak, as returned by simulate. Turn it into tables with linelist (one row per case), contacts (one row per exposure), chain_statistics or weekly_incidence.

Useful fields:

individuals: one Individual per person in the simulation: the cases and the contacts they exposed for a BranchingProcess, everyone for models with a fixed population.
cumulative_cases: the total number of cases.
extinct: whether transmission had died out when the run stopped.
current_generation: the last generation simulated.

For extension authors: transitions holds the natural-history steps of the model, and scratch is space in which an intervention can keep its own working information for the run (as Individual.state is for each person), named as the Extending EpiBranch guide describes. The simulation never reads scratch. GroupVaccination keeps its list of each group's members there.

Fields

individuals::Array{Individual{T}, 1} where T<:Real
active_ids::Vector{Int64}
current_generation::Int64
rng::Random.AbstractRNG
cumulative_cases::Int64
extinct::Bool
population_size::Any
max_infection_time::Real
attributes::Any
transitions::Vector{AbstractClinicalTransition}
scratch::Dict{Any, Any}

EpiBranch.Individual Type

julia

mutable struct Individual{T<:Real}

One person in a simulation: for a BranchingProcess, a case or a contact they exposed; for models with a fixed population, anyone in it.

id: the person's number, which is also their position in state.individuals.
parent_id: the id of their infector (0 for an index case).
generation: 0 for index cases, 1 for the people they infect, and so on.
chain_id: which index case's transmission chain they belong to.
infection_time: when they were exposed, in days since the start of the outbreak.
susceptibility, between 0 and 1: the probability that they are infected when exposed.
infectiousness, between 0 and 1: the factor applied to their onward transmission once infected (each of their contacts is infected only with this probability).
secondary_case_ids: the ids of the people they exposed. A branching process lists the contacts it created, infected or not. One whose exposure did not happen within the run has a NaN infection time: it falls after the run stopped, or a ring created it for a member who was never infected. Most contacts a run stopped before reaching were never created and are not listed. A model written against generate_offspring lists every contact it created; one written against contacts_of lists an existing person only if it infected them. The fixed-population models below list only the people they infected.
state: everything else recorded about them, such as symptom onset (:onset_time), :asymptomatic, :age, :sex, whether they were :isolated, :traced, :quarantined or :vaccinated, :test_positive, their :type in a multi-type model, whether they were :infected, and any population characteristics you add. It is read and written as a dictionary (see EpiBranch.IndividualState).
episodes: past infections, oldest first, for a person who has been infected more than once (see close_episode!). The fields above describe the current or latest infection. None of the built-in models reinfects anyone, so this is empty there.

In the fixed-population models (HomogeneousProcess, MixingProcess, NetworkProcess, RoutedNetwork, HouseholdProcess), susceptibility and infectiousness instead multiply the rate at which a pair meets. They then lower the chance of infection only over an infectious period of limited length: if the infectious period never ends, any positive value leads to infection eventually, only later on average, and a contact certain to fall within the infectious period (a fixed delay such as Dirac(2.0)) infects whatever the value.

Setting fields at simulation time

Give attributes to a ModelSpec; they are drawn for every new person when they are created. Several can be listed, and are applied in order:

julia

attributes = [
    clinical_presentation(incubation_period = LogNormal(1.6, 0.5)),
    demographics(age_distribution = Uniform(0, 90)),
    transmission_traits(susceptibility = 0.3, infectiousness = 0.9),
]

For anything else, add your own function of the random number generator and the individual, (rng, ind) -> ..., to the list.

Fields

id::Int64
parent_id::Int64
generation::Int64
chain_id::Int64
infection_time::Real
susceptibility::Real
infectiousness::Real
secondary_case_ids::Vector{Int64}
state::EpiBranch.IndividualState
episodes::Array{InfectionEpisode{T}, 1} where T<:Real

EpiBranch.IndividualState Type

julia

struct IndividualState{T<:Real} <: AbstractDict{Symbol, Any}

julia

IndividualState{T}

What an Individual's state holds: everything recorded about one person beyond the core fields, keyed by Symbol. It is an AbstractDict{Symbol, Any}, so it is read and written as a dictionary (ind.state[:onset_time], get(ind.state, :age, missing), haskey(ind.state, :traced)), and any AbstractDict{Symbol, Any} passed as state when an Individual is built is copied into one.

Each value comes back as the type it was stored with, and keys are kept in the order they were first set.

Fields

slots::Array{EpiBranch._Slot{T}, 1} where T<:Real

EpiBranch.InfectionEpisode Type

julia

struct InfectionEpisode{T<:Real}

A past infection of a person who has since been reinfected, kept on the Individual. It records that infection's infection_time, its place in the transmission tree (parent_id, generation, chain_id), the secondary cases it caused, and the person's state as it was when that infection ended (natural-history times, outcome, and what interventions recorded).

See also susceptible_again_time, the time from which a person can be infected again, and close_episode!.

Fields

infection_time::Real
parent_id::Int64
generation::Int64
chain_id::Int64
secondary_case_ids::Vector{Int64}
state::EpiBranch.IndividualState

EpiBranch.CommonRandomNumbers Type

julia

struct CommonRandomNumbers{R<:Random.AbstractRNG} <: Random.AbstractRNG

julia

CommonRandomNumbers(seed)

An AbstractRNG, built from seed, that keeps two runs on the same draws for everything routed through stream — progression delays, infectiousness and susceptibility, acceptance of an offer, and a pair's contact times — until one of them changes what actually happens to a particular person or pair. Pass it as simulate's rng in place of a plain seeded generator to compare scenarios on shared random numbers:

julia

baseline = simulate(ModelSpec(net; progression); rng = CommonRandomNumbers(12))
compared = simulate(
    ModelSpec(net; progression, interventions); rng = CommonRandomNumbers(12)
)

A draw not routed through stream — an offspring count, which node is seeded first, anything a custom attribute or model reads from rng directly — still comes from one shared sequence, exactly as for any other AbstractRNG, so running the same model twice with the same CommonRandomNumbers reproduces it exactly either way.

Fields

seed::UInt64
rng::Random.AbstractRNG
substreams::Dict{Any, Random.Xoshiro}

EpiBranch.stream Function

julia

stream(
    rng::Random.AbstractRNG,
    id,
    purpose
) -> Random.Xoshiro

julia

stream(rng, id, purpose) -> AbstractRNG

The generator to draw from for id (a person's id, or (a, b) for an ordered pair of them) and purpose (a Symbol or tag telling apart the different things drawn for the same id, such as :susceptibility and :infectiousness for one person, or :contact for one pair). The default method, for any ordinary AbstractRNG, returns rng itself: every draw still comes from the one sequence, in the order the engine asks for it, exactly as before stream existed.

CommonRandomNumbers is the opt-in alternative: it returns a generator kept just for (id, purpose), and a run built from it draws the same values for that id and purpose whatever else in the run has changed. A new transmission model or intervention that draws its own per-person or per-pair values should call stream(rng, id, purpose) for the generator to draw from, rather than rng directly, to take part.

Recording contacts

The continuous-time models (networks, households, HomogeneousProcess and MixingProcess) stop simulating contacts with a person an intervention has protected for good, such as a vaccinated contact. A contact recorder keeps simulating them and can log each one, through records_contacts. For a BranchingProcess, contacts already lists exposures that did not infect; for a continuous-time model, attach an ExposedContactRecorder to get the same, or write your own recorder type for anything else you want logged.

EpiBranch.ContactRecorder Type

julia

abstract type ContactRecorder

Asks the continuous-time models (networks, households, HomogeneousProcess and MixingProcess) to keep simulating contacts between an infector and a person whom an intervention has protected for good, such as a vaccinated contact. By default those contacts are no longer simulated, since they can no longer cause infection; keep them when you want every contact event recorded. Tracing and ring vaccination do not depend on it.

Attach one to a ModelSpec as recorder = .... A new recorder type defines records_contacts.

Fields

EpiBranch.NoContactRecorder Type

julia

struct NoContactRecorder <: ContactRecorder

The default: no contact recorder. Contacts between an infector and a person protected for good are no longer simulated (records_contacts is false for every pair).

Fields

EpiBranch.ExposedContactRecorder Type

julia

struct ExposedContactRecorder <: ContactRecorder

Logs every distinct infector-contact pair a continuous-time model's race blocks, whether or not the block is certain to recur, keeping the time of the first block. A pair logged from a block that is not certain to recur is dropped again if the same pair later infects: the block was not this pair's last outcome, so it is no longer an exposure that escaped. A blocked community introduction, where the person proposed as their own infector stands in for a source outside the population, is not a contact between two people and is not logged. Attach one to a ModelSpec as recorder = ExposedContactRecorder(); contacts then lists each pair still logged as an exposure that did not lead to infection, the same way it already does for a BranchingProcess.

Unlike a recorder whose records_contacts answers true, this one never asks the race to keep drawing a pair whose block is certain, so attaching it adds no restriction on which models can run.

Fields

Information about each case

Read what happened to a simulated case: when symptoms began, whether and when it was isolated, traced, quarantined or vaccinated, and how its infection ended.

EpiBranch.onset_time Function

julia

onset_time(ind::Individual{T}) -> Any

julia

onset_time(ind)

Time of symptom onset, in days since the start of the outbreak (NaN if the person is asymptomatic or has no onset time).

The functions that read a person's information, such as onset_time, is_isolated or is_infected, work on the people in a simulated outbreak:

julia

state = simulate(model)
count(is_isolated, state.individuals)                 # number isolated
onsets = [onset_time(i) for i in state.individuals if is_infected(i)]

EpiBranch.incubation_period Function

julia

incubation_period(ind::Individual) -> Any

julia

incubation_period(ind)

Incubation period: days from infection to symptom onset (NaN if the person is asymptomatic or has no onset time). Useful inside a generation_time function that links a case's generation time to its own incubation period.

EpiBranch.is_isolated Function

julia

is_isolated(ind::Individual) -> Bool

Whether the person was isolated in a way that counts as detecting them, so that tracing and group vaccination start from them and the line list shows it. An isolation that does not count as a detection (see EpiBranch.records_isolation), such as one after death, still stops transmission from isolation_time but leaves this false.

EpiBranch.isolation_time Function

julia

isolation_time(ind::Individual{T}) -> Any

Time from which isolation or quarantine stops the person transmitting, in days since the start of the outbreak (Inf if never), whether or not it counts as a detection (see is_isolated).

EpiBranch.isolation_release_time Function

julia

isolation_release_time(ind::Individual{T}) -> Any

Time at which the isolation or quarantine that began at isolation_time ends, in days since the start of the outbreak (Inf if it never does). It is finite when Isolation or the Quarantine of ContactTracing has a duration, so a quarantine that ended before the person was later infected does not stop their own onward transmission. Set by set_isolated!.

EpiBranch.outcome_time Function

julia

outcome_time(ind::Individual{T}) -> Any

julia

outcome_time(ind)

Time of the case's outcome, such as recovery or death, in days since the start of the outbreak: the earliest step of the natural history that ends the case (see Transition). Inf if there is none yet, or if the natural history has no such step.

EpiBranch.is_traced Function

julia

is_traced(ind::Individual) -> Bool

Whether the person was found by contact tracing.

EpiBranch.is_quarantined Function

julia

is_quarantined(ind::Individual) -> Bool

Whether the person was quarantined as a traced contact.

EpiBranch.is_vaccinated Function

julia

is_vaccinated(ind::Individual; dose_label) -> Bool

Whether the person received the vaccine dose named dose_label (by default, the single dose of a vaccination without a dose_label).

EpiBranch.immunity_time Function

julia

immunity_time(ind::Individual{T}; dose_label) -> Any

Time at which immunity from the vaccine dose named dose_label develops, in days since the start of the outbreak: the vaccination time plus the dose's delay_to_immunity (Inf if not vaccinated). A natural-history step can compare it with the time an outcome would take effect (for example onset_time(ind) for the default Death) to decide whether the dose's severity_efficacy applies: immunity that develops later gives no protection.

EpiBranch.vaccine_efficacy Function

julia

vaccine_efficacy(ind::Individual; dose_label) -> Any

The full-strength efficacy against infection drawn for this person's dose named dose_label, as it was sampled when the dose was given: the probability that an exposure is blocked under LeakyMode, or 1.0 or 0.0 for whether they responded under AllOrNothingMode. nothing if no dose with that label was recorded, which is what tells a person with no such dose apart from one given a dose of efficacy = 0.0. severity_efficacy draws no such distinction: it returns 0.0 for a person who was never dosed.

EpiBranch.severity_efficacy Function

julia

severity_efficacy(ind::Individual; dose_label) -> Float64

Vaccine efficacy against severe outcomes for this person: the proportional reduction in, for example, their probability of death once their immunity has developed (0 if not vaccinated, or if the dose has no effect on severity). Drawn once per vaccinated person, alongside efficacy, from the severity_efficacy of the dose's VaccineEffect. Use it in the probability of a natural-history step, once immunity_time has passed:

julia

Death(delay = LogNormal(2.5, 0.4),
      probability = (rng, ind) ->
          immunity_time(ind) <= onset_time(ind) ?
              0.7 * (1 - severity_efficacy(ind)) : 0.7)

julia

severity_efficacy(v::AbstractVaccination) -> Any

Efficacy of the vaccination's dose against severe disease, as given to VaccineEffect. The method for an individual returns the value drawn for that person.

EpiBranch.is_asymptomatic Function

julia

is_asymptomatic(ind::Individual) -> Bool

Whether the case never develops symptoms.

EpiBranch.is_test_positive Function

julia

is_test_positive(ind::Individual) -> Bool

Whether the case tested positive, as drawn by Isolation from its test_sensitivity.

EpiBranch.is_infected Function

julia

is_infected(ind::Individual) -> Bool

Whether the person was infected, as opposed to an exposed contact who escaped infection.

EpiBranch.susceptible_again_time Function

julia

susceptible_again_time(ind::Individual{T}) -> Any

Time from which a person who has been infected can be infected again, because their immunity has waned, in days since the start of the outbreak (Inf, never, by default). Until then EpiBranch.HostImmunity prevents reinfection. It matters only for a model that exposes people again after infection through its own contacts_of; none of the built-in models does.

Set it with a step in the model's progression, timed from the state that starts waning, for example Transition(:susceptible_again, from = :recovered, delay = Exponential(180)).

EpiBranch.individual_type Function

julia

individual_type(ind::Individual) -> Int64

Which type a person is in a multi-type branching process, as a number (1 in a single-type model).

Output

Tables and summaries of simulated outbreaks: line lists, contact lists, transmission chain statistics, containment probability, weekly incidence and realised generation intervals.

EpiBranch.linelist Function

julia

linelist(
    state::SimulationState;
    reference_date,
    infected_only
) -> Any

julia

linelist(state::SimulationState; reference_date=Date(2020, 1, 1),
         infected_only=true)

The line list of a simulated outbreak: a DataFrame with one row per case, with event times converted to calendar dates counted from reference_date (simulation day 0).

Columns always present:

Column	Meaning
`id`	case identifier
`parent_id`	`id` of the infector (0 for an index case)
`generation`	generation number (index cases are generation 0)
`chain_id`	which index case the case descends from
`date_infection`	date of infection

Further columns depend on the model and are added whenever any case has the information: for example date_onset and asymptomatic from clinical_presentation, age and sex from demographics, isolated and date_isolation from Isolation (with date_isolation_release when isolation or quarantine has a finite duration), date_reporting from Reporting or PerCaseObservation, and outcome and date_outcome from Death and Recovery. Any time recorded under a name ending in _time becomes a date_ column (:onset_time becomes date_onset); other information appears as recorded. Names starting with an underscore are internal to an intervention and never shown. To add your own date columns see event_time_metadata.

isolated, date_isolation and date_isolation_release show an isolation only when it counted as a detection (see is_isolated).

Including people who were not infected

With infected_only = false, the table has a row for every person in the simulation and an extra infected column. For a model with a fixed population, such as NetworkProcess or HouseholdProcess, that is everyone; for BranchingProcess it is the cases plus every contact they exposed who was not infected.

An uninfected person has missing for date_infection and every date that depends on infection (onset, reporting, admission, outcome, custom dates that require infection). Dates of events that can happen without infection are kept: date_trace, date_vaccination, date_immunity, and, for a traced contact put in quarantine, date_isolation (with date_isolation_release when the isolation or quarantine has a finite duration). If Isolation later replaced a quarantine with an isolation based on an onset that never happened, these two columns show the original quarantine, or missing if there was none.

Examples

julia

using Dates
model = ModelSpec(
    BranchingProcess(NegBin(2.5, 0.16), Gamma(2.0, 3.0));
    attributes = clinical_presentation(incubation_period = LogNormal(1.6, 0.5))
)
state = simulate(model; max_cases = 100)
ll = linelist(state; reference_date = Date(2024, 3, 1))
first(ll, 5)

EpiBranch.event_time_metadata Function

julia

event_time_metadata(
    _::Val{key}
) -> @NamedTuple{column::Symbol, requires_infection::Bool}

julia

event_time_metadata(::Val{key})

How a recorded event time appears as a date column in linelist. Use it when a custom intervention or progression step records its own event times and you want them shown as dates.

Returns (column = :date_name, requires_infection = true), or nothing for a column that is not a date. By default, any time recorded under a name ending in _time becomes a date_ column (:onset_time becomes date_onset) and is shown only for infected people. Times that are infinite or not numbers become missing.

An event that can happen to people who were never infected (an appointment, say) declares that its date applies to them too:

julia

EpiBranch.event_time_metadata(::Val{:appointment_time}) =
    (column = :date_appointment, requires_infection = false)

Tracing, vaccination and immunity dates, including those of named doses, are shown whether or not the person was infected. This setting changes only the line list, not the simulation.

EpiBranch.contacts Function

julia

contacts(
    state::SimulationState;
    reference_date
) -> DataFrames.DataFrame

julia

contacts(state::SimulationState; reference_date=Date(2020, 1, 1))

Who exposed whom in a simulated outbreak: a DataFrame with one row per exposure of a person by an infectious case. Filter on infected to keep only the transmission tree.

A BranchingProcess records every exposure within the run, whether or not it led to infection. A run cut short by max_cases, max_time or a stopping rule leaves out the exposures that would have come after it stopped, a case at max_generations makes no contacts, and a contact timed before its infector's own infection is no exposure and is left out too. A model written against generate_offspring records every contact it creates, infected or not; it runs one generation at a time and finishes the generation in progress when it stops, so it can run past max_cases and list exposures after max_time. A model written against contacts_of lists an uninfected exposure only of a person it creates, not of one who already existed. The models of a fixed population (HomogeneousProcess, MixingProcess, NetworkProcess, RoutedNetwork, HouseholdProcess) record only who infected whom by default, so every row has infected = true: a household member or network neighbour who was exposed but not infected does not appear. Attach an ExposedContactRecorder as the ModelSpec's recorder to have such a pair listed too, with infected = false.

Columns: from (id of the case), to (id of the person exposed), infected (whether that exposure infected them), generation (generation of the person exposed), infection_time (day of exposure, which is the day of infection when infected is true) and date_infection (the same as a date counted from reference_date).

Examples

julia

state = simulate(ModelSpec(BranchingProcess(NegBin(2.5, 0.16), Gamma(2.0, 3.0)));
    max_cases = 100)
ct = contacts(state)
tree = ct[ct.infected, :]   # infector-infectee pairs

EpiBranch.chain_statistics Function

julia

chain_statistics(
    state::SimulationState
) -> DataFrames.DataFrame

julia

chain_statistics(state::SimulationState)

Size and length of each transmission chain in one simulated outbreak, one row per chain. A chain is all the cases descended from one index case. Only infected people are counted.

Returns a DataFrame with columns chain_id, size (number of cases in the chain) and length (number of generations of onward transmission: the highest generation reached, so an index case that infects nobody has size = 1 and length = 0).

Chain length is one less than in epichains

length here, as in ChainLengths, counts generations of transmission after the index case. epichains' chain_length counts generations including the index case, so a single-case chain has length 1 there and 0 here. Add 1 to compare with epichains.

Examples

julia

model = ModelSpec(BranchingProcess(NegBin(0.8, 0.5), Gamma(2.0, 3.0)))
state = simulate(model; n_initial = 20)
chain_statistics(state)   # one row per index case

julia

chain_statistics(
    states::Vector{<:SimulationState}
) -> DataFrames.DataFrame

julia

chain_statistics(states::Vector{<:SimulationState})

Chain sizes and lengths across several simulated outbreaks, such as the output of simulate(model, n). Returns a DataFrame with columns sim_id (which simulation), chain_id, size and length, defined as for a single outbreak above.

EpiBranch.compute_trace_level! Function

julia

compute_trace_level!(
    state::SimulationState
) -> SimulationState

julia

compute_trace_level!(state::SimulationState) -> state

Record how far along a chain of contact tracing each traced person was found: trace_level 0 for the case a round of tracing started from, 1 for its traced contacts, 2 for contacts of contacts, and so on. People never traced (and cases from which no tracing started) get no level. Call it after simulate; the ! means it changes state in place, adding the level to each person's record. linelist then shows it as a trace_level column.

Level along the last tracing path, not the shortest

Every case whose tracing reaches a person records itself as their tracer, so the level is counted along the path of the last one to reach them. In a BranchingProcess that is their infector, so the level is exact. On a NetworkProcess, where a person can have several infectious neighbours, it may be longer than the shortest path to an index case.

julia

compute_trace_level!(
    states::AbstractVector{<:SimulationState}
) -> AbstractVector{<:SimulationState}

julia

compute_trace_level!(states::AbstractVector{<:SimulationState}) -> states

Add trace levels to each of several simulated outbreaks, as compute_trace_level! does for one.

EpiBranch.realised_generation_interval Function

julia

realised_generation_interval(
    ind::Individual,
    state::SimulationState
) -> Any

julia

realised_generation_interval(ind::Individual, state::SimulationState)

The realised generation interval of one case: the days from its infector's infection to its own. Index cases and people who were never infected return NaN.

The interval is shaped by the epidemic. Susceptible depletion and interventions decide which transmissions happen, so realised intervals are distributed differently from the intrinsic generation_time you give the model.

EpiBranch.realised_generation_intervals Function

julia

realised_generation_intervals(
    state::SimulationState
) -> Vector{Float64}

julia

realised_generation_intervals(state::SimulationState)

The realised generation intervals (days) of every case in one simulation except the index cases, as a vector. See realised_generation_interval.

For a model with a fixed generation_time, the intrinsic interval is that distribution itself. For a state-dependent generation_time, run the model without interventions and the realised intervals coincide with the intrinsic ones, since nothing blocks transmission, as long as every case's contacts were decided. A run stopped by max_cases or max_time leaves out the late infections of its last cases, which shortens the realised intervals; bound it by max_generations instead.

Examples

julia

using Statistics
state = simulate(ModelSpec(BranchingProcess(Poisson(2.0), Gamma(2.0, 3.0)));
    n_initial = 5, max_cases = nothing, max_generations = 6)
mean(realised_generation_intervals(state))   # compare with mean(Gamma(2.0, 3.0)) = 6

julia

realised_generation_intervals(
    states::Vector{<:SimulationState}
) -> Vector{Float64}

julia

realised_generation_intervals(states::Vector{<:SimulationState})

The realised generation intervals (days) across several simulations, pooled into one vector.

EpiBranch.containment_probability Function

julia

containment_probability(
    states::Vector{<:SimulationState};
    max_cases
) -> Float64

julia

containment_probability(states::Vector{<:SimulationState}; max_cases=NoCases())

Containment probability: the proportion of simulated outbreaks that died out, from a vector of simulations such as simulate(model, n). An outbreak still going when the simulation stopped (at max_cases, max_time or max_generations) counts as not contained.

With max_cases, an outbreak that reached that many cases also counts as not contained, even if it then died out. Pass the cap used in simulate, so that outbreaks large enough to reach it count as uncontrolled, as in ringbp.

This is the simulation estimate. probability_contain gives the closed-form value for negative binomial offspring with simple control (ind_control, pop_control), and extinction_probability the probability that transmission dies out from a single case.

Examples

julia

model = ModelSpec(
    BranchingProcess(NegBin(2.5, 0.16), Gamma(2.0, 3.0));
    attributes = clinical_presentation(incubation_period = LogNormal(1.6, 0.5)),
    interventions = [Isolation(onset_to_isolation_delay = Exponential(2.0),
        duration = Inf)]
)
states = simulate(model, 1000; max_cases = 5000)
containment_probability(states; max_cases = 5000)

EpiBranch.is_extinct Function

julia

is_extinct(
    state::SimulationState;
    by_week,
    max_cases
) -> Bool

julia

is_extinct(state::SimulationState; by_week=nothing, max_cases=NoCases())

Whether one simulated outbreak died out.

With no keywords: true if the outbreak ended with no cases left to infect anyone before the simulation stopped.
by_week = 12: true if no case has symptom onset in week 12; by_week = 12:16: none in weeks 12 to 16.
max_cases: an outbreak that reached this many cases counts as not extinct.

Weeks are 7-day blocks counted from day 0 of the simulation (week 1 is days 0 to 6). A case without an onset time, such as an asymptomatic case, is placed by its infection time instead, as in weekly_incidence.

Examples

julia

using Statistics
model = ModelSpec(
    BranchingProcess(NegBin(1.2, 0.5), Gamma(2.0, 3.0));
    attributes = clinical_presentation(incubation_period = LogNormal(1.6, 0.5))
)
states = simulate(model, 100; max_time = 140.0)
mean(is_extinct.(states; by_week = 12:20))   # share with no onsets in weeks 12-20

EpiBranch.generation_R Function

julia

generation_R(state::SimulationState) -> DataFrames.DataFrame

julia

generation_R(state::SimulationState)

The ratio of the number of cases in each generation to the number in the generation before, in one simulated outbreak. Returns a DataFrame with columns generation (g) and offspring_ratio (cases in generation g + 1 divided by cases in generation g).

This differs from the time-varying reproduction number Rt estimated from an incidence time series; the two coincide only under strong assumptions. Use it to see how transmission falls from one generation to the next in a simulation, for example through depletion of susceptibles or interventions.

In a run stopped by max_cases or max_time, the cases of the last generations have not made all their contacts, so the ratios for those generations are right-truncated and fall towards 0. A run bounded by max_generations alone (with max_cases = nothing) gives complete ratios up to the bound.

EpiBranch.weekly_incidence Function

julia

weekly_incidence(
    state::SimulationState;
    by,
    reference_date
) -> DataFrames.DataFrame

julia

weekly_incidence(state::SimulationState; by=:onset,
                 reference_date::Date=Date(2020, 1, 1))

Weekly case counts from one simulated outbreak, the epidemic curve a surveillance system would plot. Returns a DataFrame with columns week (the date of the Monday starting each week, counted from reference_date as simulation day 0) and cases. Weeks with no cases are left out.

by chooses which date places a case in a week:

:onset (default): symptom onset, as on a surveillance epidemic curve. A case without an onset time (for example an asymptomatic case from clinical_presentation) is placed by its infection time.
:infection: infection time. This is not observed in real surveillance, and the curve comes out earlier by roughly one incubation period.
:reporting: reporting time, set by Reporting or by PerCaseObservation, whichever the model uses. Cases never reported are left out.
any other recorded time, such as :admission_time.

EpiBranch.scenario_sweep Function

julia

scenario_sweep(
    params::Dict{Symbol, <:AbstractVector};
    n_sim,
    rng,
    sim_kwargs...
) -> DataFrames.DataFrame

julia

scenario_sweep(params::Dict{Symbol, Vector}; n_sim=500, rng=Random.default_rng(), sim_kwargs...)

Containment probability for every combination of scenarios. Each row of the returned DataFrame holds one combination of the listed values and its containment_probability over n_sim simulated outbreaks of a BranchingProcess.

params maps each setting to the values to try. It must include :offspring (offspring distributions) and may include :generation_time, :interventions (each value is a vector of interventions applied together), :attributes and :population_size. Those five are the only axes the sweep varies, and any other key is rejected, including settings that do change a simulation, such as :progression. Settings that apply to every run, such as max_cases, are passed once as keywords. A setting that is a sweep axis is given as a vector even when there is only one value to use.

julia

results = scenario_sweep(
    Dict(
        :offspring => [NegBin(2.5, 0.16), NegBin(1.5, 0.5)],
        # Each value is a whole stack of interventions, so one vector per delay.
        :interventions => [
            [Isolation(onset_to_isolation_delay = Exponential(d), duration = 7.0)]
                for d in [1.0, 5.0]
        ],
        :attributes => [clinical_presentation(incubation_period = LogNormal(1.6, 0.5))],
        # Without a generation time every secondary case is infected at its
        # infector's own infection time, and an isolation that starts later
        # has nothing left to prevent.
        :generation_time => [LogNormal(1.6, 0.5)],
    );
    n_sim = 200, max_cases = 500
)

Four combinations, so four rows.

Closed-form results

Results calculated exactly from the offspring distribution, without simulation: the reproduction number, extinction and epidemic probabilities, superspreading summaries (as in the superspreading R package) and the distribution of chain sizes (as in epichains).

Reproduction number, extinction and superspreading

EpiBranch.reproduction_number Function

julia

reproduction_number(offspring::HouseholdOffspring)

R_, the household reproduction number: the mean number of other households infected by one infected household early in the epidemic. Transmission between households can grow only if R_ exceeds 1. It averages the mean of each household type over the probability that a community contact reaches that type.

Every household after the first is reached through a community contact with one of its members, so its type does not depend on the type of the household that infected it. The next-generation matrix over household types therefore has rank one, and its dominant eigenvalue is this single weighted sum.

julia

reproduction_number(
    d::Distributions.Distribution{Distributions.Univariate, Distributions.Discrete}
) -> Any

julia

reproduction_number(model)
reproduction_number(offspring)

Reproduction number implied by the model's offspring distribution: the mean number of secondary cases per case in a fully susceptible population. It is not the time-varying reproduction number Rt. Given a model, each of its interventions changes the offspring distribution first, through EpiBranch.analytic_offspring_effect, so the result is the basic reproduction number only when the model has none; an intervention without such an effect raises an error.

Single-type model: the mean of the offspring distribution.
ClusterMixed offspring: that mean averaged over the chains' mixing distribution.
Multi-type model built from a next-generation matrix: R*, the largest eigenvalue of the mean next-generation matrix, whose [i, j] entry is the expected number of type-i cases infected by one type-j case. That is the matrix the model draws from, which matches the matrix supplied only if each type's offspring distribution takes its column sum as its mean (see BranchingProcess).

For single-type and multi-type models, an outbreak can grow only if the reproduction number exceeds 1. The ClusterMixed average is not such a threshold: each chain grows or not according to its own reproduction number, so chains can take off even when the average is below 1. Use extinction_probability for the probability that transmission dies out.

When not every type can infect every other, R* above 1 means some group of types can sustain transmission, but an outbreak started by a type that cannot reach that group still dies out. extinction_probability gives the probability for each type of index case.

Examples

julia

M = [1.5 0.3;
     0.3 1.0]
reproduction_number(BranchingProcess(M, R -> NegBin(R, 0.5), Exponential(5.0)))
reproduction_number(BranchingProcess(NegBin(2.5, 0.16)))  # 2.5
reproduction_number(ClusterMixed(Poisson, Gamma(2.0, 0.6)))  # 1.2

julia

reproduction_number(
    o::GeneratingFunctionOffspring
) -> Float64

julia

reproduction_number(o::GeneratingFunctionOffspring)

The law's own mean.

EpiBranch.extinction_probability Function

julia

extinction_probability(offspring::HouseholdOffspring; tol = 1e-10, max_iter = 1000)

The probability that transmission between households, started by one infected household, dies out without a major epidemic: one entry per household type (see HouseholdOffspring).

Every household infected after the first is reached by a community contact, so its type is drawn from the mixing probabilities whatever the type of the household that infected it. The probability s that transmission from such a household dies out is therefore the same for all of them and solves s = Σₙ mixing[n] · Gₙ(s), where Gₙ is the probability generating function of the offspring distribution of type n. The answer for a first household of type n is Gₙ(s), which differs across types only through how many households that first one infects.

julia

extinction_probability(
    R::Real,
    k::Real;
    tol,
    max_iter
) -> Any

julia

extinction_probability(R::Real, k::Real; tol=1e-10, max_iter=1000)

Probability that transmission from a single introduced case dies out without a major outbreak, when each case infects a negative binomial number of others with mean R and dispersion k (NegBin(R, k); smaller k means more superspreading). It is 1 when R ≤ 1. Interventions are not included; for those, see probability_contain or simulate and use containment_probability.

The result is the smallest solution of q = G(q), where G is the probability generating function of the offspring distribution, found by iteration to tolerance tol. Convergence is slow when R is close to 1: a warning is given if max_iter iterations are reached.

Examples

julia

extinction_probability(2.5, 0.16)   # strong superspreading: about 0.79
extinction_probability(2.5, 100.0)  # little superspreading: about 0.11

julia

extinction_probability(
    d::Distributions.Poisson;
    tol,
    max_iter
) -> Any

julia

extinction_probability(d::Poisson; tol=1e-10, max_iter=1000)
extinction_probability(d::NegativeBinomial; tol=1e-10, max_iter=1000)

Probability that transmission from a single introduced case dies out without a major outbreak, for a Poisson or negative binomial offspring distribution d, such as NegBin(2.5, 0.16). Computed as for extinction_probability(R, k).

julia

extinction_probability(
    model::Union{TransmissionModel, ModelSpec};
    kwargs...
) -> Any

julia

extinction_probability(model::TransmissionModel; kwargs...)
extinction_probability(spec::ModelSpec; kwargs...)

Probability that transmission from a single introduced case dies out, from the model's offspring distribution. Population characteristics and the observation model are ignored. For a single-type model the result is one number; for a multi-type model built from an offspring matrix it has one value per type of index case.

For a ModelSpec, or a process that stores its own interventions, each intervention's EpiBranch.analytic_offspring_effect is applied to the offspring distribution first. If an intervention has no such method, this throws an error that names the simulation-based alternative.

julia

extinction_probability(
    o::EpiBranch.MultiTypeOffspring;
    tol,
    max_iter
) -> Any

julia

extinction_probability(o::MultiTypeOffspring; tol=1e-10, max_iter=1000)

Probability that transmission dies out in a multi-type branching process, one value per type: element j is the probability that an outbreak started by a single type-j case dies out. Call it on a model as extinction_probability(model).

It is 1 for types whose cases cannot lead, directly or through other types, to a group of types that sustain transmission among themselves (reproduction number above 1), and for every type when reproduction_number is at most 1. A warning is given if the calculation has not converged after max_iter iterations, which happens when the reproduction number is close to 1.

How it is calculated

The result is the smallest solution in [0, 1] of q_j = G_j(Σ_i a_ij q_i), where G_j is the probability generating function of dist_fn(R_j) and a_ij = M[i, j] / R_j is the share of a type-j case's secondary cases that are of type i. The rules above assume the number of secondary cases varies between cases; a fixed offspring count such as Dirac(1) at R = 1 gives chains that never end, for which returning 1 is wrong.

julia

extinction_probability(
    o::ClusterMixed;
    tol,
    max_iter
) -> Any

julia

extinction_probability(o::ClusterMixed; tol=1e-10, max_iter=1000)

Probability that a chain started by a single index case dies out when the reproduction number varies between chains: the extinction probability of build(θ), averaged over θ drawn from mixing. Chains whose θ gives a mean of at most 1 die out with certainty (assuming the number of secondary cases is not fixed).

mixing can be any continuous distribution or a DiscreteNonParametric distribution over a set of values.

julia

extinction_probability(
    o::GeneratingFunctionOffspring;
    tol,
    max_iter
) -> Any

julia

extinction_probability(o::GeneratingFunctionOffspring; tol=1e-10, max_iter=1000)

Probability that transmission from a single introduced case dies out, for an offspring law given by its generating function. Found by the same fixed-point iteration as extinction_probability(d::Poisson).

EpiBranch.epidemic_probability Function

julia

epidemic_probability(offspring::HouseholdOffspring; kwargs...)

The probability that one infected household of each type in offspring starts transmission between households that does not die out (one minus extinction_probability).

julia

epidemic_probability(R::Real, k::Real; kwargs...) -> Any

julia

epidemic_probability(R::Real, k::Real; kwargs...)

Probability that a single introduced case leads to a major epidemic, with NegBin(R, k) offspring: one minus extinction_probability.

julia

epidemic_probability(offspring; kwargs...) -> Any

julia

epidemic_probability(offspring; kwargs...)
epidemic_probability(model; kwargs...)

Probability that a single introduced case leads to a major epidemic: one minus extinction_probability, for an offspring distribution or a model. For a multi-type model there is one value per type of index case.

Given a model, population characteristics and the observation model are ignored. Each of its interventions changes the offspring distribution first, through EpiBranch.analytic_offspring_effect; an intervention without one raises an error. For control measures with no such effect, use probability_contain, or simulate and use containment_probability.

EpiBranch.probability_contain Function

julia

probability_contain(
    R::Real,
    k::Real;
    n_initial,
    ind_control,
    pop_control,
    tol,
    max_iter
) -> Any

julia

probability_contain(R, k; n_initial=1, ind_control=0.0, pop_control=0.0)

Closed-form probability that an outbreak is contained (dies out), with NegBin(R, k) offspring, under two simple kinds of control and with several introductions. It is the closed-form counterpart of the simulation estimate containment_probability, and reduces to extinction_probability with no control and one introduction.

ind_control: the probability that each case is controlled individually, for example isolated, before it infects anyone.
pop_control: the proportional reduction in R from population-wide measures such as social distancing; the effective R is (1 - pop_control) * R.
n_initial: the number of independent introductions.

For one introduction, the containment probability q solves

julia

q = ind_control + (1 - ind_control) * G(q)

where G is the probability generating function of the offspring distribution with the effective R. For n_initial introductions it is q^n_initial.

This is a port of probability_contain (and the probability_extinct equation it builds on) in the R package superspreading (Lambert et al., https://github.com/epiverse-trace/superspreading, MIT).

Examples

julia

# isolating half of cases before they transmit, three introductions
probability_contain(2.5, 0.16; ind_control = 0.5, n_initial = 3)

julia

probability_contain(
    d::Distributions.NegativeBinomial;
    kwargs...
) -> Any

julia

probability_contain(d::Poisson; n_initial=1, ind_control=0.0, pop_control=0.0)
probability_contain(d::NegativeBinomial; n_initial=1, ind_control=0.0, pop_control=0.0)

Closed-form containment probability for a Poisson or negative binomial offspring distribution, with the same keywords as probability_contain(R, k). Poisson offspring is treated as negative binomial with very large k.

julia

probability_contain(
    model::Union{TransmissionModel, ModelSpec};
    kwargs...
) -> Any

julia

probability_contain(model::TransmissionModel; kwargs...)
probability_contain(spec::ModelSpec; kwargs...)

Closed-form containment probability for a single-type model, from its offspring distribution, with the same keywords as probability_contain(R, k).

julia

probability_contain(
    o::GeneratingFunctionOffspring;
    n_initial,
    ind_control,
    pop_control,
    tol,
    max_iter
) -> Any

julia

probability_contain(o::GeneratingFunctionOffspring; n_initial=1, ind_control=0.0, pop_control=0.0)

Closed-form containment probability for an offspring law given by its generating function, with the same keywords as probability_contain(R, k). pop_control's proportional reduction in R is applied by thinning: keeping each secondary case with probability 1 - pop_control, which rescales the law's generating function the same way it rescales a negative binomial's or a Poisson's mean.

EpiBranch.end_of_outbreak_probability Function

julia

end_of_outbreak_probability(
    R::Real,
    k::Real,
    generation_time::Distributions.Distribution,
    τ::Real
) -> Any

julia

end_of_outbreak_probability(R, k, generation_time::Distribution, τ::Real)

Probability that a cluster is over, so no further cases will occur, given that its most recent case was τ days ago. Assumes NegBin(R, k) offspring, a generation time distribution generation_time (days) and that every case is reported. Following Thompson, Morgan & Jansen (2019, Phil Trans B), it is

julia

((k + R · (1 − S(τ))) / (k + R))^k

where S(τ) is the probability that a generation time exceeds τ. It rises from the probability of a case infecting nobody at τ = 0 towards 1 as τ grows.

Use it as prob_concluded in loglikelihood(::ChainSizes, ...) or chain_size_distribution for real-time cluster data, where some clusters may still be growing.

Examples

julia

# 14 days since the last case, generation time with mean 6 days
end_of_outbreak_probability(0.8, 0.5, Gamma(2.0, 3.0), 14.0)

julia

end_of_outbreak_probability(
    offspring::Distributions.NegativeBinomial,
    generation_time::Distributions.Distribution,
    τ::Real
) -> Any

julia

end_of_outbreak_probability(offspring::NegativeBinomial, generation_time, τ)
end_of_outbreak_probability(offspring::Poisson, generation_time, τ)

End-of-outbreak probability with the offspring distribution given directly, such as NegBin(0.8, 0.5) or Poisson(0.8). For Poisson offspring (no superspreading) it is exp(−R · S(τ)).

julia

end_of_outbreak_probability(
    model::Union{BranchingProcess, ModelSpec},
    τ::Real
) -> Any

julia

end_of_outbreak_probability(model::BranchingProcess, τ)
end_of_outbreak_probability(spec::ModelSpec, τ)

End-of-outbreak probability using the offspring and generation time distributions of a BranchingProcess, τ days after the most recent case.

This assumes every case is reported. Under-reporting (a Reporting step with a probability below 1, or any deprecated PerCaseObservation) needs the recursion of Thompson, Morgan & Jansen (2019), which is not implemented, and gives an error; use a model that reports every case for the full-reporting value.

julia

end_of_outbreak_probability(
    R::Real,
    k::Real,
    generation_time::Distributions.Distribution,
    τs::AbstractVector{<:Real}
) -> Any

julia

end_of_outbreak_probability(R, k, gt, τs::AbstractVector)

End-of-outbreak probability for each of several clusters, given the days since each cluster's most recent case. Returns one value per cluster, ready to pass as prob_concluded to loglikelihood(::ChainSizes, ...).

EpiBranch.proportion_transmission Function

julia

proportion_transmission(R::Real, k::Real; prop_cases) -> Any

julia

proportion_transmission(R::Real, k::Real; prop_cases::Real=0.2)

Proportion of all transmission caused by the most infectious fraction prop_cases of cases, with NegBin(R, k) offspring. This is the "80/20 rule" measure of superspreading: with prop_cases = 0.2, the share of transmission caused by the 20% most infectious cases.

Cases are ranked by their individual reproduction number (the expected number of people they infect, Gamma-distributed with mean R and shape k). The result depends only on k, so two calls with the same k and different R give the same value.

Examples

julia

proportion_transmission(2.5, 0.16)   # share of transmission from the top 20%

julia

proportion_transmission(
    d::Distributions.NegativeBinomial;
    prop_cases
) -> Any

julia

proportion_transmission(model::BranchingProcess; prop_cases=0.2)
proportion_transmission(spec::ModelSpec; prop_cases=0.2)

Proportion of transmission caused by the most infectious fraction prop_cases of cases, using the model's offspring distribution, which must be negative binomial (or Poisson). Population characteristics in a ModelSpec are ignored.

EpiBranch.proportion_cases_individual Function

julia

proportion_cases_individual(
    R::Real,
    k::Real;
    prop_transmission
) -> Any

julia

proportion_cases_individual(R::Real, k::Real; prop_transmission::Real=0.8)

Proportion of cases responsible for a share prop_transmission of all transmission (for example the 80% in the "80/20 rule"), with NegBin(R, k) offspring, ranking cases by their individual reproduction number. The inverse of proportion_transmission. It depends only on k.

Two ways to rank cases

proportion_cases_individual ranks cases by their individual reproduction number, a continuous Gamma-distributed quantity. proportion_cases_offspring ranks them by the number of people they actually infected, a whole number. The two answer slightly different questions and can differ substantially for the same R and k. When reporting, say which one was used.

julia

proportion_cases_individual(
    d::Distributions.NegativeBinomial;
    prop_transmission
) -> Any

julia

proportion_cases_individual(d::NegativeBinomial; prop_transmission=0.8)

Proportion of cases responsible for a share prop_transmission of transmission, ranking by individual reproduction number, for a negative binomial offspring distribution.

julia

proportion_cases_individual(
    model::Union{TransmissionModel, ModelSpec};
    prop_transmission
) -> Any

julia

proportion_cases_individual(model::BranchingProcess; prop_transmission=0.8)
proportion_cases_individual(spec::ModelSpec; prop_transmission=0.8)

Proportion of cases responsible for a share prop_transmission of transmission, ranking by individual reproduction number, using the model's offspring distribution (negative binomial or Poisson).

EpiBranch.proportion_cases_offspring Function

julia

proportion_cases_offspring(
    d::Distributions.Distribution{Distributions.Univariate, Distributions.Discrete};
    prop_transmission
) -> Any

julia

proportion_cases_offspring(d::DiscreteUnivariateDistribution; prop_transmission::Real=0.8)

Proportion of cases responsible for a share prop_transmission of all transmission, ranking cases by the number of people they actually infected. This answers "what share of infected people caused 80% of onward infections?". Works for any discrete offspring distribution d with a finite mean.

Cases are ranked from the most secondary cases downwards. Where the target share is crossed part-way through cases with the same count, those cases are counted fractionally, so the target share is met exactly.

See the note in proportion_cases_individual on how the two functions differ.

julia

proportion_cases_offspring(
    R::Real,
    k::Real;
    prop_transmission
) -> Any

julia

proportion_cases_offspring(R::Real, k::Real; prop_transmission::Real=0.8)

Proportion of cases responsible for a share prop_transmission of transmission, ranking by number of people infected, with NegBin(R, k) offspring.

julia

proportion_cases_offspring(
    model::Union{TransmissionModel, ModelSpec};
    prop_transmission
) -> Any

julia

proportion_cases_offspring(model::BranchingProcess; prop_transmission=0.8)
proportion_cases_offspring(spec::ModelSpec; prop_transmission=0.8)

Proportion of cases responsible for a share prop_transmission of transmission, ranking by number of people infected, using the model's offspring distribution.

EpiBranch.proportion_cluster_size Function

julia

proportion_cluster_size(
    R::Real,
    k::Real;
    cluster_size
) -> Any

julia

proportion_cluster_size(R, k; cluster_size=10)

Proportion of all secondary cases caused by cases who each infected at least cluster_size people (superspreading events), with NegBin(R, k) offspring. With strong superspreading (small k), a large share of cases come from a few such events.

Here cluster_size is the number of secondary cases of one infector, not the size of a transmission chain as in chain_size_distribution.

It is E[X; X ≥ c] / E[X] for the offspring distribution X and c = cluster_size.

julia

proportion_cluster_size(
    d::Distributions.NegativeBinomial;
    cluster_size
) -> Any

julia

proportion_cluster_size(d::NegativeBinomial; cluster_size=10)

Proportion of secondary cases caused by cases who each infected at least cluster_size people, for a negative binomial offspring distribution.

julia

proportion_cluster_size(
    model::Union{TransmissionModel, ModelSpec};
    cluster_size
) -> Any

julia

proportion_cluster_size(model::BranchingProcess; cluster_size=10)
proportion_cluster_size(spec::ModelSpec; cluster_size=10)

Proportion of secondary cases caused by cases who each infected at least cluster_size people, using the model's offspring distribution.

EpiBranch.heterogeneous_contact_R Function

julia

heterogeneous_contact_R(
    mean_contacts::Real,
    sd_contacts::Real,
    duration::Real,
    prob_transmission::Real
) -> NamedTuple{(:R, :R_net), <:Tuple{Any, Any}}

julia

heterogeneous_contact_R(mean_contacts, sd_contacts, duration, prob_transmission)

Basic reproduction number when people differ in how many contacts they have. People with many contacts are both more likely to be infected and to infect more others, which raises R above what the average number of contacts suggests.

mean_contacts, sd_contacts: mean and standard deviation of the number of contacts one person makes per unit of time.
duration: duration of infectiousness, in the same unit of time.
prob_transmission: probability of transmission per contact.

Returns (R = ..., R_net = ...):

R: assuming everyone has the average number of contacts, prob_transmission × mean_contacts × duration.
R_net: allowing for the variation in contacts, prob_transmission × duration × (mean + variance / mean).

This assumes contacts are formed at random given each person's number of contacts, with no clustering (friends of friends being friends). A simulation on an explicit network with clustering, such as NetworkProcess, can give a different answer. It is a port of calc_network_R in the R package superspreading (Lambert et al., https://github.com/epiverse-trace/superspreading, MIT).

Examples

julia

heterogeneous_contact_R(10.0, 15.0, 5.0, 0.01)

EpiBranch.analytic_offspring_effect Function

julia

analytic_offspring_effect(
    iv::AbstractIntervention,
    offspring
) -> Union{Distributions.NegativeBinomial, Distributions.Poisson}

julia

analytic_offspring_effect(intervention, offspring) -> offspring
analytic_offspring_effect(intervention, offspring, spec) -> offspring

Return the offspring distribution a case draws from under intervention, given the one a bare process would use (a Distribution for a single-type process, a MultiTypeOffspring for a multi-type one). Define a method only when the effect is exact in closed form on its own, and return a law the analytical functions take: a Poisson, a NegativeBinomial, or a GeneratingFunctionOffspring built from the offspring's generating function or probability mass function. A population-level reduction of R by a factor p, for example, turns NegBin(R, k) into NegBin(p * R, k), as MassVaccination's method does for a campaign already complete before the simulation starts.

The three-argument form also receives the ModelSpec (or process that stores its own interventions), for an effect that needs more than the offspring law itself, such as the generation-time distribution or the incubation period. Its default falls back to the two-argument method, so an intervention whose effect does not depend on the rest of the model need only define that one; Isolation defines the three-argument form, since its effect on R depends on how its timing overlaps the generation time.

An effect with no closed form at all, such as contact tracing's or ring and group vaccination's, has no method of either arity. The generic method throws an error, which makes extinction_probability, epidemic_probability, probability_contain, reproduction_number, proportion_transmission, proportion_cases_individual, proportion_cases_offspring, proportion_cluster_size, end_of_outbreak_probability and offspring_distribution throw too when called on a ModelSpec, or a process that stores its own interventions, holding an intervention that has no method.

julia

analytic_offspring_effect(
    iso::Isolation,
    offspring::Distributions.Distribution{Distributions.Univariate, Distributions.Discrete},
    spec
) -> GeneratingFunctionOffspring{G} where G<:(EpiBranch.var"#pgf#analytic_offspring_effect##6"{var"#s185", EpiBranch.var"#isolated_pgf#analytic_offspring_effect##1"{offspring, keep_probability, delay, incubation}} where {var"#s185"<:Distributions.Distribution{Distributions.Univariate, Distributions.Discrete}, offspring, keep_probability, delay, incubation})

julia

analytic_offspring_effect(iso::Isolation, offspring::DiscreteUnivariateDistribution, spec)

The offspring law under isolation: a case drawn from offspring, thinned by the probability that each of its secondary cases still happens (see the comment above this method for the formula), averaged over spec's incubation period and this isolation's own onset_to_isolation_delay, and over whether the case is isolated at all (eligibility and test_sensitivity). Returns a GeneratingFunctionOffspring.

Needs spec's attributes to include clinical_presentation and its process to have a generation-time distribution; duration must be a fixed number (or Inf), onset_to_isolation_delay and test_sensitivity each a number or a distribution, and eligibility SymptomaticOnly or AllCases. Anything else has no closed form here, and the generic default throws, pointing at simulation instead.

EpiBranch.GeneratingFunctionOffspring Type

julia

struct GeneratingFunctionOffspring{G}

julia

GeneratingFunctionOffspring(pgf, mean)
GeneratingFunctionOffspring(pmf::AbstractVector{<:Real})

An offspring law given directly by its probability generating function pgf(s), or built from a probability mass function pmf over 0, 1, 2, ... (pmf[i] is the probability of i - 1). mean is the law's reproduction number, the mean number of secondary cases.

Return one of these from EpiBranch.analytic_offspring_effect for an effect with no simpler closed form, such as a mixture over a random isolation time: Isolation does this. extinction_probability and probability_contain accept it directly, and so does reproduction_number. proportion_transmission, proportion_cases_individual, proportion_cases_offspring and end_of_outbreak_probability need a negative binomial or Poisson's own dispersion or rate. A law known only through its generating function has none, so these raise an ArgumentError pointing at simulation instead.

Examples

julia

# Poisson(1.5) thinned by keeping each secondary case with probability 0.6
g = GeneratingFunctionOffspring(s -> exp(1.5 * (0.6 * s + 0.4 - 1)), 1.5 * 0.6)
extinction_probability(g)

# the same law from its probability mass function
using Distributions
pmf = [pdf(Poisson(0.9), n) for n in 0:200]
GeneratingFunctionOffspring(pmf)

Fields

pgf::Any
mean::Float64

Chain-size distributions

chain_size_distribution gives the probability of each final chain size (the total number of cases descending from one introduction), calculated exactly from the offspring distribution. Borel is the chain-size distribution when each case infects a Poisson number of others.

EpiBranch.chain_size_distribution Function

julia

chain_size_distribution(d::Distributions.Poisson) -> Borel

julia

chain_size_distribution(offspring::Poisson)

Exact chain size distribution, the probability that a chain started by one case has n cases in total, when each case infects a Poisson number of others: a Borel distribution. When R > 1 the probabilities of finite sizes add up to less than 1, the rest being the probability of a major outbreak.

Examples

julia

d = chain_size_distribution(Poisson(0.8))
pdf(d, 3)   # probability a chain has exactly 3 cases

julia

chain_size_distribution(
    d::Distributions.NegativeBinomial
) -> EpiBranch.GammaBorel

julia

chain_size_distribution(offspring::NegativeBinomial)

Exact chain size distribution, the probability that a chain started by one case has n cases in total, when each case infects a negative binomial number of others (a GammaBorel distribution). When R > 1 the probabilities of finite sizes add up to less than 1, the rest being the probability of a major outbreak.

Examples

julia

d = chain_size_distribution(NegBin(0.8, 0.5))
pdf(d, 3)                   # probability a chain has exactly 3 cases
sum(pdf(d, n) for n in 1:10)  # probability a chain has at most 10 cases

julia

chain_size_distribution(
    o::ClusterMixed
) -> EpiBranch.PoissonGammaChainSize

julia

chain_size_distribution(o::ClusterMixed)

Chain size distribution when the reproduction number varies between chains: exact for Poisson offspring with a Gamma-distributed rate (PoissonGammaChainSize), otherwise computed numerically (ChainSizeMixture).

julia

chain_size_distribution(
    model::TransmissionModel;
    seeds,
    prob_concluded,
    kwargs...
) -> Any

julia

chain_size_distribution(model; seeds=nothing, prob_concluded=nothing, kwargs...)
chain_size_distribution(spec::ModelSpec; seeds=nothing, prob_concluded=nothing, kwargs...)

Distribution of observed chain sizes (the number of cases in a cluster) under model, including its interventions and observation model. Use it to fit a model to cluster-size data, for example in Turing.jl, where sizes ~ ... says the observed sizes follow this distribution:

julia

@model function chain_model(sizes)
    R ~ Gamma(2, 1)
    sizes ~ chain_size_distribution(BranchingProcess(NegBin(R, 0.5)))
end

Without interventions, and with a single-type offspring distribution, this is the exact chain-size distribution (Borel or GammaBorel, adjusted for the observation model). It stays exact with prob_concluded, and with seeds unless the observation model also under-reports: there is no closed form for several index cases under under-reporting, so that one combination is estimated from simulated outbreaks and its value moves from one call to the next. With interventions, or with a multi-type or multi-window model (no single-type offspring law to fall back on), the probability of the data is estimated from simulated outbreaks instead. Keyword arguments:

seeds: per-cluster number of index cases, for multi-seed clusters. With interventions, clusters are simulated grouped by seed count, one batch of runs per distinct value in seeds.
prob_concluded: per-cluster probability the cluster is finished (its observed size is its final size), for real-time data with ongoing clusters — see end_of_outbreak_probability. With interventions, the mixture is taken against the simulated clusters' empirical P(X ≥ x) rather than an analytical tail.
n_sim, …: forwarded to the underlying simulation-based loglikelihood when the analytical fast path does not apply.

EpiBranch.Borel Type

julia

struct Borel{T<:Real} <: ChainSizeLaw

julia

Borel(μ)

Chain size distribution when each case infects a Poisson number of others with mean μ: the probability that a chain started by one case has n cases in total,

julia

P(X = n) = (μn)^(n-1) * exp(-μn) / n!  for n = 1, 2, ...

chain_size_distribution(Poisson(μ)) returns it.

When μ > 1 some chains grow without end, so the probabilities of the finite sizes add up to less than 1; the missing probability is that of a major outbreak. The formula stays valid for each n, which lets chain-size probabilities be averaged over values of μ either side of 1.

Fields

μ::Real

Other chain-size distributions

These are returned by chain_size_distribution and are rarely built by hand. GammaBorel is the chain-size distribution for negative binomial offspring. PoissonGammaChainSize is the one for Poisson offspring whose mean varies between chains following a gamma distribution.

EpiBranch.GammaBorel Type

julia

struct GammaBorel{T<:Real} <: ChainSizeLaw

julia

GammaBorel(k, R)

Chain size distribution when each case infects a negative binomial number of others with dispersion k and mean R: the probability that a chain started by one case has n cases in total. chain_size_distribution(NegBin(R, k)) returns it. It matches .nbinom_size_ll in epichains.

Argument order is (k, R)

GammaBorel(k, R) takes the dispersion first, the reverse of NegBin(R, k). GammaBorel(0.5, 0.8) is the chain size distribution for NegBin(0.8, 0.5). Building it through chain_size_distribution(NegBin(R, k)) avoids the mix-up.

When R > 1 some chains grow without end, so the probabilities of the finite sizes add up to less than 1; the missing probability is that of a major outbreak.

R = 0 is allowed, and gives the degenerate case where nobody infects anyone else and every chain has a single case. Fitting over log k needs it: once k is large enough to swamp R in floating point, mean(NegativeBinomial) rounds to exactly zero, and the fit must still get a distribution back rather than an error.

Fields

k::Real
R::Real

EpiBranch.PoissonGammaChainSize Type

julia

struct PoissonGammaChainSize{T<:Real} <: ChainSizeLaw

julia

PoissonGammaChainSize(k, R)

Chain size distribution when the reproduction number varies between chains rather than between cases: each chain draws its own rate λ from a Gamma distribution with shape k and mean R, and every case in that chain infects a Poisson(λ) number of others. It matches the gborel likelihood in epichains. Like GammaBorel, it takes k first.

This differs from GammaBorel, where each case has its own reproduction number (negative binomial offspring). chain_size_distribution returns it for ClusterMixed(Poisson, Gamma(...)).

Fields

k::Real
R::Real

A different offspring distribution for the index case

A chain whose index case has its own offspring distribution, for example a chain started by an imported case, uses IndexChainSize.

EpiBranch.IndexChainSize Type

julia

struct IndexChainSize{J<:Distributions.Distribution{Distributions.Univariate, Distributions.Discrete}, S<:Distributions.Distribution{Distributions.Univariate, Distributions.Discrete}} <: ChainSizeLaw

julia

IndexChainSize(index_offspring, offspring)

Chain size distribution when the index case has a different offspring distribution from later cases: the index case infects a number drawn from index_offspring, and every later case a number drawn from offspring. Use it when an introduced case transmits less than locally infected cases, for example because they arrive part-way through their infectious period or are quarantined on arrival.

Only chains started by a single index case are supported (seeds == 1 in ChainSizes).

Note

The probability of a chain of n cases sums over the number j the index case infected, each starting its own chain:

P(N = n) = P(J = 0) 1{n = 1} + Σ_{j ≥ 1} P(J = j) P(chains from j cases have n - 1 cases)

This needs the chain size distribution of offspring from several starting cases, which exists for Poisson, NegativeBinomial and ClusterMixed(Poisson, ::Gamma). For other ClusterMixed offspring, pdf and logpdf give an error once n is large enough to need it.

Examples

julia

d = IndexChainSize(Poisson(0.3), NegBin(0.8, 0.5))
pdf(d, 5)
loglikelihood(ChainSizes([1, 2, 5, 1]), d)

Fields

index_offspring::Distributions.Distribution{Distributions.Univariate, Distributions.Discrete}
dist::Distributions.Distribution{Distributions.Univariate, Distributions.Discrete}

Household reproduction number

The household reproduction number R* is the expected number of other households infected by one infected household. These functions give R*, the distribution of the number of households one household infects, and the final size of the epidemic within a single household.

EpiHouseholds.household_offspring Function

julia

household_offspring(model; global_rate, n_samples = 10_000, rng, tol, max_offspring)

Spread between households in the early phase of an epidemic: how many other households one infected household goes on to infect. From the result, reproduction_number gives the household reproduction number R* (the epidemic can grow between households only if R* > 1) and extinction_probability the chance that transmission between households dies out. model is a HouseholdProcess, or a ModelSpec with one.

global_rate is the rate (per day) at which an infectious person makes contacts outside their household. Early in an epidemic each such contact reaches a susceptible person in a household not yet infected, so the number of households one infected household infects is Poisson with mean global_rate times the total time its members spend infectious. The result is a HouseholdOffspring with that distribution for each household type, and the probability that a community contact reaches each type (proportional to the number of people in households of that type).

The household sizes come from the model, as do the natural history and interventions: isolation shortens each case's infectious period in the community and lowers R* as it lowers transmission in a simulation. Each newly infected household starts with one index case, a member chosen at random, so the process must have no external_hazard.

Example

julia

using EpiBranch, EpiHouseholds, Distributions

model = ModelSpec(HouseholdProcess(fill(4, 1000), Exponential(1 / 0.5));
    progression = [Transition(:recovered; from = :infection, rate = 1 / 4,
        terminal = true)])

offspring = household_offspring(model; global_rate = 0.15)
reproduction_number(offspring)        # R*
extinction_probability(offspring)     # one per household type
household_offspring_law(offspring)    # distribution for a household reached by a community contact

How it is calculated

The epidemic within each household is calculated exactly where a closed form exists, and by simulating households otherwise. An exponential contact interval and an exponential infectious period, with no interventions, give an exact result. With a shared kernel, n_samples households of each size are simulated. A covariate kernel (infector, susceptible) -> Distribution has no closed form: the whole model is simulated as many times as it takes to reach n_samples households (at least once), so each household's epidemic follows the covariates of its actual members. Households whose kernels agree pair for pair, in member order, share a type; the rest are types of their own. When every household is its own type, as with a continuous individual covariate, each type's distribution comes from those few runs and is rough, while the overall distribution and R* average over all of them. The step from household epidemic to number of households infected is exact either way, so simulation error enters only through the within-household epidemic.

With a shared kernel and no interventions the mean is exact whenever the infectious period is a single delay in the progression, because the mean total time infectious is then the mean final size times the mean infectious period. An intervention cuts the window at a time that depends on the rest of the case's timeline, so with one in the model the mean is simulated. The exception is a large household with weak transmission and a random infectious period, where the exact calculation loses accuracy and the mean comes from the simulated households instead.

With a shared kernel the households of each size are simulated with interchangeable members, so the natural history and interventions must treat all members alike. A progression or attributes that read an individual's covariates need the kernel given as a function, even a constant one such as (i, j) -> Exponential(3.0).

Each household is simulated on its own timeline, so interventions cannot be wrapped in Scheduled. To see what switching a policy on does, calculate this twice, once without the intervention and once with it (unwrapped): the two R* values are the reproduction numbers before and after the switch.

tol bounds the probability left in the tail beyond the returned range of household counts, and max_offspring caps that range.

EpiHouseholds.HouseholdOffspring Type

The spread of infection between households, treating each infected household as one "case": the distribution of how many other households one infected household infects, from which the household reproduction number R* and the probability that transmission between households dies out follow. Returned by household_offspring.

Households are grouped into types that share the same within-household epidemic. With one contact-interval distribution shared by every pair, the type is the household size. With a kernel that varies by pair, households are the same type only when their kernels agree pair for pair; otherwise each is a type of its own. The fields hold one entry per type, in increasing size order:

sizes: the household size of each type.
households: the ids of the model's households of each type.
mixing: the probability that a community contact reaches a household of each type. A contact reaches a person, and with them their household, so this is proportional to the number of people living in households of that type.
laws: for each type, the distribution of the number of households it infects, as a DiscreteNonParametric.
means: the mean of each of those distributions, exact where the within-household calculation is, and taken from the simulated households otherwise (see household_offspring).
global_rate: the community contact rate used.

reproduction_number gives R*, extinction_probability the probability that transmission started by one infected household of each type dies out, and household_offspring_law the offspring distribution itself as a Distributions.jl distribution.

EpiHouseholds.household_offspring_law Function

julia

household_offspring_law(offspring::HouseholdOffspring)
household_offspring_law(offspring::HouseholdOffspring, size)

The distribution of the number of other households infected by one infected household, as a Distributions.jl distribution: for a household of the given size, or, with no size, for a household reached through a community contact, which is the distribution every household after the first follows. When several types share a size, as they can with a covariate kernel, the distribution for that size mixes them, weighted by how likely a contact is to reach each.

The distribution without a size treats households as the cases of a single-type branching process, so it can be passed to BranchingProcess to simulate chains of infected households.

EpiHouseholds.household_final_size Function

julia

household_final_size(size, kernel, window; initial_infectives = 1)

The final-size distribution of the epidemic within a single household of size members: how many of them are ultimately infected, counting the initial_infectives it starts with.

kernel is the within-household contact-interval distribution, the time from a case becoming infectious to it making infectious contact with a given household-mate, as in HouseholdProcess. window is the length of a case's infectious period, a distribution or a fixed number. Contact happens only if it falls inside the window, and contact with someone already infected is wasted.

The distribution comes from Ball's (1986) triangular recursion, which is exact for any kernel and window. It runs on the probability that a specified set of m susceptibles all escape one case, E[S(L)ᵐ] for a window length L and a kernel survival function S: contacts from one case are independent given its window, so this one quantity is all the final size needs from the kernel.

The recursion subtracts terms far larger than the probabilities they leave, so in a large household with weak transmission it can lose accuracy, sooner for a random window (whose escape probabilities are integrated numerically) than for a fixed one; the function then throws.

Examples

julia

using EpiHouseholds, Distributions

# a household of four, exponential contact intervals, a fixed six-day window
d = household_final_size(4, Exponential(8.0), 6.0)
mean(d)     # mean number infected, the index case included
pdf(d, 4)   # probability the whole household is infected

Inference

Fit models to outbreak data. Offspring counts, chain sizes and chain lengths are fitted by likelihood; individual infection times from households or contact networks are fitted by a pairwise survival likelihood.

Data types

EpiBranch.OffspringCounts Type

julia

struct OffspringCounts

julia

OffspringCounts(data)

Observed numbers of secondary cases: how many people each case infected, one count per case. Pass it to loglikelihood with an offspring distribution, such as NegBin(R, k), to estimate R and the dispersion k.

Examples

julia

data = OffspringCounts([0, 1, 2, 0, 3, 1, 0, 0, 5, 0])
loglikelihood(data, NegBin(0.8, 0.5))

# maximum-likelihood estimate of R for fixed k = 0.5, over a grid
Rs = 0.05:0.05:3.0
R_hat = Rs[argmax([loglikelihood(data, NegBin(R, 0.5)) for R in Rs])]

For R and k together, maximise loglikelihood with an optimiser such as Optim.jl, or fit in Turing.jl with offspring_distribution.

Fields

data::Vector{Int64}

EpiBranch.ChainSizes Type

julia

struct ChainSizes

julia

ChainSizes(data; seeds = ones(Int, length(data)))

Observed chain sizes: the total number of cases in each transmission chain (cluster). Pass it to loglikelihood with an offspring distribution or a model, or fit in Turing.jl with chain_size_distribution.

data: the size of each cluster.
seeds: the number of index cases (separate introductions) in each cluster, 1 by default.

Every cluster is assumed to be over, so its size is its final size. For real-time data where some clusters may still grow, pass prob_concluded (the probability each cluster is over) to loglikelihood; see end_of_outbreak_probability. For data recorded only once a cluster reaches a given size (for example only clusters of two or more cases), use a MinimumSize observation model.

Examples

julia

# one index case per cluster
data = ChainSizes([1, 1, 3, 1, 5])
loglikelihood(data, NegBin(0.8, 0.5))

# clusters with several introductions
data = ChainSizes([3, 5, 10, 2]; seeds = [1, 2, 1, 1])

# only clusters of two or more cases are recorded
data = ChainSizes([2, 3, 5, 2])
loglikelihood(data, observe(chain_size_distribution(NegBin(0.8, 0.5)), MinimumSize(2)))

Fields

data::Vector{Int64}
seeds::Vector{Int64}

EpiBranch.ChainLengths Type

julia

struct ChainLengths

julia

ChainLengths(data)

Observed chain lengths: the number of generations of onward transmission in each chain, 0 for a chain of a single case. Pass it to loglikelihood with an offspring distribution or a model, or fit in Turing.jl with chain_length_distribution.

Chain length is one less than in epichains

epichains' chain_length counts generations including the index case, so a single-case chain has length 1 there and 0 here. Subtract 1 from epichains chain lengths before using them here.

Examples

julia

data = ChainLengths([0, 1, 0, 2, 1])
loglikelihood(data, Poisson(0.5))

Fields

data::Vector{Int64}

EpiBranch.IncidenceCounts Type

julia

struct IncidenceCounts{D<:VecOrMat{Int64}, S<:Union{Int64, Vector{Int64}}}

julia

IncidenceCounts(data; seeds = 1)

Observed incidence: the number of new cases in each successive time step (day, week, or whatever unit data is recorded in) of a single outbreak, starting one step after the seeds index cases. Pass it to loglikelihood with an offspring distribution and a generation-time distribution, or with a model that already has its own, to fit a branching process to a case or admission time series without simulating it.

data: the incidence count for each step, in calendar order. For a multi-type process (MultiTypeOffspring), a vector holds the total over types, and a matrix with one row per step and one column per type holds the counts by type.
seeds: the number of index cases introduced at step 0, 1 by default. For a multi-type process, a vector with the number of each type.

Examples

julia

data = IncidenceCounts([1, 2, 4, 5, 9, 12, 18])
loglikelihood(data, NegBin(1.3, 0.5), LogNormal(1.6, 0.5))

# children and adults, one column each, started by one adult
by_type = IncidenceCounts([0 1; 1 1; 1 3; 2 4]; seeds = [0, 1])

Fields

data::VecOrMat{Int64}
seeds::Union{Int64, Vector{Int64}}

Likelihood and fitting

loglikelihood works on each kind of outbreak data, so parameters can be estimated by maximum likelihood or in a Bayesian model:

julia

loglikelihood(OffspringCounts(data), Poisson(0.5))
loglikelihood(OffspringCounts(data), NegBin.(μ, 0.5))  # one distribution per case
loglikelihood(ChainSizes(data), NegBin(0.8, 0.5))
loglikelihood(ChainLengths(data), Poisson(0.5))
loglikelihood(ChainSizes(data), model)   # interventions/observation read from model
loglikelihood(IncidenceCounts(data), NegBin(0.8, 0.5), LogNormal(1.6, 0.5))

The dot in NegBin.(μ, 0.5) applies NegBin to each element of μ, like vectorised R code. With one distribution per case, each count is evaluated against its own distribution, which allows case-level covariates such as NegBin.(exp.(X * β), k). For data that list only cases with at least one secondary case, pass the zero-truncated distribution truncated.(offspring, 1, Inf); Distributions.truncated works with a single distribution or with one per case.

For maximum-likelihood estimation, maximise loglikelihood with an optimiser such as Optim.jl, or use Turing's maximum_likelihood; the same Turing model containing data ~ chain_size_distribution(model) works for both. See the chains tutorial for examples.

Distributions for Bayesian fitting

These return the distribution of the data under a model, which can be written on the right-hand side of ~ in a Turing.jl model (data ~ ...). Where no exact formula exists, the distribution is estimated by simulation. The chain_size_distribution entry above lists its exact methods.

EpiBranch.chain_length_distribution Function

julia

chain_length_distribution(
    model::TransmissionModel;
    kwargs...
) -> EpiBranch._ChainLengthLaw{M, @NamedTuple{}} where M<:TransmissionModel

julia

chain_length_distribution(model; kwargs...)

Distribution of observed chain lengths (generations of onward transmission in a cluster, 0 for a single case) under model. Use it to fit a model to chain-length data, for example in Turing.jl:

julia

@model function chain_model(lengths)
    R ~ Beta(1, 1)
    lengths ~ chain_length_distribution(BranchingProcess(Poisson(R)))
end

Keywords n_sim, n_initial, max_cases, max_generations, max_time, stopping_rules and rng control the simulation used when no exact formula exists. Interventions come from the model; add them with a ModelSpec. See ChainLengths for the length convention.

EpiBranch.offspring_distribution Function

julia

offspring_distribution(
    model::TransmissionModel
) -> Distributions.Distribution

julia

offspring_distribution(model)
offspring_distribution(spec::ModelSpec)

Offspring distribution of model: the distribution of the number of secondary cases per case. For a BranchingProcess it is the distribution you gave as offspring. Use it to fit a model to counts of secondary cases, for example in Turing.jl:

julia

@model function offspring_model(counts)
    R ~ Beta(1, 1)
    counts ~ offspring_distribution(BranchingProcess(Poisson(R)))
end

Incidence time-series likelihood

A closed-form, moment-matched likelihood for case or admission counts per step (day, week, ...), built from the offspring and generation-time distributions rather than simulated. See its "Limits" note for what it ignores.

EpiBranch.incidence_moments Function

julia

incidence_moments(
    offspring::Union{GeneratingFunctionOffspring, Distributions.DiscreteUnivariateDistribution, Distributions.UnivariateDistribution{Distributions.Discrete}, Distributions.DiscreteDistribution{Distributions.Univariate}},
    generation_time::Distributions.Distribution,
    n::Integer;
    seeds,
    reporting_probability,
    reporting_delay,
    tol,
    max_lag
) -> NamedTuple{(:mean, :cov), <:Tuple{Any, LinearAlgebra.Symmetric}}

julia

incidence_moments(offspring, generation_time::Distribution, n::Integer; seeds = 1,
    reporting_probability = 1, reporting_delay = nothing)

Mean vector and covariance matrix of the incidence counts a branching process produces at steps 1:n, started by seeds index cases at step 0, with offspring law offspring (a Poisson, NegativeBinomial, GeneratingFunctionOffspring, or any other discrete distribution EpiBranch knows the generating function of) and a generation_time distribution giving the time from one case's infection to each of its secondary infections, in the same unit as a step of n (so a weekly series needs a generation time in weeks).

With reporting_probability below 1, each case is reported independently with that probability, and the counts are of reported cases. A reporting_delay (a distribution or a fixed number of steps, measured from infection, in the same unit) dates each report that many steps after the infection, so the counts are by report date; the seeds' own reports count too once a delay moves them into an observed step. Both default to every case counted in the step it is infected. The result is exact for independent reporting, and reporting_probability can be differentiated like the offspring parameters.

Returns a (mean, cov) named tuple: mean[t] is the expected count at step t, and cov is their n × n covariance matrix, both ignoring susceptible depletion (see loglikelihood for the likelihood this feeds and its limits). Matching a MvNormal(mean, cov) or a moment-matched lognormal to these is what the NormalApproximation and LognormalApproximation of that likelihood do; call this directly for a forecast instead, or to compare against simulated incidence (seeds simulated introductions, summed into steps the same way data is).

Examples

julia

m = incidence_moments(NegBin(1.3, 0.5), LogNormal(1.6, 0.5), 20; seeds = 3)
m.mean[end]
sqrt(m.cov[end, end])  # standard deviation of the last step's count

julia

incidence_moments(
    offspring::EpiBranch.MultiTypeOffspring,
    generation_time::Distributions.Distribution,
    n::Integer;
    seeds,
    reporting_probability,
    reporting_delay,
    tol,
    max_lag
)

julia

incidence_moments(offspring::MultiTypeOffspring, generation_time::Distribution,
    n::Integer; seeds, reporting_probability = 1, reporting_delay = nothing)

Mean and covariance of the incidence counts of each type of a multi-type branching process at steps 1:n, started by seeds[j] index cases of type j at step 0. generation_time is shared by every type, and reporting_probability and reporting_delay act as in the single-type method.

Returns a (mean, cov) named tuple: mean[t, j] is the expected count of type j at step t, and cov is the covariance matrix of vec(mean)'s entries, steps first and then types, so cov[t + (j - 1) * n, t′ + (j′ - 1) * n] is the covariance between type j at step t and type j′ at step t′.

Examples

julia

offspring = EpiBranch.MultiTypeOffspring([1.2 0.4; 0.3 0.9], R -> NegBin(R, 0.5))
m = incidence_moments(offspring, LogNormal(1.6, 0.5), 20; seeds = [0, 1])
m.mean[end, :]  # expected count of each type on the last step

EpiBranch.IncidenceApproximation Type

julia

abstract type IncidenceApproximation

julia

IncidenceApproximation

How loglikelihood(::IncidenceCounts, ::Distribution) turns the mean and covariance from incidence_moments into a likelihood: NormalApproximation matches a multivariate normal, and LognormalApproximation a multivariate lognormal. Dispatch on a new subtype to match something else.

Fields

EpiBranch.NormalApproximation Type

julia

struct NormalApproximation <: IncidenceApproximation

julia

NormalApproximation()

Match a MvNormal with the mean and covariance from incidence_moments. Defined for any observed counts, including zeros; see LognormalApproximation for the alternative that respects non-negativity.

Fields

EpiBranch.LognormalApproximation Type

julia

struct LognormalApproximation <: IncidenceApproximation

julia

LognormalApproximation()

Match a multivariate lognormal with the mean and covariance from incidence_moments: the log-counts are matched to a MvNormal by the usual lognormal moment equations, logpdf adding the log-Jacobian -Σ log(data) for the change of variables. Needs every observed count to be strictly positive, since a lognormal puts no mass at 0; use NormalApproximation for data that includes zeros.

Fields

Observation models

An observation model describes how outbreaks are detected or reported: for example each case is reported with some probability, or only chains of at least a given size are seen. Add one to a model with observation = ....

EpiBranch.ObservationModel Type

julia

abstract type ObservationModel

How true cases become observed data: which cases are detected, how long reporting takes, or which outbreaks are recorded at all. Pass one to a model as observation = .... Built in: NoObservation (every case seen), MinimumSize (only chains of at least a given size are recorded) and the deprecated PerCaseObservation (each case detected with some probability, after a reporting delay; use a Reporting step instead).

To write a new one, define observe (for the likelihood) and apply_observation! (for simulation) for it; see "Adding an observation model" in the Extending EpiBranch guide.

Fields

EpiBranch.NoObservation Type

julia

struct NoObservation <: ObservationModel

No observation model: latent cases are observed exactly.

Fields

EpiBranch.PerCaseObservation Type

julia

struct PerCaseObservation{P, D, F} <: ObservationModel

julia

PerCaseObservation(; detection_prob = 1.0, delay = Dirac(0.0),
                   from = :onset_time)

Deprecated: report cases with a Reporting step in the model's progression instead. simulate, progression_loglik, the chain-size likelihood and end_of_outbreak_probability all read the reporting probability from that step.

Under-reporting and reporting delay: each case is reported independently with probability detection_prob, delay days after symptom onset (or after from). Each simulated case records :reported and :reporting_time.

detection_prob: the probability a case is reported, a number, a distribution, or a function of the random number generator and the individual, (rng, ind) -> ... (for example to make reporting depend on age). It must be above 0 and at most 1; leave the observation model out rather than setting it to 0.
delay: days from from to report, as a number, distribution or such a function.
from: the time reporting is measured from. Symptom onset by default, because surveillance follows onset; from = ind -> ind.infection_time measures from infection. A case without an onset time (for example an asymptomatic case from clinical_presentation) is measured from its infection time.

detection_prob = 1.0, delay = Dirac(0.0) (the defaults) means every case is seen at once; detection_prob = ρ, delay = Dirac(0.0) is under-reporting alone.

Note

The closed-form chain-size results (observe, ThinnedChainSize) need detection_prob to be a single number and give an error otherwise. For reporting that varies between cases, use the simulation-based likelihood.

Examples

julia

# 60% of cases reported, on average 3 days after onset
PerCaseObservation(detection_prob = 0.6, delay = Gamma(3.0, 1.0))

Report cases with this or with a Reporting step in the progression, not both: they would each write the case's report, and simulate stops with an error if a model has both.

Fields

detection_prob::Any
delay::Any
from::Any

EpiBranch.MinimumSize Type

julia

struct MinimumSize <: ObservationModel

julia

MinimumSize(min_size)

Only chains of at least min_size cases are recorded, as when only clusters of two or more cases are investigated. min_size is at least 1. The exact likelihood conditions chain sizes on being at least min_size (TruncatedChainSize), and the simulation-based likelihood leaves out simulated chains below it, so both describe the same recorded data. simulate and chain_statistics still return every chain.

It acts on whole chains rather than on individual cases, so it leaves the simulated cases themselves unchanged. A model has one observation model, so MinimumSize cannot be combined with PerCaseObservation.

Fields

min_size::Int64

EpiBranch.observe Function

julia

observe(base, _::NoObservation) -> Any

julia

observe(base_distribution, obs::ObservationModel)

Turn the true distribution of chain sizes into the distribution a surveillance system would see under the observation model obs. The result is itself a distribution, so it can be used in a likelihood in place of the true one.

NoObservation: returned unchanged.
PerCaseObservation (deprecated): under-reporting, each case detected with probability detection_prob (ThinnedChainSize).
MinimumSize: only chains of at least min_size cases are recorded (TruncatedChainSize).

To add a new observation model, define this together with apply_observation!.

Examples

julia

# only chains of at least 2 cases are recorded
observe(chain_size_distribution(NegBin(0.8, 0.5)), MinimumSize(2))

EpiBranch.ChainSizeLaw Type

julia

abstract type ChainSizeLaw <: Distributions.Distribution{Distributions.Univariate, Distributions.Discrete}

julia

ChainSizeLaw

Parent type of every chain-size distribution: Borel, GammaBorel, PoissonGammaChainSize, ChainSizeMixture, IndexChainSize, TruncatedChainSize and ThinnedChainSize. One loglikelihood(::ChainSizes, ::ChainSizeLaw) method covers all of them, so observe(law, obs) can be passed straight to loglikelihood for any observation model, built in or user-defined, as long as the law it returns subtypes ChainSizeLaw and has a logpdf method.

Fields

EpiBranch.ThinnedChainSize Type

julia

struct ThinnedChainSize{D<:Distributions.Distribution{Distributions.Univariate, Distributions.Discrete}, T<:Real} <: ChainSizeLaw

julia

ThinnedChainSize(base, detection_prob)

Distribution of observed chain sizes under under-reporting: each case in a chain whose true size follows base is detected independently with probability detection_prob, and the observed size is the number detected. A chain with nobody detected leaves no trace in the data, so the distribution is conditioned on at least one case being detected. Usually built for you by observe with a PerCaseObservation.

logpdf(d, obs) adds logpdf(base, n) + logpdf(Binomial(n, p), obs) over n >= obs until the remaining tail is negligible, then subtracts log(1 - P(0 detected)), where P(0 detected) = Σ_n P(base = n) (1 - p)^n is summed the same way. Only logpdf on the base is needed, so any base distribution works without a method of its own. P(0 detected) depends on base and detection_prob alone, so it is worked out once when the distribution is built. Its field takes whatever element type base's parameters give logpdf, rather than Float64, so a Dual can pass through it when something differentiates through base, such as fitting an offspring parameter with ForwardDiff.

Fields

base::Distributions.Distribution{Distributions.Univariate, Distributions.Discrete}
detection_prob::Float64
log_prob_any_detected::Real

EpiBranch.TruncatedChainSize Type

julia

struct TruncatedChainSize{D<:Distributions.Distribution{Distributions.Univariate, Distributions.Discrete}} <: ChainSizeLaw

julia

TruncatedChainSize(base, min_size)

Chain size distribution for clusters recorded only once they reach min_size cases (for example when only clusters of two or more cases are investigated): the distribution base, conditioned on the size being at least min_size,

julia

P(N = n | N ≥ min_size) = P(N = n) / P(N ≥ min_size)

and zero below min_size. observe builds it from a MinimumSize observation model. Clusters with several index cases are conditioned using their own number of index cases.

Fields

base::Distributions.Distribution{Distributions.Univariate, Distributions.Discrete}
min_size::Int64

Reproduction number varying between chains

EpiBranch.ClusterMixed Type

julia

struct ClusterMixed{F, D<:Distributions.Distribution}

julia

ClusterMixed(build, mixing)

Offspring where the reproduction number varies from chain to chain rather than from case to case, for example between settings or clusters: each chain draws a value θ from mixing, and every case in that chain draws its number of secondary cases from build(θ). Use it in a BranchingProcess, or directly with loglikelihood and chain_size_distribution.

build is a function of θ returning an offspring distribution, or the Poisson family itself. Poisson offspring with a Gamma-distributed rate has an exact chain size distribution, which is used automatically; other combinations are evaluated numerically by integrating over mixing.

Examples

julia

# Poisson offspring whose rate is Gamma distributed between chains
# (shape 2, mean 0.8): exact chain size distribution
o = ClusterMixed(Poisson, Gamma(2.0, 0.4))
loglikelihood(ChainSizes([1, 2, 1, 5]), o)

# negative binomial offspring (k = 0.5) with R varying between chains:
# evaluated numerically
o = ClusterMixed(R -> NegBin(R, 0.5), Gamma(2.0, 0.3))
loglikelihood(ChainSizes([1, 1, 3, 2]), o)

Fields

build::Any
mixing::Distributions.Distribution

EpiBranch.ChainSizeMixture Type

julia

struct ChainSizeMixture{F, D<:Distributions.Distribution} <: ChainSizeLaw

julia

ChainSizeMixture(build, mixing)

Chain size distribution when the reproduction number varies between chains (ClusterMixed) and no exact formula is available: the chain size probabilities of build(θ), averaged over θ drawn from mixing. chain_size_distribution returns it.

The average is computed by numerical integration over the central 99.8% of mixing (its 0.001 to 0.999 quantiles).

Fields

build::Any
mixing::Distributions.Distribution

Pairwise survival likelihood

The likelihood of data on who was exposed to whom and when, used to estimate the contact interval or transmission rate from household or contact-network studies with individual infection times. The contact interval is the time from the start of the infector's infectious period to a contact that would infect if nothing intervened.

EpiBranch.pairwise_surv_loglik Function

julia

pairwise_surv_loglik(
    kernel,
    data::PairwiseSurvivalData
) -> Any

julia

pairwise_surv_loglik(kernel, data::PairwiseSurvivalData) -> Float64

Log-likelihood of exposure data in survival form, summed over who could have infected whom: use it to estimate the contact interval from rows of PairwiseSurvivalData. Each infected person contributes the log of the total hazard from their possible infectors at their infection time, and every row subtracts the cumulative hazard over its exposure interval:

julia

ll = Σ_susceptible log Σ_{event rows} hazard(stop)
     − Σ_rows [cumhazard(stop) − cumhazard(start)]

People never infected are right-censored: they contribute only the probability of escaping infection. kernel is the contact-interval distribution shared by every row, or a function r -> Distribution of the row number for covariates. The result is differentiable in the kernel's parameters, so it can be maximised with Optim.jl or added to a Turing model with @addlogprob!.

julia

pairwise_surv_loglik(
    kernel,
    data::InfectionLayer,
    layout::ContactPairsLayout;
    external_hazard,
    susceptibility
) -> Any

julia

pairwise_surv_loglik(kernel, data::InfectionLayer, layout::ContactPairsLayout;
                     external_hazard = 0.0, susceptibility = nothing) -> Real
pairwise_surv_loglik(kernel, data::InfectionLayer; external_hazard = 0.0,
                     susceptibility = nothing) -> Real

Log-likelihood of an outbreak in households or on a contact network, given who was infected when and who could have infected whom: use it to estimate the contact interval (and so the transmission rate) from such data. The data are an InfectionLayer, for example from household_infections or network_infections. The likelihood sums over who infected whom, so the transmission tree need not be known.

Each person accumulates hazard from every possible infector while that infector is infectious and they are still uninfected, and each infected person adds the log of the total hazard at their infection time. If an infected person (other than one the likelihood conditions on) faces zero total hazard at their infection time, because no possible infector can infect them then and there is no infection from outside, the infection times are impossible under the model and the result is -Inf.

kernel is the contact interval, in days from the start of the infector's infectious period. It can be one distribution shared by every pair, a function (infector, susceptible) -> Distribution of the two people's numbers for covariates, a PairKernel that also uses the infector's infection time (and, with a state, each person's record), or a vector of vectors parallel to a contact list (kernel[i][k] for person i's k-th listed contact).

external_hazard is infection from outside the households or network: a constant rate per person per day, or a distribution of the time of infection from outside. Either is shared by everyone, or differs between people: a vector with one rate or distribution per person, or a function j -> rate or distribution of person j's number, for example to give each arm of a trial its own community hazard. A function is taken to be in force, so the likelihood explains index cases as described next even if it returns zero for everyone. Such infections happen only in the first data.obs_end days. With an external hazard, index cases are explained like any other case; without one, the likelihood conditions on them. Each person is exposed to the external hazard until the earlier of their infection and data.obs_end, so a person infected after obs_end must have been infected by a possible infector. Spread within households or the network continues after obs_end: a person never infected is exposed over each possible infector's whole infectious period.

susceptibility is an effect on people's own risk of infection, such as a candidate VaccineEffect or a model's interventions; nothing (the default) leaves every hazard as the kernel gives it. Through susceptibility_components the effect says how it scales every hazard a person faces, from possible infectors and from outside alike. A VaccineEffect reads each person's immunity time from the record's host_times, which household_infections (in EpiHouseholds) and network_infections (in EpiNetwork) record for a model with vaccination; its efficacy must be a number. Under LeakyMode, a vaccinated person's hazards from their immunity time on are multiplied by 1 - efficacy, or by 1 - efficacy * waning(dt) with waning. Under AllOrNothingMode, their contribution is the mixture efficacy * Lᵖ + (1 - efficacy) * Lᵘ: Lᵖ the likelihood of their escapes and infection if fully protected from their immunity time on (zero if they were infected after it), and Lᵘ if unprotected. Whether the vaccine protects someone is decided once per person and holds for all their exposures, so the escape from all their possible infectors sits inside the mixture.

Everything is cut at followup_end(data): a person infected after it counts as uninfected until then, and no exposure counts after it. This gives the same value as first cutting the data there: later infections unobserved, and removal times and obs_end capped at it.

In inference, build the layout once with compile_contact_pairs and reuse it while the imputed infection times change; its external setting must agree with external_hazard. The two-argument form builds a layout on each call.

Example

Estimate the mean contact interval from a simulated household outbreak:

julia

using EpiBranch, EpiHouseholds, Distributions, StableRNGs
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.(grid))]   # close to the true mean of 4 days

Gradients

The likelihood is differentiable in the kernel's parameters. A Gamma kernel or community hazard needs a reverse-mode automatic differentiation backend such as Mooncake, because ForwardDiff cannot differentiate its cumulative hazard. Weibull and Exponential work with either.

A vanishing community hazard is not the no-community case

The two condition on different things, and the log-likelihood jumps between them at α = 0. With external_hazard = α > 0 an index case infected at time t contributes log(α) − αt, which falls to -Inf as α → 0, because a model that allows infection from outside has to explain the index cases it saw. At exactly external_hazard = 0 index cases are conditioned on and contribute nothing, leaving a finite value. A likelihood ratio between "some community transmission" and "none" therefore cannot be read off by letting α approach zero: evaluate the two models separately.

The jump is at that one point. Near it, the log-likelihood is k log α − αT up to terms free of α, where k counts the cases only the community can explain and T is the total time people are exposed to it. In log α this is a straight line of slope k.

julia

pairwise_surv_loglik(windows, data::RoutedNetworkInfections; external_hazard = 0.0) -> Real

Log-likelihood of a routed-network outbreak under windows, a collection of RouteWindows in the same order as data's routes, giving each route's kernel and which states (including EpiBranch.INTERVENTION_REMOVAL) end it.

Each susceptible's survival is the sum of every route's cumulative hazard from every possible infector on it (routes are independent competing risks), and an infected person's event term is the log of the summed hazard, over every route and every possible infector, at their infection time — the marginal density of whichever route and infector's draw came first in the race. A route whose until does not list EpiBranch.INTERVENTION_REMOVAL keeps exposing its susceptibles through a stretch that an intervention such as Isolation or Quarantine removed the infector for on another route; one that does list it has that stretch, read from data's host_times as the single-route likelihood reads it, taken out of its own exposure.

external_hazard is community infection shared across every route, as for NetworkProcess's likelihood.

Example

julia

using EpiNetwork, EpiBranch, Distributions
household = [[2, 3], [1, 3], [1, 2], Int[]]
community = [[4], Int[], Int[], [1]]
windows = [
    RouteWindow(:household; until = (:recovered,), kernel = Exponential(4.0), reach = household),
    RouteWindow(:community; until = (:recovered,), kernel = Exponential(20.0), reach = community),
]
progression = [Transition(:recovered; from = :infection, delay = 6.0, terminal = true)]
state = simulate(ModelSpec(RoutedNetwork(windows); progression); n_initial = 2)
data = routed_network_infections(state, RoutedNetwork(windows))
pairwise_surv_loglik(windows, data)

EpiBranch.PairwiseSurvivalData Type

julia

struct PairwiseSurvivalData{T<:Real}

julia

PairwiseSurvivalData(sus, start, stop, event)

Exposure data in survival-analysis form for pairwise_surv_loglik: one row per susceptible person and possible infector, giving the time over which the person was exposed to that infector and whether they were infected at the end of it. Row r covers the exposure interval (start[r], stop[r]] (days) for susceptible person sus[r], and event[r] is true if they were infected at stop[r]. The rows record only exposure and infection, without who lives where or the order of infections.

Fields

sus::Vector{Int64}
start::Vector{T} where T<:Real
stop::Vector{T} where T<:Real
event::Vector{Bool}

EpiBranch.InfectionLayer Type

julia

abstract type InfectionLayer

julia

InfectionLayer

The record of who was infected when in an outbreak, together with who could have infected whom (household members, or neighbours in a contact network). This is the data pairwise_surv_loglik scores to estimate the contact interval. A subtype holds, for each person i (numbered 1:n):

infection_time[i]: the infection time in days, NaN if never infected;
infectious_time[i]: when they became infectious;
removal_time[i]: when they stopped being infectious (recovery, death or isolation);
is_index[i]: whether they were infected from outside the households or network;

and obs_end, the day after which no more infections from outside arrive (only used when there is a community hazard); spread within the households or network continues after it. A subtype may also hold host_times, a named tuple of further per-person time vectors such as onset_time, which a PairKernel state function or a vaccine effect reads in the likelihood as it reads individual.state in simulation. missing marks a person without that time; NaN is a recorded value, as for the onset of an asymptomatic case in a simulation.

Follow-up ends at followup_end: a followup_end field when the subtype has one, and Inf otherwise. The likelihood ignores everything after it. A person infected later counts as uninfected until then, exposure stops there, and a case still infectious at the end of follow-up can keep a removal time of Inf. Simulated outbreaks run to the end and need no end of follow-up.

In real data the infection times are usually unobserved; they are known exactly after a simulation and imputed (data augmentation) in inference. What is observed, such as onsets and test results, comes from the natural history and is scored separately by progression_loglik; the sum of the two is the full log-likelihood of the augmented data. There is no closed-form likelihood of the onsets alone, because the unobserved infection times cannot be integrated out exactly.

household_infections (in EpiHouseholds) and network_infections (in EpiNetwork) read a simulated outbreak into an InfectionLayer. Each infected person's infectious period starts at the process's from state and ends at the earliest of its until states and the time the model's interventions take them out of transmission, such as isolation or quarantine after tracing; these are the same periods the simulation used. Any other change to transmission must be built into the kernel passed to the likelihood, except changes to a person's own susceptibility, which the model's interventions declare through susceptibility_components. loglikelihood(data, model) for households and networks checks the model's interventions with infection_likelihood_compatible; where it refuses, call pairwise_surv_loglik with a kernel that includes the extra effects. Passing followup_end to household_infections or network_infections evaluates the outbreak as if observation had stopped then.

To define a new kind of contact structure, a subtype also defines contact_structure; compile_contact_pairs and pairwise_surv_loglik then work on it with no further methods. HouseholdInfections and NetworkInfections are the worked examples.

Fields

EpiBranch.compile_contact_pairs Function

julia

compile_contact_pairs(
    membership::AbstractVector{<:Integer},
    is_index::AbstractVector{Bool},
    infected::AbstractVector{Bool};
    external
) -> ContactPairsLayout

julia

compile_contact_pairs(membership::AbstractVector{<:Integer}, is_index, infected; external = false)
compile_contact_pairs(contacts::AbstractVector{<:AbstractVector{<:Integer}}, is_index, infected; external = false)
compile_contact_pairs(data::InfectionLayer; external = false)

List once who could have infected whom, as a ContactPairsLayout to reuse across likelihood evaluations during inference.

The contact structure is either a membership vector, where people sharing a label can all infect one another (such as households), or a list of contacts, where contacts[i] lists the people i can infect (a directed network; list each contact both ways for an undirected one). A person's possible infectors are their household members or the people listing them. A contact listed twice counts as two separate contacts.

infected is true for each person infected in every set of infection times the layout will be used with (their infection time may still be imputed). Only infected people can infect others. is_index marks people infected from outside: without a community hazard the likelihood conditions on them and they appear only as infectors. With external = true every infection is explained, and each person gets an extra row for infection from outside. The single-argument form reads the contact structure from data with contact_structure and takes as infected everyone with a non-NaN infection_time.

EpiBranch.ContactPairsLayout Type

julia

struct ContactPairsLayout

julia

ContactPairsLayout

The list of who could have infected whom that the pairwise likelihood is evaluated over, built once with compile_contact_pairs and reused while infection times change during inference. Each row is one ordered (susceptible, possible infector) pair, plus, when there is a community hazard, one row per person for infection from outside. Rows whose exposure periods do not overlap are kept and skipped when evaluating, so one layout works for every set of infection times with the same contact structure and the same people infected.

component gives each person's group in the contact structure (their household, or their connected part of a contact network), numbered 1:ncomponents. pairwise_surv_loglik_by_component uses it to split the likelihood by group, for a sampler that updates infection times one group at a time.

Fields

sus::Vector{Int64}
infector::Vector{Int64}
contact_index::Vector{Int64}
is_ext::Vector{Bool}
sus_unique::Vector{Int64}
sus_row_ranges::Vector{UnitRange{Int64}}
sus_row_order::Vector{Int64}
no_rows::Vector{Int64}
external::Bool
nhosts::Int64
component::Vector{Int64}
ncomponents::Int64

EpiBranch.PairKernel Type

julia

struct PairKernel{F, S, C, W<:Tuple}

julia

PairKernel(callback; state = nothing, calendar = nothing, watches = nothing)

A contact interval that depends on who the two people are or when the contact happens: for example, faster transmission from older cases, or a lockdown that halves transmission from day 30. callback returns the contact-interval distribution (in days, measured from the start of the infector's infectious period), and calendar optionally scales the rate of transmission by calendar day. It works with NetworkProcess, HouseholdProcess and the InfectionLayer forms of pairwise_surv_loglik, so the same kernel can be simulated and fitted.

With state = nothing, callback(context::PairContext) receives a PairContext (the two people's numbers and the infector's infection time) and returns a distribution, or a (profile = ..., calendar = ...) named tuple. Use this when the kernel reads only fixed characteristics and the infector's infection time:

julia

kernel = PairKernel(context -> Exponential(exp(0.1 * context.infector_infection_time)))

With state given, callback(context, source, target) also receives a record for each of the two people. In simulation, state is a function of the individual that builds this record. For the likelihood, pass a vector of records indexed by person number as state, or extract them from a simulation with record_kernel; the callback is the same in both.

watches lists every key of individual.state that the state function reads, as a tuple of Symbols (see EpiBranch.watched_records). When one of these values changes for a person during a simulation, for example when they are vaccinated, the contacts not yet made that depend on them are redrawn under the new rate. A change to a key the function reads but watches leaves out stops the simulation with an error. List every key the function reads, even if nothing in the current model sets it. watches is required when state is a function; pass () when the function reads no state key (for example one that looks people up by ind.id). Anything that can change during the outbreak must be read from individual.state, which is the only place changes are followed, and not from a field such as ind.susceptibility. With no state, or with a vector of records, nothing can change, and passing a non-empty watches is an error.

calendar, a Steps or Seasonal schedule, a callable t -> multiplier, or any other type with a calendar_multiplier method, multiplies the transmission rate by the schedule's value on the calendar day (the start of the infector's infectious period plus the time since). A pair whose schedule depends on a person's record returns it from the callback instead, as (profile = ..., calendar = ...), which overrides the shared calendar for that pair.

julia

PairKernel((ctx, source, target) -> Gamma(2.0, 1.5);
           calendar = Steps([30.0], [1.0, 0.25]),   # shared policy: rate falls to a quarter on day 30
           state = ind -> (age = ind.state[:age],),  # optional per-person record
           watches = (:age,))

PairKernel((ctx, source, target) ->
               isfinite(target.date) ?
               (profile = Exponential(2.5), calendar = Steps([target.date], [1.0, 0.25])) :
               Exponential(2.5);
           state = ind -> (date = get(ind.state, :policy_time, Inf)::Float64,),
           watches = (:policy_time,))

The rate at which one person infects another is the contact-interval distribution's hazard at the time since the infector became infectious, times the schedule's value on that calendar day. For a step schedule such as Steps both simulation and the likelihood compute this exactly.

Records must keep their history

The callback must describe a hazard that only changes forwards in time: recording an event at day t must not change the hazard before t. A final "vaccinated" flag is not enough; keep the vaccination date, so the rate before it is unchanged. Neither the callback nor the state function may change anything, and the kernel may read a person's information only through state, since those records are all the likelihood is given.

A simulation in which no watched record changes draws outbreaks from the same distribution as an ordinary kernel. In a household model, a kernel with watches simulates all households together on one timeline, so that a policy can depend on cases in other households; this gives the same distribution of outbreaks from a different sequence of random numbers.

A plain distribution remains the simplest kernel and needs no PairKernel.

Fields

callback::Any
state::Any
calendar::Any
watches::Tuple

EpiBranch.watched_records Function

julia

watched_records(kernel) -> Tuple{}

julia

watched_records(kernel)

The keys of individual.state that kernel's transmission rate depends on, as a tuple of Symbols. The default, (), is a kernel whose rates are fixed for the whole outbreak.

When one of these values changes for a person during a simulation (a vaccination date being set, say), the contacts already drawn from a kernel that lists that key are drawn again under the new rate, given the exposure that has already happened. A kernel that lists nothing is never redrawn.

Declare every key the kernel reads, whether or not anything in a given model writes it, and read anything that can move from individual.state rather than from a field of the individual, which this cannot name. A PairKernel reports its watches, and () once its records are extracted; any other kernel type declares its own method.

EpiBranch.Steps Type

julia

struct Steps{T<:Real}

julia

Steps(breaks, values)

A change in transmission on given calendar days, such as a lockdown or school closure that scales everyone's rate of infecting others. breaks are the days (since the start of the outbreak) on which the multiplier changes and values the multipliers: values[1] before breaks[1], values[k+1] from breaks[k] (inclusive) up to breaks[k+1], and values[end] from breaks[end] onwards. breaks must be finite and strictly increasing, and every value non-negative; a value of zero stops transmission until the next break, or for good if it is the last value.

Pass it as the calendar of a PairKernel, which multiplies the contact-interval hazard by the value in force on each calendar day.

Example

Transmission halves from day 30:

julia

using EpiBranch, Distributions
lockdown = Steps([30.0], [1.0, 0.5])
kernel = PairKernel(ctx -> Exponential(2.0); calendar = lockdown)

Fields

breaks::Vector{T} where T<:Real
values::Vector{T} where T<:Real

EpiBranch.Seasonal Type

julia

struct Seasonal{T<:Real}

julia

Seasonal(; amplitude, peak_day, period = 365.2425)

A smoothly oscillating multiplier on calendar time: 1 + amplitude * cos(2π * (t - peak_day) / period). amplitude must be between 0 and 1, which keeps the multiplier non-negative; peak_day is the calendar time of the highest multiplier, and period defaults to a year.

Used as the calendar a PairKernel multiplies its contact-interval hazard by, on the calendar-time axis. It declares EpiBranch.SmoothCalendar, so its cumulative hazard is integrated by quadrature rather than exactly as Steps's is. amplitude and peak_day share a type parameter so the likelihood can be differentiated with respect to either.

Fields

amplitude::Real
peak_day::Real
period::Float64

EpiBranch.record_kernel Function

julia

record_kernel(
    k,
    state::SimulationState
) -> PairKernel{F, Nothing} where F

julia

record_kernel(kernel, state::SimulationState)

Prepare a PairKernel for fitting to a simulated outbreak: returns a kernel with the same callback and calendar whose state is the vector of each person's record at the end of state. Each record is copied, so later changes to the simulation do not alter it. Other kernels are returned unchanged.

Only the record the state function returns at the end of the run is kept; values overwritten during the run are lost. Where rates change, the function must return event dates or full histories. The likelihood then gives the transmission part along these histories; the probabilities of the characteristics and intervention events themselves need their own terms if they belong in the joint likelihood. During inference, build PairKernel(callback; state = records) from the current values of unobserved records at each evaluation.

EpiBranch.PairContext Type

julia

struct PairContext{T<:Real}

julia

PairContext(infector, susceptible, infector_infection_time)

What a PairKernel function knows about one infector and the person they may infect, in simulation and in the likelihood alike: the two people's numbers in the population (infector, susceptible) and the infector's infection time (days since the start of the outbreak).

The susceptible person's own infection time is left out because it is not yet known when their contact interval is drawn. To use fixed characteristics such as age or household, look them up by person number in a table the function refers to.

Fields

infector::Int64
susceptible::Int64
infector_infection_time::Real

EpiBranch.LayerHost Type

julia

struct LayerHost{T, S}

julia

LayerHost

One person of an InfectionLayer, as a PairKernel state function sees them when the likelihood is evaluated. It mirrors a simulated individual: its population id, its infection_time (NaN if never infected), and a state holding the person's recorded event times under their keys. state reads like an individual's: state[key] and get(state, key, default) give the recorded time, including a recorded NaN, and a person whose entry is missing has none, so get returns the default and state[key] throws. Reading a key the infection record does not hold throws an ArgumentError, so a state function cannot silently fall back to a default for a time the likelihood was never given.

Fields

id::Int64
infection_time::Any
state::Any

Household data

EpiHouseholds.household_infections Function

julia

household_infections(state, model::ModelSpec; obs_end = model.process.obs_end,
                     followup_end = Inf, host_times = ()) -> HouseholdInfections

Collect who was infected when from a household outbreak simulated from model (infection times, start and end of each infectious period, index cases), in the form the pairwise likelihood needs. The infectious periods are the ones the simulation used, as described for InfectionLayer. Besides the infectious periods, the event times listed below are recorded, and loglikelihood applies interventions that change susceptibility, such as vaccination, from them. Any other effect on transmission must be built into the kernel passed to the likelihood. A bare HouseholdProcess is accepted too (its infectious period starts at :infection, and it has no interventions).

host_times names further per-person event times to record, such as (:onset_time,), read from each person's state (missing where a person has none), for a PairKernel to use. The times the model's interventions need (see susceptibility_host_times), such as a vaccinee's :immunity_time, are recorded as well.

EpiHouseholds.HouseholdInfections Type

julia

HouseholdInfections(household_of, infection_time, infectious_time, removal_time, is_index;
                    obs_end = Inf, followup_end = Inf, host_times = (;))

Household outbreak data for the pairwise likelihood: for each person, their household, when they were infected, when their infectious period started and ended, and whether they were an index case. Household members are each other's possible infectors. household_of gives each person's household; the other per-person vectors, obs_end, followup_end and host_times are as described for InfectionLayer. Read one from a simulation with household_infections, or build it from data, imputing unobserved infection times in inference.

EpiHouseholds.HouseholdPairsLayout Type

julia

HouseholdPairsLayout

The list of who could have infected whom in a household population, for the pairwise likelihood: each row is one ordered (susceptible, household member) pair. It is another name for EpiBranch's ContactPairsLayout, built from households. Build it with compile_household_pairs.

EpiHouseholds.compile_household_pairs Function

julia

compile_household_pairs(household_of, is_index, infected; external=false)
compile_household_pairs(data::HouseholdInfections; external=false, condition_on=RecruitedIndex())

List once who could have infected whom in a household population, where household members are each other's possible infectors; this is compile_contact_pairs for households, with the arguments described there. Evaluate the result with pairwise_surv_loglik(kernel, data, layout; external_hazard).

condition_on is a ConditionOn rule saying how to condition on the index case: which member of each household the likelihood takes as given, when there is no community hazard (external = false; with one, every infection is explained and the rule has no effect). RecruitedIndex, the default, conditions on the recruited index case, data.is_index, fixed when the data are read. EarliestInfected conditions on the household member with the earliest infection_time in data, worked out on every call. Use it when the recruited index case need not be the first member infected, which is expected in real recruited households; otherwise an imputed infection time earlier than the recruited index case's makes the likelihood impossible (-Inf). Because the conditioned member can change between calls, build a fresh layout for EarliestInfected on every evaluation rather than reusing one across imputed infection times.

EpiHouseholds.ConditionOn Type

julia

ConditionOn

How the household likelihood conditions on the index case: which member of each household is taken as given rather than explained by transmission within the household. RecruitedIndex and EarliestInfected are the two supplied; a rule of your own needs only a condition_mask method.

EpiHouseholds.RecruitedIndex Type

julia

RecruitedIndex()

Condition each household on its recruited index case, data.is_index, as recorded.

EpiHouseholds.EarliestInfected Type

julia

EarliestInfected()

Condition each household on whichever member has the earliest infection_time (on a tie, the lowest person number). This is worked out from data on every call, so the conditioned member can change as imputed infection times change.

Network data

EpiNetwork.network_infections Function

julia

network_infections(state, model::ModelSpec{<:NetworkProcess};
                   obs_end = model.process.obs_end, followup_end = Inf,
                   host_times = ()) -> NetworkInfections
network_infections(state, process::NetworkProcess) -> NetworkInfections

Collect who was infected when from a network outbreak simulated from model (infection times, start and end of each infectious period, index cases), in the form the pairwise likelihood needs, with the model's network as the contact structure. The infectious periods are the ones the simulation used, as described for InfectionLayer. Besides the infectious periods, the event times listed below are recorded, and loglikelihood applies interventions that change susceptibility, such as vaccination, from them. Any other effect on transmission must be built into the kernel passed to the likelihood. A bare NetworkProcess is accepted too (its infectious period starts at :infection, and it has no interventions).

EpiNetwork.NetworkInfections Type

julia

NetworkInfections(contacts, infection_time, infectious_time, removal_time, is_index;
                  obs_end = Inf, followup_end = Inf, host_times = (;))

Network outbreak data for the pairwise likelihood: for each person, when they were infected, when their infectious period started and ended, and whether they were infected from outside the network. contacts[i] lists the people i can infect, as for NetworkProcess, and a person's possible infectors are those who list them. The per-person vectors, obs_end, followup_end and host_times are as described for InfectionLayer. Read one from a simulation with network_infections, or build it from data, imputing unobserved infection times in inference.

EpiNetwork.routed_network_infections Function

julia

routed_network_infections(state, model::ModelSpec{<:RoutedNetwork};
                          obs_end = model.process.obs_end, followup_end = Inf)
    -> RoutedNetworkInfections
routed_network_infections(state, process::RoutedNetwork) -> RoutedNetworkInfections

Collect who was infected when from a routed-network outbreak simulated from model, with every route's own infectious window read back exactly as the race closed it, including a lapsing removal (such as an Isolation with a finite duration) stored in the record's host_times the same way network_infections stores it, for pairwise_surv_loglik to read back out on whichever routes list EpiBranch.INTERVENTION_REMOVAL.

EpiNetwork.RoutedNetworkInfections Type

julia

RoutedNetworkInfections(route_names, contacts, infectious_time, removal_time,
                        infection_time, is_index; obs_end = Inf,
                        followup_end = Inf, host_times = (;))

Routed-network outbreak data for the pairwise likelihood: an infection record per person, with a separate infectious window and contact list for each route, in the order route_names gives. contacts[r][i] lists the people i can infect on route r, as RouteWindow's reach; infectious_time[r] and removal_time[r] are when route r opened and closed for each person, the same windows RoutedNetwork's race used. The remaining arguments are as for NetworkInfections: infection_time, is_index, obs_end, followup_end and host_times are shared across routes, since they describe the person rather than any one route. Read one from a simulation with routed_network_infections.

Distribution helpers

EpiBranch.NegBin Function

julia

NegBin(R::Real, k::Real) -> Distributions.NegativeBinomial

julia

NegBin(R, k)

A negative binomial offspring distribution with mean R (the reproduction number) and dispersion k, so the variance is R + R²/k. Smaller k means more superspreading: k below 1 is strong superspreading, and as k grows the distribution approaches a Poisson. Returns a NegativeBinomial from Distributions.jl.

Warning

NegativeBinomial(r, p) from Distributions.jl takes a number of successes and a success probability, not a mean and dispersion. NegativeBinomial(R, k) therefore gives a different distribution, without any error when k is at most 1. Use NegBin(R, k).

Examples

julia

NegBin(2.5, 0.16)   # R = 2.5, k = 0.16 (k as estimated for SARS)

EpiBranch.scale_distribution Function

julia

scale_distribution(
    d::Distributions.Poisson,
    factor::Real
) -> Distributions.Poisson

julia

scale_distribution(d, factor::Real)

Multiply the mean of an offspring distribution, the reproduction number, by factor, keeping the same family and dispersion: for example to model a control measure that cuts transmission by a fixed proportion. Defined for Poisson (Poisson(R * factor)) and the negative binomial (same dispersion k, mean R * factor); any other distribution raises an error.

Examples

julia

scale_distribution(NegBin(2.5, 0.16), 0.5)   # R = 1.25, k = 0.16

EpiBranch.incubation_linked_generation_time Function

julia

incubation_linked_generation_time(
;
    presymptomatic_fraction,
    omega
) -> EpiBranch.var"#15#16"{Float64, Float64}

julia

incubation_linked_generation_time(; presymptomatic_fraction=0.3, omega=2.0)

A generation time linked to each case's own incubation period, as in Hellewell et al. (2020), for use as the generation_time of a BranchingProcess. Incubation periods come from clinical_presentation.

presymptomatic_fraction: the share of transmission that happens before the infector's symptom onset.
omega: the spread, in days, of the generation time around the incubation period.

Each case's generation time is drawn from a skew-normal distribution centred on its incubation period, with scale omega and the skew chosen so that presymptomatic_fraction of generation times are shorter than the incubation period. Negative values are excluded, which makes the realised presymptomatic share approximate: close when the incubation period is long compared with omega, less so for short ones. Cases without an incubation period (for example asymptomatic cases) are centred on 5 days.

Examples

julia

model = BranchingProcess(
    NegBin(2.5, 0.16),
    incubation_linked_generation_time(presymptomatic_fraction=0.3)
)

Default placeholders

The values shown as defaults in signatures when a setting is left out, such as max_cases = NoCases().

EpiBranch.NoPopulation Type

julia

struct NoPopulation

Placeholder meaning no limit on the population size (an unbounded population).

Fields

EpiBranch.NoAttributes Type

julia

struct NoAttributes

Placeholder meaning no population characteristics are drawn.

Fields

EpiBranch.NoTypeLabels Type

julia

struct NoTypeLabels

Placeholder meaning the types of a multi-type model have no names.

Fields

EpiBranch.NoAgeDistribution Type

julia

struct NoAgeDistribution

Placeholder meaning no age distribution is given (ages are uniform over the age range).

Fields

EpiBranch.NoCases Type

julia

struct NoCases

Placeholder meaning no case cap when deciding whether an outbreak was contained.

Fields

EpiBranch.NoGenerationTime Type

julia

struct NoGenerationTime

Placeholder meaning no generation time: a BranchingProcess without timing, used for chain sizes and lengths only.

Fields

For extension authors

The entries below are for writing a new transmission model, intervention, transition or observation model, and are not needed to run the models above. Extending EpiBranch explains how they fit together.

Writing a transmission model

A transmission model defines initialise_state to set up its starting population. The helpers listed after it create the SimulationState and add the index cases.

A model that runs its own simulation loop has two more to build on. simulate_once runs one outbreak with the modelling layers passed in rather than read off the model, and sellke_race! settles one continuous-time race and brings the state's running totals back into line afterwards.

EpiBranch.TransmissionModel Type

julia

abstract type TransmissionModel

The parent type of all transmission models: BranchingProcess, HomogeneousProcess, MixingProcess, NetworkProcess, RoutedNetwork and HouseholdProcess. Needed only to write a new kind of model; the Extending EpiBranch guide lists what a new model defines.

Fields

EpiBranch.draw_offspring Function

julia

draw_offspring(
    rng::Random.AbstractRNG,
    offspring::Distributions.Distribution,
    individual,
    state::SimulationState
) -> Any

Draw the number of secondary cases of one case from an offspring distribution.

julia

draw_offspring(rng, offspring, individual, state) -> Any

Draw the number of secondary cases of one case from an offspring function, written either as (rng, individual) -> n or as (rng, individual, state) -> n. The second form can read the state of the whole outbreak, for example to reduce transmission once the case count passes a threshold.

julia

draw_offspring(
    rng::Random.AbstractRNG,
    offspring::EpiBranch.MultiTypeOffspring,
    individual,
    state::SimulationState
) -> Any

julia

draw_offspring(rng, offspring::MultiTypeOffspring, individual, state)

Draw the number of secondary cases of each type for a case of type j: a total from dist_fn(R_j), split at random between types in proportion to column j of the next-generation matrix.

julia

draw_offspring(
    rng::Random.AbstractRNG,
    offspring::ClusterMixed,
    individual,
    state::SimulationState
) -> Any

julia

draw_offspring(rng, offspring::ClusterMixed, individual, state)

Draw the number of secondary cases of one case when the reproduction number varies between chains: θ is drawn once per chain, at the index case, and every case in the chain draws from build(θ).

EpiBranch.reads_settled_state Function

julia

reads_settled_state(_) -> Bool

julia

reads_settled_state(x) -> Bool

Whether an offspring specification or a generation time x can read anything settling a case writes: a time a Transition records, a value an intervention sets, or the state of the whole outbreak, such as its case count. A branching process asks this before a ring of contact tracing reaches a person whose infection is not yet decided. If neither the window's offspring nor its generation time reads settled state, the person's contacts are drawn then and the ring grows through them, whether or not the person is infected. Otherwise they are drawn only once the person settles as a case, so an uninfected ring member has none.

The default is true, the safe answer for a function, which may read anything. A Distribution, a MultiTypeOffspring (which reads only the case's type, set when it is created), a ClusterMixed (which reads its parent's draw) and no generation time answer false. An offspring specification of your own that draws from nothing settling writes declares EpiBranch.reads_settled_state(::YourOffspring) = false.

EpiBranch.simulate_once Function

julia

simulate_once(model::HouseholdProcess, sim_opts; interventions, attributes,
              progression, observation, recorder, rng)

Simulate a household outbreak exactly in continuous time; this is the model the pairwise likelihood describes. The natural history and interventions come from the caller (a bare process, or a ModelSpec), and the start of the infectious period (from) is taken from the progression. Returns an EpiBranch SimulationState; linelist(state) turns it into a table with one row per case.

With no external hazard each household starts with one index case at day 0 and spreads only within the household. With an external_hazard (a positive rate or a distribution of the time of infection from outside), community introductions happen over the first model.obs_end days, which must be finite. Chosen initial_cases are infected at day 0 regardless, and any external hazard still acts on everyone else from day 0.

julia

simulate_once(
    model::TransmissionModel,
    sim_opts::SimOpts;
    interventions,
    attributes,
    progression,
    observation,
    recorder,
    rng
)

julia

simulate_once(model, sim_opts; interventions, attributes, progression,
              observation, recorder, rng)

Run a single outbreak simulation against a built SimOpts, with the modelling layers (interventions, attributes, progression, observation, recorder) supplied explicitly rather than read off model. simulate reads them off a bare process, or off a ModelSpec, and wraps this in a retry loop when a condition is given; this method itself takes no condition.

The built-in models run on the continuous-time race and each define a method, reading the layers off the arguments and deriving any window state they need (the infectious window's from, say, via infectious_from) from the composed progression. The default method runs the generation-based engine, for a model written against generate_offspring or contacts_of. recorder is read only by the race; the generation-based engine never drops a pair for a standing block, so it has nothing to ask one.

julia

simulate_once(
    model::BranchingProcess,
    sim_opts::SimOpts;
    interventions,
    attributes,
    progression,
    observation,
    recorder,
    rng
)

julia

simulate_once(model::BranchingProcess, sim_opts; interventions, attributes,
              progression, observation, recorder, rng)

Simulate a branching process on the continuous-time race. Each case creates its contacts when it settles, one route per Infectiousness window, and the race decides every contact in order of infection time.

julia

simulate_once(
    model::HomogeneousProcess,
    sim_opts::SimOpts;
    interventions,
    attributes,
    progression,
    observation,
    recorder,
    rng
)

julia

simulate_once(model::HomogeneousProcess, sim_opts; interventions, attributes,
              progression, observation, recorder, rng)

Simulate the homogeneous population on the continuous-time race, as mass action among its members (_MassAction), with the modelling layers supplied by the caller (a bare process, or a ModelSpec). The infectious window's from state is resolved here from the composed progression.

julia

simulate_once(
    model::MixingProcess,
    sim_opts::SimOpts;
    interventions,
    attributes,
    progression,
    observation,
    recorder,
    rng
)

julia

simulate_once(model::MixingProcess, sim_opts; interventions, attributes,
              progression, observation, recorder, rng)

Simulate the structured population as mass action between its mixing groups, with the modelling layers supplied by the caller (a bare process, or a ModelSpec). The infectious window's from state is resolved here from the composed progression; the mixing groups are read off each individual's state, so the attributes composing the model must set model.mixing_by's keys before the race runs.

julia

simulate_once(model::NetworkProcess, sim_opts; interventions, attributes,
              progression, observation, recorder, rng)

Simulate a network outbreak exactly in continuous time, with the natural history and interventions supplied by the caller (a bare process, or a ModelSpec). The start of the infectious period (from) is taken from the progression. Returns an EpiBranch SimulationState; linelist(state) turns it into a table.

With no external hazard, sim_opts.n_initial different people are infected at day 0 and the outbreak spreads along the network. With an external_hazard (a positive rate or a distribution of the time of infection from outside), community introductions happen over the first model.obs_end days, which must be finite (otherwise everyone would eventually be infected from outside). Chosen initial_cases are infected at day 0 regardless, and any external hazard still acts on everyone else from day 0.

EpiBranch.generate_offspring Function

julia

generate_offspring(model, parent, state) -> count

How many contacts the case parent makes in this generation, before susceptibility and interventions decide which of them are infected: one number, or one number per type in a multi-type model. Defined by a new transmission model in which every contact is a new person. For such a model the simulation calls it once for each infectious case, then creates that many contacts and times each from the model's generation time, so the model itself does not create people, time infections or apply interventions. The built-in branching process draws its counts from its offspring distribution directly.

A model whose contacts are people who already exist, such as a contact network or households, defines contacts_of and a collect_exposures method that returns gather_by_target instead. No model in the package is written that way: the network and household models run the continuous-time race instead.

EpiBranch.contacts_of Function

julia

contacts_of(
    model::TransmissionModel,
    parent,
    state::SimulationState
)

julia

contacts_of(model, node, state) -> iterable of (contact, infection_time)

The people an infectious case node can reach in this generation, each with the time they would be infected, as (contact, infection_time) pairs.

Only needed when writing a new transmission model that spreads over a fixed set of people, such as a contact network (the case's network neighbours) or households (the other household members). Such a model defines this function for its own type, and a collect_exposures method that returns gather_by_target, so that a person exposed by several cases in one generation is infected at most once.

Models in which every contact is a new person, such as a branching process, define generate_offspring instead. In both cases the simulation applies interventions, natural history and the infection decision itself. Calling contacts_of on a model that does not define it raises an error explaining what to define.

EpiBranch.collect_exposures Function

julia

collect_exposures(
    model::TransmissionModel,
    state::SimulationState
) -> Tuple{Array{Individual{T}, 1} where T<:Real, Array{EpiBranch._OneEdge{T}, 1} where T<:Real, Union{SubArray{Any, 1, P, Tuple{UnitRange{Int64}}, true} where P<:(Array{Individual{T}, 1} where T<:Real), SubArray{E, 1, P, Tuple{UnitRange{Int64}}, true} where {E<:Individual, P<:(Array{Individual{T}, 1} where T<:Real)}}, EpiBranch._AllNew}

julia

collect_exposures(model, state) -> (targets, edges, minted, is_new)

The people each infectious case exposes in this generation, and when. Only needed when writing a new transmission model.

Returns the distinct exposed people (targets); for each of them, the (infector_id, infection_time) exposures reaching them (edges) and whether they were created in this generation (is_new); and the contacts newly created in this generation (minted).

The default suits models in which every contact is a new person: it asks each infectious case how many contacts it makes via generate_offspring, then creates and times that many contacts. Models whose contacts are people who already exist and can be exposed by more than one case in a generation (networks, households) define a method of this function that returns gather_by_target, which calls contacts_of and groups the exposures by person.

EpiBranch.gather_by_target Function

julia

gather_by_target(
    model::TransmissionModel,
    state::SimulationState
) -> Tuple{Array{Individual{T}, 1} where T<:Real, Array{Vector{var"#s183"}, 1} where var"#s183"<:Tuple{Int64, Real}, Union{SubArray{Any, 1, P, Tuple{UnitRange{Int64}}, true} where P<:(Array{Individual{T}, 1} where T<:Real), SubArray{E, 1, P, Tuple{UnitRange{Int64}}, true} where {E<:Individual, P<:(Array{Individual{T}, 1} where T<:Real)}}, Vector{Bool}}

julia

gather_by_target(model, state) -> (targets, edges, minted, is_new)

A version of collect_exposures for models in which several cases can expose the same person in one generation (networks, households). It collects every exposure of a person and lets the simulation decide once whether, and by whom, they are infected. Newly created contacts are each treated separately, so a model may mix new people with existing ones.

EpiBranch.model_generation_time Function

julia

model_generation_time(model::TransmissionModel) -> Any

julia

model_generation_time(model)

The generation time a model uses to time each case's secondary cases, read by the default collect_exposures. Returns the model's generation_time field. A new model that times its contacts another way (for example with a separate delay per transmission route) defines this function for its own type instead of having a generation_time field.

EpiBranch.transmission_risks Function

julia

transmission_risks(_::TransmissionModel) -> Tuple{}

julia

transmission_risks(model) -> iterable of risk sources

Ways in which the transmission model itself can stop a contact becoming infected, applied alongside susceptibility, infectiousness and the interventions. Only needed when writing a new model. For example, a metapopulation model in which transmission between two places succeeds only with some probability returns a source here, so that probability acts on every contact while the contact is still made (and can still be traced). Each source defines competing_risk. Returns none by default.

EpiBranch.race_groups Function

julia

race_groups(
    model::TransmissionModel,
    kernel
) -> Union{Tuple{Vector{Int64}}, Vector{Vector{Int64}}}

julia

race_groups(model, kernel)

For model authors: the groups of people a continuous-time model simulates separately, as disjoint groups of population ids, each simulated in its own call to _sellke_race! with its own random number stream. HouseholdProcess defines this to simulate each household on its own. A new kernel type can define a method for a given model to choose a different grouping.

EpiBranch.sellke_race! Function

julia

sellke_race!(
    state::SimulationState,
    members::AbstractVector{Int64},
    rng::Random.AbstractRNG,
    sim_opts::SimOpts;
    kwargs...
)

julia

sellke_race!(state, members, rng, sim_opts; kwargs...) -> state

Public wrapper around _sellke_race!: resolves the race's max_time, max_generation and stopping_rules from sim_opts's stopping rules, then reconciles state's aggregate bookkeeping (cumulative_cases, max_infection_time, extinct) from the per-individual state the race wrote directly, before returning state. A model with its own simulation loop calls this rather than _sellke_race! itself, so it never has to resolve these or reconcile the bookkeeping by hand. kwargs are _sellke_race!'s own (seed!, targets or routes, from, until, interventions, contacts, risks, introduction, watches, recorder), everything but the three this wrapper sets.

A model with more than one natural race partition (a household process, over its households, via race_groups) calls _sellke_race! directly in its own loop instead, with max_time, max_generation and stopping_rules resolved once outside it, and reconciles state once after the loop: this wrapper's per-call reconciliation scans every individual, which would cost O(races × population) run inside such a loop rather than the O(population) a single-race model gets from calling this once. See HouseholdProcess in the EpiHouseholds package.

EpiBranch.infectious_from Function

julia

infectious_from(progression) -> Symbol

julia

infectious_from(progression) -> Symbol

The state a case's infectious window opens at, derived from progression: :infectious when the progression produces it (a latent period), otherwise :infection. A model whose from field can override the derivation reads that first and falls back to this only when it is nothing, as something(model.from, infectious_from(progression)).

EpiBranch.make_contact! Function

julia

make_contact!(
    state::SimulationState,
    parent::Individual,
    infection_time::Real;
    type_idx
) -> Individual

julia

make_contact!(state, parent, infection_time; type_idx = NoTypeLabels())

Add one new contact of the case parent, exposed at infection_time (days), to the outbreak, and return it. Its population characteristics are drawn when it is created; interventions act on it later in the generation.

Only needed when writing a new transmission model. The simulation calls it for every contact drawn from a branching process's offspring distribution, or from generate_offspring in a model that defines it. A model that creates new people inside its contacts_of calls it directly, returning each contact with its infection time. Such a model also needs a collect_exposures method that returns gather_by_target: without one the simulation asks it for generate_offspring and never reaches contacts_of. A model whose contacts are all new and share one generation time is simpler written as generate_offspring; contacts_of earns its keep where the model times each contact itself, as here, where every route out of a case has its own delay:

julia

function EpiBranch.collect_exposures(m::MyModel, state::SimulationState)
    return gather_by_target(m, state)
end

function EpiBranch.contacts_of(m::MyModel, parent, state::SimulationState)
    map(m.routes) do route
        t = parent.infection_time + rand(state.rng, route.delay)
        (make_contact!(state, parent, t), t)
    end
end

The simulation does everything else: it applies interventions to the case and its contacts, decides whether each contact is infected, draws the natural history of new cases and updates the case counts.

EpiBranch.susceptible_fraction Function

julia

susceptible_fraction(
    state::SimulationState{<:Any, <:Any, NoPopulation}
) -> Float64
susceptible_fraction(
    state::SimulationState{<:Any, <:Any, NoPopulation},
    extra_infected::Int64
) -> Float64

julia

susceptible_fraction(state::SimulationState, extra_infected::Int = 0) -> Float64

Share of the population still susceptible at this point of the outbreak: 1.0 when the population is unbounded (no population_size), otherwise one minus the share already infected. extra_infected counts people infected in the current generation who are not yet in the case count.

For extension authors: a package can add its own kind of structured population by defining a new type for population_size and a method of this function for it.

EpiBranch.initialise_state Function

julia

initialise_state(
    model::TransmissionModel,
    sim_opts::SimOpts,
    interventions,
    transitions,
    attributes,
    rng::Random.AbstractRNG
) -> SimulationState

julia

initialise_state(model, sim_opts, interventions, transitions, attributes, rng) -> SimulationState

Set up the outbreak at time 0, before any transmission: the population and its index cases, as a SimulationState. The default, used by branching processes, creates sim_opts.n_initial index cases. A model with a fixed population that is used up as people are infected (a network, households) defines its own method, usually by building the population with new_state and add_individuals! and infecting the index cases with seed!.

EpiBranch.new_state Function

julia

new_state(
    model::TransmissionModel,
    transitions,
    attributes,
    rng::Random.AbstractRNG
) -> SimulationState

julia

new_state(model, transitions, attributes, rng) -> SimulationState

An empty outbreak for model, before anyone exists or is infected: a SimulationState at generation 0 holding the model's population size, the population characteristics (attributes) and the natural-history transitions. Used when writing a new transmission model: its initialise_state starts from this and adds people with add_individuals! and index cases with seed!.

EpiBranch.add_individuals! Function

julia

add_individuals!(
    state::SimulationState,
    n::Integer,
    interventions;
    n_types,
    setup,
    infection_time
) -> Union{Vector{Any}, Vector{E} where E<:Individual}

julia

add_individuals!(state, n, interventions; n_types = 1, setup = (ind, i) -> nothing,
                 infection_time = NaN)

Add n people to the outbreak's population and return them. Used by a new transmission model's initialise_state to build its population.

Each person starts with the given infection_time: NaN by default, since nobody is infected yet. seed! sets the index cases' infection time, and the model sets everyone else's if they are infected later. Pass infection_time = 0 when everyone added is an index case, so that population characteristics drawn at creation see the time of their infection.

For each person, the population characteristics are drawn first, then setup(ind, i) runs (to record model-specific information such as a network node or household id), then a random type is drawn for multi-type models, and finally each intervention sets up its own information on that person.

EpiBranch.seed! Function

julia

seed!(
    state::SimulationState,
    ids,
    interventions,
    transitions
) -> SimulationState

julia

seed!(state, ids, interventions, transitions) -> state

Make the people with the given ids the index cases of the outbreak, infected at time 0. Used by a new transmission model's initialise_state after add_individuals!.

It sets their infection time to 0, marks them infected, sets their symptom onset from any incubation period, and sets the outbreak's case count and list of infectious cases. It also checks, on the first index case, that the information the interventions and natural history need is present, so a missing population characteristic gives a clear error.

EpiBranch.close_episode! Function

julia

close_episode!(
    ind::Individual,
    episode::InfectionEpisode
) -> Individual

julia

close_episode!(ind::Individual, episode::InfectionEpisode)

Keep a past infection of ind in ind.episodes when they are reinfected. episode is usually recorded with InfectionEpisode before the new infection replaces ind's current fields. ind.secondary_case_ids is emptied, so the new infection counts only its own secondary cases.

Stopping rules

EpiBranch.should_stop Function

julia

should_stop(
    _::AbstractStoppingRule,
    _::SimulationState
) -> Bool

julia

should_stop(rule::AbstractStoppingRule, state::SimulationState) -> Bool

Whether rule ends the simulation at this point of the outbreak. Define a method of this function for a new stopping rule. Default: false.

EpiBranch.time_bound Function

julia

time_bound(_::AbstractStoppingRule) -> Float64

julia

time_bound(rule::AbstractStoppingRule) -> Real

The time limit (days) a stopping rule sets, or Inf if it sets none (the default). The mass-action, network and household models end the whole run at the earliest such limit among its rules, before any candidate infection past it is drawn. Define it alongside should_stop for a rule that, like MaxTime, should also end those runs.

EpiBranch.generation_bound Function

julia

generation_bound(_::AbstractStoppingRule) -> Int64

julia

generation_bound(rule::AbstractStoppingRule) -> Int

The generation limit a stopping rule sets, or typemax(Int) if it sets none (the default). The mass-action, network and household models settle cases in time order, not generation order, so state.current_generation cannot be trusted to climb monotonically as should_stop is polled; those models instead drop a settled case's own onward proposals once its generation reaches this bound. Define it alongside should_stop for a rule that, like MaxGenerations, should also cap such a run's depth.

EpiBranch.honoured_without_should_stop Function

julia

honoured_without_should_stop(
    _::AbstractStoppingRule
) -> Bool

julia

honoured_without_should_stop(rule::AbstractStoppingRule) -> Bool

Whether the mass-action, network and household models apply rule in full through time_bound and generation_bound alone, with no need to poll should_stop after each settled case. True for Extinction, MaxTime and MaxGenerations. Any other rule keeps the default false, and those models poll should_stop for it instead — including a rule that bounds time or generation depth but also tests something else, such as a case count, since only the bound is captured without asking.

Writing an intervention

EpiBranch.initialise_individual! Function

julia

initialise_individual!(
    _::AbstractIntervention,
    individual,
    state
)

Record this intervention's starting information on a person when they enter the simulation (for example "not yet isolated"). Default: does nothing.

Also called again, before resolve_individual!, on a case settling a new episode after a reinfection, so per-episode state it draws (such as Isolation's :test_positive) is fresh for that episode rather than left over from an earlier one.

julia

initialise_individual!(
    v::AbstractVaccination,
    individual,
    state
)

Record a person as not yet vaccinated with this dose, unless the model's attributes already gave them the dose (for example from a campaign before the simulation starts), in which case that dose and its date are kept.

For an all-or-nothing vaccine, such an earlier dose's efficacy is turned into full or no protection here, once, as for a dose given during the simulation. A recorded efficacy of 0 or 1 is kept as it is.

julia

initialise_individual!(
    _::AbstractClinicalTransition,
    individual,
    state
)

Set the step's default state on a newly created case. Does nothing by default.

EpiBranch.resolve_individual! Function

julia

resolve_individual!(
    _::AbstractIntervention,
    individual,
    state
)

julia

resolve_individual!(intervention, case, state)

Decide what the intervention does to a case before it transmits, for example when it is isolated. Default: does nothing.

Both engines call it once for each infection, as the case is settled: after its infection time is fixed and its clinical transitions have run, and before it infects anyone. It is never called on an uninfected individual, even a contact that keep_active keeps in the simulation to generate contacts of its own. While it runs, state.max_infection_time is the case's infection time, state.current_generation is the case's generation and state.cumulative_cases counts the case and every case settled before it. A Scheduled gate is evaluated against these values.

The continuous-time models settle cases in order of infection time. The generation engine settles a generation's cases together, once its infections are decided. There state.cumulative_cases counts the cases of earlier generations and those of the case's own generation infected no later than it; a case of a later generation is not counted, even one infected earlier.

julia

resolve_individual!(
    _::AbstractClinicalTransition,
    individual,
    state
)

Draw whether and when this step happens to the case, and record it. Does nothing by default.

EpiBranch.apply_post_transmission! Function

julia

apply_post_transmission!(
    iv::AbstractIntervention,
    state,
    new_contacts
)

julia

apply_post_transmission!(intervention, state, new_contacts)

Deprecated: implement trace_contacts! for what an intervention does to the contacts a case reached, and intervention_actions for actions with a delivery time. Both engines call those once per case.

The generation engine used to hand this hook a whole generation's new contacts at once. It still does for an intervention that defines a method of its own, with a deprecation warning, so an intervention written against it keeps working there. Called on an intervention without one, it traces each contact from its parent_id and then admits the intervention's actions for the batch, as the built-in methods it replaces did.

EpiBranch.on_infection_settled! Function

julia

on_infection_settled!(
    _::AbstractIntervention,
    individual,
    state,
    rng
)

julia

on_infection_settled!(intervention, individual, state, rng)

Continuous-time models only: called on each case as soon as its infection time is fixed, before its onset, clinical course or other interventions are worked out from it.

Use it for an effect that depends on when the person was actually infected. For example, RingVaccination can vaccinate someone before they are infected; whether post-exposure protection stops the infection is decided here, once the infection time is known. Draw random numbers from the rng passed in, not from state. Default: does nothing.

EpiBranch.trace_contacts! Function

julia

trace_contacts!(
    _::AbstractIntervention,
    state,
    infector,
    contacts
)

julia

trace_contacts!(intervention, state, infector, contacts[, not_before])

Act on the contacts that infector reached and that it can still affect. Default: no-op.

Both engines call it once per case, after the case is settled and before any of these contacts is decided infected or not, and before any intervention's intervention_actions for the same case. While it runs, state.max_infection_time is the infector's infection time, as in resolve_individual!. The continuous-time models call it as the race settles the case, with its neighbours or household members not yet settled. The generation engine calls it for each case of the previous generation, with the contacts it made in this one. A contact several cases reached is passed to each of them. The generation engine also calls it on an uninfected individual that keep_active kept generating contacts, with the state as the generation left it.

not_before[i], when given, is the earliest time contacts[i] can be sought: when that person became a contact of infector. A contact met only at a funeral cannot be traced before the funeral, while a household member is a contact from the start and has not_before = -Inf. On RoutedNetwork it is the time of the linking route's contacts_from state. The race calls this method whenever the model yields (id, time) contacts and the four-argument method otherwise; by default the five-argument method calls the four-argument one, so an intervention that does not time its action from the contact need not handle it.

Contacts are collected only once some intervention in the stack returns true from traces_contacts. Once they are, this method is called on every intervention, so declaring false keeps the walk from starting but does not keep this method dormant in a stack that holds something which does trace.

Both engines also call it a second time per case, for an intervention declaring traces_source: once with infector the case itself and contacts its own infector, the single-element collection (source,), so a backward-tracing intervention such as SourceTracing can act on the person who infected this case rather than the people it infected. An index case or a community introduction has no infector, so this second call is skipped for it.

A backward-tracing intervention still receives the first, forward call too, with case's real contacts: declaring traces_source only adds the second call, it does not withdraw the intervention from the first. Check that the contacts offered actually is case's own infector, for example only(contacts).id == case.parent_id, before acting, so that call is a no-op on a stack's ordinary forward round; SourceTracing does this.

julia

trace_contacts!(
    ct::ContactTracing,
    state,
    infector,
    contacts
)
trace_contacts!(
    ct::ContactTracing,
    state,
    infector,
    contacts,
    not_before
)

Trace the contacts of infector, applying the same eligibility, probability, delay and action to each pair. Each contact is traced no earlier than its not_before time when one is given. Both engines call this once per case.

julia

trace_contacts!(st::SourceTracing, state, case, contacts)
trace_contacts!(
    st::SourceTracing,
    state,
    case,
    contacts,
    not_before
)

Trace case back to its own infector, named in contacts as both engines offer it (see traces_source): the same eligibility, probability, delay and action ContactTracing would use for a forward trace, applied instead to the single contact offered here.

EpiBranch.traces_contacts Function

julia

traces_contacts(iv::AbstractIntervention) -> Bool

julia

traces_contacts(intervention) -> Bool

Whether intervention implements trace_contacts!. The continuous-time race gathers a case's reachable contacts only when some intervention says true, so an intervention that does not trace costs nothing. The default reads the methods themselves: true for a type with a trace_contacts! method of its own, so an intervention written outside the package need not declare it. A type with a trace_contacts! method that should not count, such as one that only opts out of tracing as CapacityConstrained's does, defines traces_contacts(::T) = false itself.

EpiBranch.traces_source Function

julia

traces_source(_::AbstractIntervention) -> Bool

julia

traces_source(intervention) -> Bool

Whether intervention traces a case back to its own infector rather than forward to the people it infected: both engines then call its trace_contacts! a second time per case, as documented there. Default: false, so writing trace_contacts! for forward tracing (as ContactTracing does) never also triggers a backward call; a backward-tracing intervention such as SourceTracing declares traces_source(::T) = true itself. Unlike traces_contacts, this is not inferred from the methods an intervention defines, since both directions are written through the same trace_contacts! verb.

EpiBranch.has_effect Function

julia

has_effect(iv::AbstractIntervention) -> Bool

julia

has_effect(intervention) -> Bool

Whether intervention can change a simulation or a likelihood. ModelSpec asks this of each intervention it is given and warns when one answers false, so an inert intervention does not run unnoticed.

The default reads the method tables: true for a type with a method of its own, at the number of arguments the engine calls it with, for one of the hooks the package dispatches on an intervention. An intervention that reaches the model some other way, through a seam a package of its own defines for instance, defines EpiBranch.has_effect(::MyIntervention) = true to say so. The question is asked once per intervention type, so the answer should not depend on the instance.

EpiBranch.supplies_contacts Function

julia

supplies_contacts(_::TransmissionModel) -> Bool

julia

supplies_contacts(model) -> Bool

Whether a continuous-time model knows each case's contacts (household members or network neighbours), so that contact tracing has people to act on (see trace_contacts!). True for the household and network processes; false by default, and for the mass-action models (HomogeneousProcess, MixingProcess), whose random mixing has no list of contacts. A model returning true must pass a contacts function to _sellke_race!.

EpiBranch.keep_active Function

julia

keep_active(
    _::AbstractIntervention,
    state,
    targets,
    is_new
) -> Union{Tuple{}, Vector{Int64}}

julia

keep_active(intervention, state, targets, is_new) -> iterable of Int

The ids of the targets a case reached that a ring of contact tracing grows from, or that should go on generating contacts of their own although they are not infected. Default: none.

Every built-in model asks it during tracing, straight after trace_contacts! for each case or ring member the ring grows from, with the contacts that call received: return those the ring is traced onward from. On a branching process their own contacts are then drawn, if they can be before the person settles (see reads_settled_state). It is not asked when no intervention traces.

A model written against generate_offspring or contacts_of asks it instead once per case, right after the case's trace_contacts! and intervention_actions, with the same contacts and the same state. An intervention still implementing the deprecated apply_post_transmission! is asked there once per generation, about the whole batch, after that hook has run. The contacts it returns stay active for one more generation, so their contacts are offered to them through the same hooks in turn, and the InfectiousSource risk source (a default) stops them infecting anyone.

is_new[i] flags which of targets were freshly created, on a model run one generation at a time; every built-in model passes false for all of them. A contact returned here that turns out to be infected is active as a case anyway.

julia

keep_active(
    ct::ContactTracing,
    state,
    targets,
    is_new
) -> Union{Tuple{}, Vector{Int64}}

With depth > 1, keep following uninfected traced contacts for one more generation so tracing can reach their contacts. They do not infect those contacts. Returns the ids of those uninfected ring members (none for depth == 1).

EpiBranch.competing_risk Function

julia

competing_risk(
    _::AbstractIntervention,
    parent,
    contact,
    state
) -> Union{Nothing, EpiBranch._StretchRisks{Risk, Vector{Any}, Float64}, EpiBranch._StretchRisks{Risk{_A, Float64, _B}, Array{Tuple{T, T}, 1}, Float64} where {_A, _B, T}}

julia

competing_risk(intervention, parent, contact, state)
    -> Union{Nothing, Risk, collection of Risk}

How this intervention can prevent the infection of contact by parent (the infector): one or more Risks giving the probability and the time window, or nothing if it has no effect on this infection. Default: nothing.

An intervention that acts in more than one way, such as ring vaccination reducing both the contact's susceptibility and the infector's onward transmission, returns several risks, each applied independently: a tuple of them, or any other collection that can be iterated and indexed, as the built-in removals return a vector with one risk per stretch of removal.

A community introduction on a continuous-time model has no infector, and the person being introduced stands in for one, so a risk that reads the infector sees the contact itself. Return nothing when parent === contact if that is not what your risk means, as RingVaccination's onward-transmission risk does.

Both engines resolve risks after the infector's trace_contacts! and intervention_actions have run, so a risk can read state that other interventions have written on the contact (e.g. :vaccination_time set by tracing-driven vaccination). The generation engine resolves a generation's contacts once every case of the previous generation has been through those hooks; the continuous-time models resolve each infection as they propose it.

julia

competing_risk(
    iso::Isolation,
    parent,
    contact,
    state
) -> Union{Nothing, EpiBranch._StretchRisks{Risk, Vector{Any}, Float64}, EpiBranch._StretchRisks{Risk{_A, Float64, _B}, Array{Tuple{T, T}, 1}, Float64} where {_A, _B, T}}

An isolated infector does not infect a contact while isolated, from its isolation time until its isolation_release_time (never released when duration is Inf). Each such contact is prevented with probability 1 - post_isolation_transmission.

EpiBranch.Risk Type

julia

struct Risk{T, P, R}

julia

Risk(event_time, block_probability, release_time)

How an intervention can prevent one infection: from event_time until release_time (days), an infection that would happen is prevented with probability block_probability. Every risk in force at the time of the infection is drawn independently, and the contact is infected only if none of them prevents it.

Each field is a number or a function (rng, parent, contact, state) -> Real, for example vaccine efficacy that depends on age, or a quarantine's duration. parent is the infector.

Leave event_time = -Inf (the default) for protection that does not start at a particular time, such as reduced susceptibility. Leave release_time = Inf (the default) for protection that never lapses, such as isolation or quarantine with an infinite duration.

Returned by competing_risk.

Fields

event_time::Any
block_probability::Any
release_time::Any

EpiBranch.HostSusceptibility Type

julia

struct HostSusceptibility

Partial susceptibility of the exposed person: a contact whose susceptibility is below 1 escapes infection with probability 1 - susceptibility. Applied to every exposure in a branching process, like an intervention's competing_risk.

Fields

EpiBranch.InfectorInfectiousness Type

julia

struct InfectorInfectiousness

Reduced infectiousness of the infector: a case whose infectiousness is below 1 fails to infect each contact with probability 1 - infectiousness. Applied to every exposure in a branching process, like an intervention's competing_risk.

Fields

EpiBranch.InfectiousSource Type

julia

struct InfectiousSource

Only infected people transmit. A contact who was not infected can still be followed up for its own contacts (for example to trace contacts of contacts, see keep_active), but none of those contacts is infected through them. Has no effect when everyone being followed is infected.

Fields

EpiBranch.AbortedInfection Type

julia

struct AbortedInfection

Stopped infections: a case whose infection was stopped before symptom onset (abort_infection!, as post-exposure vaccination with RingVaccination does) infects nobody from that time on, whether or not the intervention that stopped it is still active.

Fields

EpiBranch.HostImmunity Type

julia

struct HostImmunity

Immunity after infection: a person who has been infected cannot be infected again until susceptible_again_time, which is Inf (never) unless a natural-history step sets it. On the generation engine this matters only for a model that exposes an infected person again through its own contacts_of; the continuous-time models (HomogeneousProcess, MixingProcess, HouseholdProcess, the network processes) apply it to every settled member on their own, once their progression can reach :susceptible_again, with no change needed to how they name contacts. Either way, a model allows reinfection by adding a step to its progression that sets :susceptible_again_time; the earlier infection is then kept as a past episode (see close_episode!) when the person is reinfected.

Fields

EpiBranch.infectious_removal_time Function

julia

infectious_removal_time(
    _::AbstractIntervention,
    _::Individual
) -> Any

julia

infectious_removal_time(intervention, individual) -> Real

The time (days) at which this intervention ends individual's infectious period for good. Every built-in model, the BranchingProcess as much as HomogeneousProcess, MixingProcess and the network and household models, ends a case's infectious period at the earliest such time across interventions: Isolation at the isolation time (perfect isolation only, post_isolation_transmission = 0) and ContactTracing at a quarantined contact's trace time, when the isolation or quarantine has no end (duration = Inf). The default Inf (no removal) suits an intervention that reduces transmission per contact, such as a leaky vaccine, which acts through competing_risk instead. A model written against generate_offspring or contacts_of does not read it.

julia

infectious_removal_time(
    ct::ContactTracing,
    ind::Individual
) -> Any

A quarantined contact stops transmitting from its quarantine time. A contact recorded with FlagOnly is removed only once Isolation isolates it. A quarantine that is never released ends the person's infectious period. With a finite duration (see Quarantine) each of their contacts during quarantine is blocked and they transmit again once released, as with Isolation.

EpiBranch.risk_applies Function

julia

risk_applies(_::AbstractIntervention, route) -> Any

julia

risk_applies(intervention, route) -> Bool

Whether this intervention's competing_risk acts on infections through a given transmission route. route is a RouteWindow, or nothing for an infection from outside the population. A model with a single route names it :transmission. Each Infectiousness window of a BranchingProcess is a route named :transmission, whose until holds the window's own states and EpiBranch.INTERVENTION_REMOVAL; it names no from. The default is true, so protection such as vaccination follows a person across routes. Scheduled and CapacityConstrained answer for the intervention they contain.

Isolation and ContactTracing act only on routes listing EpiBranch.INTERVENTION_REMOVAL in until, and not on infections from outside the population. Other interventions may select routes by any property of the route. This does not affect protection set by the model itself, or a model written against generate_offspring or contacts_of, which applies every risk to every contact.

EpiBranch.standing_block Function

julia

standing_block(source) -> Bool

julia

standing_block(source) -> Bool

For intervention authors: whether a block that is certain for an infector and contact (probability 1) stays in force for good once it applies, for example permanent immunity. When it does, a continuous-time simulation (households, networks, HomogeneousProcess, MixingProcess) stops proposing further contacts between the pair. false by default.

The Risk a source returns cannot say this itself. competing_risk reads the current state, so a block that is certain at one contact may have lifted by the next (a ward that reopens, a campaign that ends, a quarantine that expires), and a Risk holding plain numbers looks the same in both cases. Declare true only for a source whose certain block, once in force for a pair, never lifts. A simulation in which a case's infectious period may never end (no until state that every case reaches) needs that declaration to finish; without it, a certain block raises an ArgumentError instead of silently dropping transmission that could still happen.

EpiBranch.intervention_time Function

julia

intervention_time(
    _::AbstractIntervention,
    _::Individual
) -> Any

julia

intervention_time(intervention, individual)

Time (days) at which this intervention acts on a person, such as their isolation time. Scheduled runs an intervention that defines this whether or not its schedule is active at the case's or contact's own clock, checks the schedule again at intervention_time, and undoes the effect with reset! if that falls outside it.

Default: -Inf. An intervention relying on the default is instead gated before it runs: with no named time, reset! has nothing to undo afterwards.

EpiBranch.reset! Function

julia

reset!(_::AbstractIntervention, _::Individual)

julia

reset!(intervention, individual)

Undo this intervention's effect on a person. Scheduled calls it when the effect would have happened before start_time. Default: does nothing.

EpiBranch.abort_infection! Function

julia

abort_infection!(ind::Individual, time::Real)

julia

abort_infection!(ind, time)

Stop the person's infection at time (days), before symptom onset, as post-exposure prophylaxis does; RingVaccination uses it for post-exposure vaccination. For intervention authors: call it from an intervention once the person has an infection time. If it is called more than once, the earliest time counts.

In every transmission model the person stays a case but infects nobody from time on, even after the intervention that stopped the infection is no longer active, and all their transmission routes end then. They have no symptom onset (:onset_time is NaN, :asymptomatic stays false), so nothing that starts at onset happens. Natural-history steps that would take effect at or after time are undone (see resolve_transitions!).

In a branching process, an intervention acting on contacts before it is decided whether they are infected (in trace_contacts!, or in an action it has admitted) sees each contact's earliest exposure time as its provisional infection time. If the contact turns out not to be infected, or to be infected only at or after time, the stop has no effect and the onset is restored.

Raises an error unless time is after the infection time and, for a person with an incubation period, before symptom onset.

EpiBranch.infection_aborted_time Function

julia

infection_aborted_time(ind::Individual{T}) -> Any

julia

infection_aborted_time(ind)

Time at which the person's infection was stopped before symptom onset, for example by post-exposure vaccination, in days since the start of the outbreak (Inf if it was not). Set by abort_infection!.

EpiBranch.set_isolated! Function

julia

set_isolated!(
    ind::Individual,
    time::Real;
    release_time,
    key
)

Record that a person is isolated from time (days) until release_time. Used by interventions that isolate or quarantine people. release_time is required: pass Inf for isolation that never ends, so that indefinite isolation is always a stated choice.

A transmission route that isolation should end lists EpiBranch.INTERVENTION_REMOVAL in its until, which also allows for leaky isolation. :isolated in an until instead refers to a Transition(:isolated, …) step in the natural history.

key names the record_removal! key the stretch goes under, the shared key by default. A removal that only blocks a fraction of the hazard passes its own key instead, apart from the shared one, so a likelihood never mistakes one for the other.

EpiBranch.clear_isolated! Function

julia

clear_isolated!(ind::Individual)

Undo a person's isolation, reversing set_isolated!.

Isolation and contact tracing rules

EpiBranch.is_eligible_for_isolation Function

julia

is_eligible_for_isolation(
    _::IsolationEligibility,
    individual,
    state
) -> Bool

julia

is_eligible_for_isolation(eligibility, individual, state) -> Bool

Whether this case can be isolated under the eligibility rule (true by default).

EpiBranch.records_isolation Function

julia

records_isolation(
    _::IsolationEligibility,
    individual,
    state,
    isolation_time
) -> Any

julia

records_isolation(eligibility, individual, state, isolation_time) -> Bool

Whether an isolation at isolation_time counts as a detection of the case. The case stops transmitting from isolation_time either way; this decides only whether is_isolated reports it, and so whether tracing, group vaccination, line lists and detection counts see it. By default an isolation at or after the case's own outcome_time (recovery or death) is not counted, since isolating someone who has already recovered or died is not a detection.

Define a method to count a late detection anyway. Post-mortem detection is the usual reason, as with an Ebola death found at burial, which triggers tracing:

julia

struct DetectAfterOutcome <: EpiBranch.IsolationEligibility end
EpiBranch.is_eligible_for_isolation(::DetectAfterOutcome, ind, state) =
    !is_asymptomatic(ind)
EpiBranch.records_isolation(::DetectAfterOutcome, ind, state, t) = true

EpiBranch.is_eligible Function

julia

is_eligible(
    _::TraceEligibility,
    infector,
    contact,
    state
) -> Bool

julia

is_eligible(eligibility, infector, contact, state) -> Bool

Whether contact of the case infector is traced under this rule (true by default).

EpiBranch.trigger_time Function

julia

trigger_time(_::TraceEligibility, infector, state) -> Any

julia

trigger_time(eligibility, infector, contact, state) -> Float64
trigger_time(eligibility, infector, state) -> Float64

The time (in days) from which contacts of the case infector are traced under this eligibility rule; ContactTracing adds its isolation_to_trace_delay on top. By default this is the case's isolation time, or Inf (never) if the isolation does not count as a detection (see is_isolated). OnSymptomOnset uses the onset time instead, so suspicion-based tracing starts at symptom onset rather than waiting for isolation or confirmation.

Combined rules take their time from the conditions they combine. a & b starts at the latest time among its conditions, so OnSymptomOnset() & !OnIsolation() traces from onset; a | b starts at the earliest time among the conditions that hold. A negation that holds, such as !OnIsolation() for a case not yet isolated, has no time of its own: it does not move the time of an &, and a rule with no timed condition starts at the case's isolation time (or earlier if a timed branch of an | holds). A combined rule that does not hold gives Inf.

Rules that are logically equal usually give the same time, but not always, because untimed negations and TraceEveryone() (timed at isolation) can add or remove times. For example, S & TraceEveryone() with S = OnSymptomOnset() starts at isolation while S alone starts at onset.

For extension authors: a rule timed from another event of the case defines a three-argument method; one timed from the contact defines a four-argument method (the four-argument form falls back to the three-argument one). ContactTracing calls the four-argument form, and only once is_eligible holds. A single rule's trigger_time does not itself check eligibility. A combined rule containing a condition that reads the contact must be timed with the four-argument form.

EpiBranch.traces Function

julia

traces(_::TraceRate, infector, contact, state, rng) -> Any

julia

traces(rate, infector, contact, state, rng) -> Bool

Whether this eligible contact is successfully traced (randomly, using rng).

EpiBranch.draw_trace_delay Function

julia

draw_trace_delay(
    _::TraceDelay,
    infector,
    contact,
    state,
    rng
) -> Any

julia

draw_trace_delay(delay, infector, contact, state, rng) -> Float64

Draw the delay in days before this contact is reached.

EpiBranch.apply_trace! Function

julia

apply_trace!(
    _::TraceAction,
    contact,
    state,
    trace_time,
    rng
)

julia

apply_trace!(action, contact, state, trace_time, rng)

Act on contact, traced at trace_time (days); changes contact in place.

Vaccination

EpiBranch.vaccine_effect Function

julia

vaccine_effect(v::AbstractVaccination) -> VaccineEffect

The VaccineEffect of a vaccination: what its dose does. Every vaccination intervention defines this, and the package reads efficacy, delay, mode and dose label through it.

EpiBranch.realised_efficacy Function

julia

realised_efficacy(_::LeakyMode, eff, rng) -> Any

julia

realised_efficacy(mode::AbstractEffectMode, eff, rng) -> Real

Turn a vaccinated person's drawn efficacy eff into the probability that each of their exposures is prevented once immune. Called once, when the dose is given; the result applies at every exposure. LeakyMode keeps eff. AllOrNothingMode draws once whether the person is protected, giving 1.0 (fully protected) with probability eff and 0.0 (no protection) otherwise. A new mode defines this method.

EpiBranch.realise_prior_dose! Function

julia

realise_prior_dose!(
    mode::AbstractEffectMode,
    individual,
    label,
    state
)

julia

realise_prior_dose!(mode::AbstractEffectMode, individual, label, state)

For a dose the model's attributes gave a person before the simulation started, apply realised_efficacy to its recorded efficacy, so it gets the same per-person draw as a dose given during the simulation. Changes individual.state in place.

A new mode usually needs no method of its own: the default applies the mode's realised_efficacy to the recorded efficacy, which is what AllOrNothingMode needs. LeakyMode does nothing, since its recorded efficacy is already the value it uses.

EpiBranch.supports_waning Function

julia

supports_waning(_::AbstractEffectMode) -> Bool

julia

supports_waning(mode::AbstractEffectMode) -> Bool

Whether waning can be used with this vaccine mode. Waning reduces the per-exposure protection of a leaky vaccine; an all-or-nothing vaccine has no such protection to reduce, so AllOrNothingMode returns false and VaccineEffect rejects waning with it. Default: true. A new all-or-nothing style mode returns false to get the same check.

Interventions with limited resources

EpiBranch.InterventionAction Type

julia

struct InterventionAction{I, T, F}

julia

InterventionAction(individual, time, effect!)

A proposed action of an intervention (such as giving one dose) on individual at day time, which Scheduled and CapacityConstrained can accept or turn down before it happens.

For extension authors: effect!(individual, time, state) performs the action once accepted; proposing an action must not record it. The effect may record a future date, but acceptance happens on the current simulation clock. An effect limited by capacity sets its intervention's capacity_key only when it uses up capacity.

Fields

individual::Any
time::Any
effect!::Any

EpiBranch.intervention_actions Function

julia

intervention_actions(
    _::AbstractIntervention,
    state,
    candidates
) -> Any

julia

intervention_actions(intervention, state, candidates)

The actions (such as doses) an intervention proposes for candidates, as InterventionActions, so that each one can be checked before it happens. An intervention may widen the candidates, for example to everyone in their groups. Return nothing, the default, when the intervention has no actions.

Both engines call it once per case, after every intervention's trace_contacts! for that case, so a producer can read what tracing wrote. The race, which runs every built-in model, passes the case and the contacts it can still reach, or on a fixed-population model everyone not yet settled for an intervention declaring continuous_actions, and admits only actions dated no earlier than the case's infection time. A branching process also calls it each time it creates a contact, with that contact alone as the candidate, for an intervention declaring needs_contacts = false, and admits that contact's actions at any finite time, since nothing about the new person is settled yet. The generation engine, which runs a model written against generate_offspring or contacts_of, passes the contacts the case made and admits every action with a finite time; once a generation's cases are settled, it also calls it again for each settled case in turn, with that case alone as the candidate, for an intervention that declares settles_generation_batch, so a trigger that only becomes true on settlement (such as a case's own lab confirmation) still reaches its group. Scheduled and CapacityConstrained admit the actions through their own gates.

EpiBranch.action_draw! Function

julia

action_draw!(sample, individual, key) -> Any

julia

action_draw!(sample, individual, key)

Draw sample() once for this person and action key, and return the same value whenever the action is proposed again, so a person who was turned down (for example by a coverage draw) stays turned down. Use different keys for different visits or doses. The stored draws do not appear in the line list.

EpiBranch._action_cache Function

julia

_action_cache(individual) -> Any

julia

_action_cache(individual)

For extension authors: the per-person store of random draws made when proposing actions (see action_draw!). An intervention can also keep its own records there, such as which proposal a decision was made under. It does not appear in the line list.

EpiBranch.apply_actions! Function

julia

apply_actions!(iv, state, candidates)

julia

apply_actions!(intervention, state, candidates)

Propose the intervention's actions for candidates and perform those that its Scheduled and CapacityConstrained settings accept. Both check each action's own day: a schedule against its start and end, and capacity against the budget of the period that day falls in. A person turned down may be offered the action again later, with the same random draws. The intervention itself leaves out actions already performed.

EpiBranch.continuous_actions Function

julia

continuous_actions(_::AbstractIntervention) -> Any

julia

continuous_actions(intervention) -> Bool

Whether this intervention can act once a case's infection is settled, rather than only when it makes a contact. The race offers every intervention its actions per settled case, and offers one declaring this, on a fixed-population model, everyone not yet settled rather than only the contacts the case can reach. The generation engine, by default (see settles_generation_batch), also offers it the whole batch once a generation is settled. Return true only if proposing actions needs no infection time of a contact not yet infected and never changes a case already settled. The default is false.

EpiBranch.continuous_candidates Function

julia

continuous_candidates(
    iv::AbstractIntervention,
    state,
    current,
    traced,
    is_pending,
    members,
    contacts
) -> Any

julia

continuous_candidates(intervention, state, current, traced, is_pending, members, contacts)

The people to offer intervention's actions to once current's infection is settled. For an intervention declaring continuous_actions, the default is current and every member of members that is_pending still admits, which suits any intervention but costs the size of the race, however few people it actually targets. For any other, the default is current and the contacts it can still reach, those tracing reached from it where there was tracing, as for RingVaccination. GroupVaccination writes its own method, offering the members of current's group instead.

traced is the set tracing reached, or nothing if there was no tracing. A ring wider than one step reaches people who are not current's direct contacts, so a list built from current's contacts alone would miss them. is_pending(id) reports whether id belongs to this race and has not yet been settled, so an action may still reach it. The defaults (walking current's direct contacts when there is no traced to use) need members and contacts; an intervention narrowing its own candidates usually needs only current, traced and is_pending.

Write a method of this function for an intervention declaring continuous_actions whose candidates are not the whole population: an intervention written outside the package otherwise gets the default above, however few people it targets.

EpiBranch.acts_on_population Function

julia

acts_on_population(_::AbstractIntervention) -> Bool

julia

acts_on_population(intervention) -> Bool

Whether this intervention's actions belong to the whole population rather than to a case's contacts or group. On the race, intervention_actions is then called once, before any case settles, with every member as a candidate, and for each contact a branching process creates, instead of per settled case. Return true only if proposing an action needs nothing a case's settlement would supply: no infection time, no exposure, no settled state. The default is false.

MassVaccination declares it: each person's dosing time comes from their own eligibility_time, not from who infected whom, so the whole rollout can be decided upfront. A dose dated later in the simulation still protects correctly, because the per-contact competing risk is checked against the recorded dose time when a candidate infection is popped, not when it was proposed.

EpiBranch.needs_contacts Function

julia

needs_contacts(_::AbstractIntervention) -> Bool

julia

needs_contacts(intervention) -> Bool

Whether this intervention's actions need the model to name each case's contacts. The race narrows intervention_actions's candidates to a case and the people it can still reach, through supplies_contacts; the homogeneous and mixing models name none, so an intervention that needs them is never offered them there. The default is true, matching RingVaccination, whose candidates come from walking the contact graph. GroupVaccination declares false: its candidates are a case's group, found through its group_key, not through any contact. A branching process offers an intervention declaring false each contact as it creates it, too, so a group's dose reaches a member created after it.

EpiBranch.settles_generation_batch Function

julia

settles_generation_batch(iv) -> Any

julia

settles_generation_batch(intervention) -> Bool

Whether _settle_generation! should also offer intervention's actions for the whole settled batch, on top of the contacts each source is already offered mid-generation (_case_hooks!). Defaults to continuous_actions. RingVaccination is the one exception: a contact reaches it, pre-infection, through _case_hooks! already, which is where post_exposure_efficacy needs to reconsider a dose, against an exposure not yet complicated by that contact's own clinical resolution. Offering it again from the settled batch would reconsider a dose after this same round has already resolved the case's onset and isolation from it.

EpiBranch.may_revise Function

julia

may_revise(
    _::AbstractIntervention,
    prior_trigger,
    new_trigger
) -> Any

julia

may_revise(intervention, prior_trigger, new_trigger) -> Bool

Whether a dose already accepted under prior_trigger may be moved to new_trigger, found later. The default is false: an accepted dose keeps its date. GroupVaccination and RingVaccination allow the move when new_trigger is earlier, since in continuous-time models a group's earliest eligible case, or a contact's earliest trace, may only be found after a later one.

EpiBranch.is_settled Function

julia

is_settled(state, ind::Individual) -> Bool

julia

is_settled(state, ind) -> Bool

For intervention authors, in the network and household models: whether ind's infection, and what interventions do to them, is final for this run, so nothing found later will change it. false for a person whose infection is still undecided, and for the case currently being processed (so its own interventions can still be revised).

EpiBranch.persistent_competing_risks Function

julia

persistent_competing_risks(_::AbstractIntervention) -> Bool

julia

persistent_competing_risks(intervention) -> Bool

Whether protection already given (for example to a vaccinated person) continues after a Scheduled campaign has stopped. The default is false; vaccination returns true, and Scheduled and CapacityConstrained answer for the intervention they contain. Scheduled still controls who newly receives the intervention.

For extension authors: an intervention returning true must base its competing_risk on the effects it has recorded and their dates, returning nothing before any effect is recorded.

EpiBranch.capacity_key Function

julia

capacity_key(iv::AbstractIntervention) -> Any

julia

capacity_key(intervention) -> Symbol

For extension authors: the key in ind.state that marks a person as having used one unit of intervention's capacity (for vaccination, received the dose with its dose_label). CapacityConstrained counts it to measure how much capacity has been used. Defined for RingVaccination, MassVaccination and GroupVaccination. Define it, and capacity_time_key, to make a new intervention capacity-limited.

EpiBranch.capacity_time_key Function

julia

capacity_time_key(iv::AbstractIntervention) -> Any

julia

capacity_time_key(intervention) -> Symbol

For extension authors: the key in ind.state holding the day the capacity_key was set (for vaccination, the vaccination day). CapacityConstrained places usage, including a dose given by another intervention sharing the same capacity, in the period this date falls in, so a dose a delay dates later than its admission is charged to the day it is actually given.

Start triggers

EpiBranch.is_triggered! Function

julia

is_triggered!(trigger, state::SimulationState) -> Bool
is_triggered!(trigger, state::SimulationState, action) -> Bool

Whether trigger's condition holds now. A trigger must keep no state of its own between calls, because the same trigger object is reused across the repeated simulations of simulate(model, n) (serially, or across threads under parallel = true): derive what it needs from state, or keep working values in state.scratch, which belongs to the run.

The engine asks the gate about 1.7 times per individual, so a trigger that walks the population on every call makes a run quadratic in the number of individuals. Read a running total instead, as Infections reads cumulative_cases.

The three-argument form is asked when Scheduled admits an InterventionAction (such as a dose), with action the one proposed. Its default calls the two-argument form, ignoring action, so an existing trigger keeps working unchanged. A trigger that reads action.individual, for example a group attribute, gates that action for the person it is for rather than the whole population at once.

EpiBranch.trigger_can_lapse Function

julia

trigger_can_lapse(_::AbstractTrigger) -> Bool

julia

trigger_can_lapse(trigger) -> Bool

Whether trigger's condition, once true, can turn false again. The default is the conservative true: a trigger written outside the package is assumed able to close again until it says otherwise, as an opaque Scheduled predicate already is. Infections and ReportedCases both return false, since neither cumulative_cases nor the number of reports already due can fall.

Answer about the accepted timeline. Admitting a scheduled action asks the gate with that action's own time in place of the simulation clock, which settles whether that action is admitted rather than whether the window closes. The distinction is worth the care: a schedule that can close speaks for none of the isolation releases its intervention recorded; true where false belongs costs a likelihood term rather than buying safety.

EpiBranch.required_transitions Function

julia

required_transitions(
    _::AbstractTrigger
) -> Union{Tuple{}, Tuple{UnionAll}}

julia

required_transitions(trigger) -> Tuple of transition types

The transition types a trigger needs the model's progression to hold. Default: (), nothing required.

A trigger reads Individual.state, and the keys it reads are written by the progression: a model composed without the transition that writes them leaves the gate shut for the whole run. ModelSpec asks this of every start_after trigger and warns when the model does not hold what the trigger names. ReportedCases names Reporting.

A key an attributes function supplies, such as :onset_time from clinical_presentation, cannot be named this way: declare it with required_fields instead.

Writing a transition

EpiBranch.resolve_transitions! Function

julia

resolve_transitions!(state::SimulationState, individual)

julia

resolve_transitions!(state, individual)

Draw a case's natural history: the time of each step in the model's progression (becoming infectious, symptom onset, hospitalisation, outcome) and which outcome it reaches first. The results are recorded on individual.state (:infectious_time, :onset_time, :outcome, :outcome_time and whatever else each step records).

The built-in simulation calls this for every new case. Only a new transmission model that runs its own simulation loop needs to call it, once per case, after the case's population characteristics and intervention information are set. The steps come from the model's progression, placed on the outbreak when it is built with new_state.

If the infection was stopped before symptom onset (see abort_infection!), its natural history ends at that time. Steps that take effect before then stand; any step that would take effect at or after it is undone, and steps that follow from an undone one do not happen. Undoing restores each entry a step set on individual.state, so a transition must record its results by setting entries there rather than by changing a vector or dictionary it finds there.

EpiBranch.transition_time Function

julia

transition_time(
    rng,
    individual,
    start_time,
    delay;
    probability,
    purpose
) -> Any

julia

transition_time(rng, individual, start_time, delay; probability = nothing,
                 purpose = :progression)

Draw the time of a clinical event that follows start_time (days since the start of the simulation) by delay days, or nothing if the event does not happen. For writing new progression steps.

delay is a number, a distribution or a function of the random number generator and the individual, (rng, ind) -> .... probability, if given, is the chance the event happens (a number or such a function). When the starting event never happened (start_time is not finite) the result is nothing.

Each supplied probability uses one random draw, even when it is 0 or 1, and probability = nothing uses none. A case whose starting event never happened uses none either, the anchor being checked first; keep both in mind if reproducing a run's random numbers matters. The caller records the returned time and decides which of several terminal events comes first.

purpose names this step for stream: under CommonRandomNumbers its draws then keep their place for this individual regardless of what else the run draws. The probability gate and the delay draw use separate streams, (purpose, :selected) and (purpose, :delay): reshaping probability alone does not move the delay. Steps sharing a name share a stream, so give a new built-in step one of its own.

EpiBranch.censor_after_outcome! Function

julia

censor_after_outcome!(
    _::AbstractClinicalTransition,
    individual
)

julia

censor_after_outcome!(transition::AbstractClinicalTransition, individual)

Reset transition's own event to "did not occur" if the time it recorded falls after individual's :outcome_time: once a case has died or recovered, nothing else can happen to it. Default: no-op.

Called once the simulation engine has resolved every transition for an individual, for every transition with is_terminal(transition) == false, right after :outcome/:outcome_time are set — a terminal transition's own candidate is exactly what decided the outcome (or lost to an earlier one), so it is never censored this way. The built-ins (Hospitalisation, Reporting, and a non-terminal Transition) reset their own flag to false and time to Inf, which is what transition_loglik then reads as a censored (rather than occurred) event, keeping a simulated individual's state and its likelihood in agreement. A custom non-terminal transition overrides this alongside the flag/time keys its own resolve_individual! writes.

EpiBranch.transition_loglik Function

julia

transition_loglik(
    t::AbstractClinicalTransition,
    individual
) -> Any

julia

transition_loglik(t::AbstractClinicalTransition, individual) -> Float64

The log-likelihood of what a step did to one case: the log-probability that the step happened (or did not), plus the log-density of its delay when it did. progression_loglik calls it once per step per case. It is 0.0 when the case never reached the state the step starts from.

Defined for the built-in steps (Transition, Reporting, Hospitalisation, Death, Recovery). A custom step used with progression_loglik needs its own method that reads back what its resolve_individual! recorded. Its delay must be a Distribution or a fixed number, since a function of the individual has no density.

EpiBranch.transition_term Function

julia

transition_term(
    probability,
    delay,
    ind,
    anchor,
    occurred;
    censor
) -> Any

julia

transition_term(probability, delay, individual, anchor, occurred; censor = infection_aborted_time(individual))

The log-probability that a step happened or did not happen to a case, for use inside transition_loglik; when it happened, add the log-density of its delay.

Use it rather than reading probability directly. It gives the right answer in two cases a custom step would otherwise get wrong: outcomes split with exclusive_probabilities, which contribute the probability of the outcome the case drew, and steps undone because the infection was stopped early (for example by post-exposure prophylaxis), which count as not yet having happened by the time the infection was stopped (infection_aborted_time). A step that does not end the case passes its own censor, which also covers the case's outcome time, so an event that had not happened by the time the case died or recovered counts as censored there too.

Population characteristics and recorded contacts

EpiBranch.GroupAttribute Type

julia

struct GroupAttribute{D}

A population characteristic drawn once per group and shared by every member of that group, such as a household's reporting probability. The group is whatever the individual has under group_key, set by groups or an earlier entry in attributes. Each group's value is drawn when its first member is created and kept for the rest of the outbreak.

Create one with group_attribute.

Fields

key::Symbol
group_key::Symbol
propensity::Any
cache::Dict{Any, Any}

EpiBranch.records_contacts Function

julia

records_contacts(
    _::ContactRecorder,
    parent,
    contact,
    state,
    t
) -> Bool

julia

records_contacts(recorder::ContactRecorder, parent, contact, state, t) -> Bool

Whether a continuous-time model should keep simulating contacts between the infector parent and contact after an intervention has protected contact for good (see EpiBranch.standing_block). Asked on every proposal a risk blocks, not only ones certain to recur, with t the time of that proposal, so a recorder can log each blocked contact; only the answer given for a proposal whose block is certain decides whether the pair's draws continue. Default false.

Returning true for a certain block restricts which models can run: a model whose infectious period never ends would simulate contacts for ever, so it raises an error instead, as described in EpiBranch.standing_block.

Writing an observation model

An observation model defines observe for the closed-form likelihood and apply_observation! for simulation.

EpiBranch.apply_observation! Function

julia

apply_observation!(_::NoObservation, state, rng) -> Any

julia

apply_observation!(obs::ObservationModel, state, rng)

Mark which simulated cases are reported, and when, under an observation model. simulate calls it once the outbreak has been simulated; it changes state in place. Under PerCaseObservation each case records :reported (true/false) and :reporting_time. NoObservation and MinimumSize leave the cases unchanged.

To add a new observation model, define this together with observe.

Pairwise likelihood building blocks

EpiBranch.calendar_multiplier Function

julia

calendar_multiplier(s::Steps, t::Real) -> Real

julia

calendar_multiplier(schedule, t)

The multiplier a calendar schedule applies to transmission on day t (non-negative; 1 leaves it unchanged). PairKernel multiplies its contact-interval hazard by it.

Every schedule defines this. A new schedule is a type with this method and, depending on its calendar_shape, either next_calendar_break (piecewise constant, the default) or a calendar_shape method declaring it smooth. Steps is the piecewise-constant schedule the package provides.

EpiBranch.next_calendar_break Function

julia

next_calendar_break(s::Steps, t::Real) -> Any

julia

next_calendar_break(schedule, t)

The first day strictly after t on which a piecewise-constant schedule changes value, or Inf if it stays constant from t onwards.

Simulation and the likelihood use it to find where each constant stretch of the schedule ends. A smooth schedule has no breaks and does not define it.

EpiBranch.calendar_shape Function

julia

calendar_shape(schedule) -> EpiBranch.SmoothCalendar

julia

calendar_shape(schedule)

Whether a calendar schedule changes in steps or smoothly: PiecewiseConstantCalendar, the default, which needs next_calendar_break, or SmoothCalendar, which a smooth schedule declares for its own type. A seasonal multiplier, for example:

julia

struct DayOfWeekRamp{T <: Real}
    slope::T
end
EpiBranch.calendar_multiplier(s::DayOfWeekRamp, t) = 1 + s.slope * mod(t, 7)
EpiBranch.calendar_shape(::DayOfWeekRamp) = EpiBranch.SmoothCalendar()

Seasonal is the smooth schedule the package provides, and a plain callable t -> multiplier is read as a smooth schedule too, with no calendar_shape method of its own needed.

EpiBranch.PiecewiseConstantCalendar Type

julia

struct PiecewiseConstantCalendar

julia

PiecewiseConstantCalendar()

Marks a calendar schedule as constant between the days next_calendar_break reports, like Steps. This is the default calendar_shape, and its effect on transmission is computed exactly.

Fields

EpiBranch.SmoothCalendar Type

julia

struct SmoothCalendar

julia

SmoothCalendar()

Marks a calendar schedule whose multiplier changes continuously, such as a seasonal curve, as its calendar_shape. Its effect on transmission is computed numerically, so the multiplier should be smooth. It is assumed never to stop transmission for good.

Fields

EpiBranch.pair_kernel Function

julia

pair_kernel(
    k::Distributions.Distribution{Distributions.Univariate, Distributions.Continuous},
    i,
    j,
    infection_time
) -> Distributions.Distribution{Distributions.Univariate, Distributions.Continuous}

julia

pair_kernel(kernel, infector, susceptible, infector_infection_time)
pair_kernel(kernel, infector, susceptible, infector_infection_time, infectious_time)
pair_kernel(kernel, infector, susceptible, infector_infection_time, infectious_time,
    state)

The contact-interval distribution for one infector and one susceptible person. A shared distribution is returned unchanged; a function receives the two people's numbers; a PairKernel receives a PairContext and, when it has a state, each person's record. The five-argument form adds the start of the infector's infectious period, which a PairKernel with a calendar schedule needs. The six-argument form also passes the SimulationState, from which a PairKernel with a state function reads both people's records; simulation must use it, and the shorter forms, meant for the likelihood, throw for such a kernel. Every other kernel returns what the five-argument form does.

EpiBranch.pairwise_surv_loglik_by_component Function

julia

pairwise_surv_loglik_by_component(
    kernel,
    data::InfectionLayer,
    layout::ContactPairsLayout;
    external_hazard,
    susceptibility
) -> Vector

julia

pairwise_surv_loglik_by_component(kernel, data::InfectionLayer, layout::ContactPairsLayout;
                                  external_hazard = 0.0,
                                  susceptibility = nothing) -> Vector{<:Real}
pairwise_surv_loglik_by_component(kernel, data::InfectionLayer;
                                  external_hazard = 0.0,
                                  susceptibility = nothing) -> Vector{<:Real}

pairwise_surv_loglik split by household (or by connected part of a contact network): entry c is the log-likelihood of group c of layout, and sum(pairwise_surv_loglik_by_component(...)) == pairwise_surv_loglik(...). layout.component gives each person's group and layout.ncomponents their number.

A group containing an infected person whom no possible infector can explain gets -Inf, as does the total, but the other groups keep their finite values. A sampler that imputes infection times one household at a time can therefore accept or reject each proposal on its own entry without rebuilding the layout.

The arguments are as in pairwise_surv_loglik; see also PairwiseReduction for other groupings.

EpiBranch.infection_likelihood_compatible Function

julia

infection_likelihood_compatible(component) -> Bool

julia

infection_likelihood_compatible(component) -> Bool

Whether the pairwise likelihood can account for this intervention (or other model component) from the recorded infectious periods alone: true when it changes transmission only by starting or ending people's infectious periods, changes a person's susceptibility only through susceptibility_components, or does not affect infection at all. The default is false, so loglikelihood(data, model) refuses a model with an intervention that has not declared this. A removal that only reduces transmission, such as leaky isolation, is still compatible when it declares removal_block_probability below 1; the likelihood scales its stretches' hazard rather than excluding them. Changes to infectiousness and to the contact interval still need a kernel that includes them, passed to pairwise_surv_loglik directly.

Declaring true is a promise by the intervention's author; nothing checks it. The infection likelihood takes the infection record as given and does not include the probability of clinical outcomes, of who received interventions, of people's characteristics, or of observation.

EpiBranch.susceptibility_components Function

julia

susceptibility_components(effect, host) -> Any

julia

susceptibility_components(effect, host) -> components or nothing

How an intervention such as vaccination changes one person's risk of infection in the pairwise likelihood (the susceptibility keyword of pairwise_surv_loglik). host is a LayerHost: the person's id, infection_time, and the recorded event times under host.state.

It returns nothing, the default, when effect leaves the person's risk unchanged. Otherwise it returns weight => modifier pairs, each modifier a HazardScaling or nothing (no change), with weights that are probabilities summing to one. The person's likelihood contribution is the weighted mixture of the likelihood of their escapes and infection under each modifier. One component describes an effect every exposure shares; several describe an unobserved state drawn once per person that then governs all their exposures together.

For a vaccination's VaccineEffect under LeakyMode this is one component, 1 => HazardScaling(τ, 1 - efficacy) from the person's immunity time τ (with waning, the factor is dt -> 1 - efficacy * waning(dt)). Under AllOrNothingMode it is two: efficacy => HazardScaling(τ, 0.0) for someone the vaccine fully protects and 1 - efficacy => nothing for someone it does not. The immunity time is read from :immunity_time (:immunity_time_<label> for a labelled dose), and a person without one is unchanged. Every AbstractVaccination returns its VaccineEffect's answer, and an InterventionWrapper that of the intervention it wraps.

For a collection, such as a model's interventions, the one non-nothing answer among them is used, and an ArgumentError is raised if more than one changes the same person. To support a new kind of effect, define a method for it and susceptibility_host_times for the event times it reads.

EpiBranch.susceptibility_host_times Function

julia

susceptibility_host_times(component) -> Tuple{Symbol}

julia

susceptibility_host_times(component) -> Tuple of Symbols

The event times (keys of individual.state) that component's susceptibility_components reads, such as a vaccinee's immunity time. household_infections and network_infections record them in the infection record's host_times whenever the model includes component. The default is (). An AbstractVaccination reads its dose's immunity time, and an InterventionWrapper the times of the intervention it wraps.

EpiBranch.binding_release Function

julia

binding_release(component) -> Bool

julia

binding_release(component) -> Bool

For intervention authors: whether the release_time this intervention reports on a Risk can be relied on. A block it says ends on day t is not still in force after t, and a block it reports as never ending has not lifted by the next contact. The default is false, the cautious answer.

A continuous-time simulation (households, networks, HomogeneousProcess, MixingProcess) asks this where a pair could go on having contacts indefinitely, to tell a block that ends the pair's transmission from one it must keep testing contacts against. Without the declaration such a block raises an error instead of silently dropping transmission that could still happen, for the reason standing_block gives.

The built-in isolation and quarantine declare it, because the periods record_removal! records are only ever added to. A Scheduled intervention that can end withdraws the declaration, since ending it lifts a block before the release it reported.

EpiBranch.removal_gap_host_times Function

julia

removal_gap_host_times(component) -> Tuple{Symbol}

julia

removal_gap_host_times(component) -> Tuple of Symbol

For an intervention that removes a case from transmission for a while and then releases them (isolation that ends, quarantine that expires): the individual.state keys under which it records those periods. household_infections and network_infections add them to the infection record's host_times, and the likelihood removes each recorded period from the exposure of everyone the case could have infected. The default is (), for an intervention that never releases anyone.

Both built-in removals record through record_removal!: Isolation uses the reserved :_removal_stretches key for perfect isolation with a duration that can end, and ContactTracing uses the quarantine's own key for any Quarantine, whatever its duration. A period with no release ends the exposure where it starts. An intervention of your own names the key it records under, whose value is a vector of (start, release) pairs. A wrapper that can lift the removal part-way through a period, such as a Scheduled with an end time, names none. The infectious period then ends at the first removal only if that removal never releases. A removal with a release is left out of the record, and in simulation acts contact by contact instead.

EpiBranch.removal_gap_host_times_withdrawn Function

julia

removal_gap_host_times_withdrawn(
    component
) -> Union{Tuple{}, Tuple{Symbol}}

julia

removal_gap_host_times_withdrawn(component) -> Tuple of Symbol

The keys removal_gap_host_times leaves unnamed on purpose, because component can withdraw the removal part-way through a stretch it (or the intervention it wraps) still records through record_removal! with a release that can lapse. The default is (): most components either name every key they record under or never record one with a release at all.

Every key a running simulation ever records with a finite release is checked against the model's removal_gap_host_times and this trait when the infection layer is built; one named by neither raises an error, since nothing would otherwise say whether leaving it out was deliberate or a forgotten key (see the extending guide's removal_gap_host_times section).

EpiBranch.removal_block_probability Function

julia

removal_block_probability(component) -> Float64

julia

removal_block_probability(component) -> Real

The fraction of the hazard component blocks over each stretch it names in removal_gap_host_times: 1 (the default) for a removal that stops transmission outright, such as a quarantine or perfect isolation, so its stretches are taken out of the exposure entirely. Below 1 for one that only reduces it, such as leaky isolation's 1 - post_isolation_transmission: the likelihood scales the hazard by this fraction over the stretch instead of excluding it, and adds the log of the surviving fraction to the hazard of an infection that falls inside one, rather than treating that infector as ruled out.

Declaring a value below 1 is a promise by the intervention's author, as infection_likelihood_compatible is: nothing checks that the stretches the component records actually reduce transmission by this much. Two components naming the same key must agree on this fraction, since the likelihood cannot tell their stretches apart once merged; give each its own key otherwise.

household_offspring (in EpiHouseholds) weights each recorded stretch by this fraction too, when it works out a case's community-infectious person-time.

EpiBranch.record_removal! Function

julia

record_removal!(
    ind::Individual{T},
    start::Real,
    release::Real;
    key
)

julia

record_removal!(ind, start, release; key = :_removal_stretches)

Record that ind could not transmit from start until release (days; Inf if never released), for example while isolated or quarantined. set_isolated! records under the default key, shared by isolation and quarantine.

An intervention of your own that removes people in its own way passes its own key, and names it in removal_gap_host_times. The likelihood then leaves out the same periods of exposure that the simulation blocked, so that simulate and loglikelihood agree. Every key that ever records a stretch with a finite release here is checked against removal_gap_host_times when the infection layer is built, and the likelihood stops with an error for one that no component names or withdraws.

EpiBranch.removal_stretches Function

julia

removal_stretches(ind::Individual) -> Any
removal_stretches(ind::Individual, key::Symbol) -> Any

julia

removal_stretches(ind, key = :_removal_stretches)

Every period during which ind could not transmit, such as isolation or quarantine, as sorted, non-overlapping (start, release) pairs in days; a release of Inf means it never ended. Recorded by record_removal!.

isolation_time and isolation_release_time give only the current isolation. This gives the whole history, which the likelihood needs: for example a person quarantined, released, and isolated again later. Do not change the returned vector.

EpiBranch.HazardScaling Type

julia

struct HazardScaling{S<:Real, F}

julia

HazardScaling(start, factor)

A change in one person's risk of infection from a given day, such as protection from a vaccine: from day start on, every hazard of infection they face, from each possible infector and from the community alike, is multiplied by factor. factor is a non-negative number, or a function dt -> Real of the days since start for protection that changes over time, such as waning. Before start the hazard is unchanged.

One component of what susceptibility_components returns.

Fields

start::Real
factor::Any

EpiBranch.contact_structure Function

julia

contact_structure(data::InfectionLayer) -> Vector{Int64}

julia

contact_structure(data::InfectionLayer)

Who could have infected whom in data, in a form compile_contact_pairs accepts: a membership vector (people sharing a label can all infect one another, as in households) or a list of contacts (contacts[i] lists the people i can infect, as in a directed contact network). An InfectionLayer subtype defines this.

EpiBranch.followup_end Function

julia

followup_end(data::InfectionLayer) -> Any

julia

followup_end(data::InfectionLayer)

The day observation of data stopped. The pairwise likelihood covers infections and exposure up to it and ignores everything after it. The default reads a followup_end field when the InfectionLayer subtype has one, and is Inf otherwise; a subtype that stores it elsewhere defines a method.

EpiBranch.host_times Function

julia

host_times(data::InfectionLayer) -> Any

julia

host_times(data::InfectionLayer)

Each person's recorded event times in data beyond the start and end of their infectious period (such as onset_time), as a named tuple of vectors, with missing for a person without that time. A PairKernel state function or a vaccine effect reads these in the likelihood as it reads individual.state in simulation. The default reads a host_times field when the InfectionLayer subtype has one, and is empty otherwise; a subtype that stores them elsewhere defines a method.

EpiBranch.PairwiseReduction Type

julia

abstract type PairwiseReduction

julia

PairwiseReduction

How the pairwise likelihood is split into groups, for example by household, age stratum or spatial patch, so that each group's log-likelihood is reported separately. A subtype defines ngroups (how many groups) and group (which group a person's terms add into), and nothing else. pairwise_surv_loglik puts everything into one group; pairwise_surv_loglik_by_component groups by household or connected part of the network. A new grouping is written the same way from outside the package and evaluated with pairwise_reduce.

Fields

EpiBranch.ngroups Function

julia

ngroups(reduction::EpiBranch.PairwiseReduction) -> Int64

julia

ngroups(reduction::PairwiseReduction) -> Int

How many groups reduction splits the pairwise likelihood into. A PairwiseReduction subtype defines this.

EpiBranch.group Function

julia

group(reduction::EpiBranch.PairwiseReduction, host) -> Int64

julia

group(reduction::PairwiseReduction, host::Int) -> Int

Which of reduction's 1:ngroups(reduction) groups the likelihood terms of person host add into. A PairwiseReduction subtype defines this.

EpiBranch.pairwise_reduce Function

julia

pairwise_reduce(
    reduction::EpiBranch.PairwiseReduction,
    kernel,
    data::InfectionLayer,
    layout::ContactPairsLayout;
    external_hazard,
    susceptibility
) -> Any

julia

pairwise_reduce(reduction::PairwiseReduction, kernel, data::InfectionLayer,
                 layout::ContactPairsLayout; external_hazard = 0.0,
                 susceptibility = nothing) -> Vector{<:Real}

The pairwise log-likelihood split into the groups of reduction, a PairwiseReduction: returns one log-likelihood per group, ngroups(reduction) in all. pairwise_surv_loglik and pairwise_surv_loglik_by_component are this call with their own groupings; a new grouping (by stratum, by spatial patch) calls it directly with its own PairwiseReduction subtype. The other arguments are as in pairwise_surv_loglik.

EpiHouseholds.condition_mask Function

julia

condition_mask(rule::ConditionOn, data::HouseholdInfections) -> AbstractVector{Bool}

Which person in each household rule conditions on: a vector shaped like data.is_index, true for each household's conditioned member and false elsewhere. One method per rule.

Internal functions

These are documented for people working on the package itself. Three of them are part of the extension surface rather than internals, and are declared public in Julia's sense: EpiBranch.get_generation_time, EpiBranch.population_size and EpiBranch.required_fields. The rest are helpers the package may change at any time.

EpiBranch.get_generation_time Function

julia

get_generation_time(
    gt::Distributions.Distribution,
    individual
) -> Distributions.Distribution

julia

get_generation_time(gt, individual)

The generation time distribution for one case. A distribution is shared by every case. A function of the individual, (ind) -> distribution, gives each case its own, and can read anything recorded about the case, such as its incubation period or other population characteristics. That way the generation time and symptom onset of a case can be linked rather than drawn independently. Use incubation_period to read the incubation period inside such a function.

Called only for a case with at least one contact to time, after the number of contacts is drawn. A case making none is never asked for a distribution, so the function need not handle a case lacking whatever it reads.

EpiBranch._advance_generation! Function

julia

_advance_generation!(
    model::TransmissionModel,
    state::SimulationState,
    interventions::Vector{<:AbstractIntervention}
)
_advance_generation!(
    model::TransmissionModel,
    state::SimulationState,
    interventions::Vector{<:AbstractIntervention},
    intervention_tuple
)

Simulate one generation of transmission. Contacts are made and timed (collect_exposures); tracing and vaccination act on those contacts (_intervene!); then each contact is infected or not, and each new case is settled, which draws its natural history and lets the interventions act on it (_resolve!, _settle_generation!). In a branching process a contact is created together with its infection time; in a fixed population the contacts already exist.

EpiBranch._settle_generation! Function

julia

_settle_generation!(
    state::SimulationState,
    interventions,
    cases,
    cases_before
)

julia

_settle_generation!(state, interventions, cases, cases_before)

Settle a generation's new cases: every case's clinical transitions (resolve_transitions!), then every case's intervention hooks (resolve_individual!), then, once every case's own state is settled, each intervention declaring settles_generation_batch is offered its actions for every case in turn (_offer_continuous_actions!), so a trigger that only becomes true on settlement (such as a case's own lab confirmation) still reaches its group. The generation engine has no time order beyond the generation. Each case, including in the settled-batch offer, sees the clock at its own infection time and, as its case count, the cases_before cases of earlier generations plus those of its own generation infected no later than it. Cases sharing an infection time are counted in the order they appear in cases. The clock, count and generation are restored to the generation's totals once every case has been offered, which is what the stopping rules read.

EpiBranch._offer_continuous_actions! Function

julia

_offer_continuous_actions!(
    state::SimulationState,
    interventions,
    cases,
    cases_before,
    rank
)

julia

_offer_continuous_actions!(state, interventions, cases, cases_before, rank)

Once cases are settled, offer each intervention declaring settles_generation_batch its actions one case at a time, with the clock at that case's own infection time, so a wrapper such as CapacityConstrained bills capacity to the triggering case's own period rather than the generation's. Each call is scoped to the single case: an intervention such as GroupVaccination finds the rest of a triggered group through its own population-wide index, not through the candidates it is passed, so a single-case candidate still finds everything a whole-batch call would.

EpiBranch._settle_case! Function

julia

_settle_case!(
    state::SimulationState,
    case,
    interventions,
    count
)

julia

_settle_case!(state, case, interventions, count)

Run the hooks that settle one case on the continuous-time race, once its infection time is fixed and before it transmits: its clinical transitions (resolve_transitions!), then each intervention's resolve_individual!, with the state at the case as _at_case! sets it.

EpiBranch._at_case! Function

julia

_at_case!(state::SimulationState, case, count)

julia

_at_case!(state, case, count)

Set the state a case's settling hooks read: the clock (state.max_infection_time) at the case's infection time, state.cumulative_cases to count and state.current_generation to the case's generation.

EpiBranch._intervene! Function

julia

_intervene!(
    state::SimulationState,
    interventions::Vector{<:AbstractIntervention},
    targets::Vector{<:Individual},
    edges::AbstractVector{<:AbstractVector{<:Tuple}},
    minted,
    is_new
) -> Vector{Int64}

Phase 2 — contact-level interventions act on this generation's exposures. Newly minted contacts get their intervention state initialised; each target is given a provisional parent (its earliest exposing edge); then each individual of the active set is offered the contacts it reached, through the per-case hooks (_act_per_case!), before infection is resolved. An intervention still implementing the deprecated batch hook apply_post_transmission! is handed the whole batch instead, at its place in the list, and then asked keep_active about the whole batch. Returns the ids the interventions asked to keep active.

EpiBranch._act_per_case! Function

julia

_act_per_case!(
    state::SimulationState,
    interventions,
    targets,
    edges,
    is_new
) -> Vector{Int64}

julia

_act_per_case!(state, interventions, targets, edges, is_new) -> Vector{Int}

Offer each individual of the active set the contacts it reached this generation, through the per-case hooks the continuous-time race calls when it settles a case: every one of interventions' trace_contacts!, then their intervention_actions, admitted through their wrappers, then their keep_active, whose answers are returned. The caller leaves out interventions on the deprecated batch hook. A target reached by several individuals is offered to each of them, as a neighbour is on the race. Each infected case is also offered to interventions declaring traces_source (see trace_contacts!'s docstring), whether or not it reached any contact of its own this generation — unlike the forward hooks, which have nothing to act on when it did not.

Cases are visited in order of infection time, each with the state its settling hooks saw (_at_case!): the clock at its infection time and its place in the case count. An uninfected individual kept active by keep_active sees the state as the generation left it. The state is restored afterwards.

EpiBranch._case_hooks! Function

julia

_case_hooks!(
    state::SimulationState,
    case,
    interventions,
    contacts,
    is_new,
    kept
)

julia

_case_hooks!(state, case, interventions, contacts, is_new, kept)

Offer case's contacts to interventions' trace_contacts!, then their intervention_actions, admitted through their wrappers, then their keep_active, appending its answers to kept. This is the per-case hook _act_per_case! calls for each source, and the one a continuous-time race calls when it settles a case.

EpiBranch._resolve! Function

julia

_resolve!(
    model::TransmissionModel,
    state::SimulationState,
    interventions::Vector{<:AbstractIntervention},
    targets::Vector{<:Individual},
    edges::AbstractVector{<:AbstractVector{<:Tuple}},
    is_new;
    ...
)
_resolve!(
    model::TransmissionModel,
    state::SimulationState,
    interventions::Vector{<:AbstractIntervention},
    targets::Vector{<:Individual},
    edges::AbstractVector{<:AbstractVector{<:Tuple}},
    is_new,
    reinfections;
    ...
)
_resolve!(
    model::TransmissionModel,
    state::SimulationState,
    interventions::Vector{<:AbstractIntervention},
    targets::Vector{<:Individual},
    edges::AbstractVector{<:AbstractVector{<:Tuple}},
    is_new,
    reinfections,
    kept;
    intervention_tuple
)

Step 3 of a generation: decide which exposed people are infected, and record the new cases. An exposure need not lead to infection. The infector's infectiousness, the contact's susceptibility, anything the model adds and any interventions can each prevent it, and infection follows only if none of them does. A person exposed by several cases is infected if any of the exposures transmits, at the time of the earliest one that does. When someone who was infected before is reinfected (possible only once their immunity has waned, see HostImmunity), the earlier infection is kept as a past episode; see close_episode!.

EpiBranch._decide_infected Function

julia

_decide_infected(
    state::SimulationState,
    contact::Individual,
    model_risks,
    interventions,
    infected_so_far::Int64
) -> Bool

Decide whether contact is infected by the case that exposed it, and return true if so. An index case is always infected. Otherwise, in a finite population the contact first escapes with probability one minus the susceptible_fraction. Then each thing that can prevent infection is checked in turn: the contact's susceptibility and the infector's infectiousness, anything from the model's transmission_risks, and the interventions. The first one that prevents it decides.

EpiBranch.population_size Function

julia

population_size(_::TransmissionModel) -> NoPopulation

The number of people who can be infected in model, or NoPopulation() if unlimited (the default).

EpiBranch._create_individual Function

julia

_create_individual(
    state::SimulationState,
    parent_id::Int64,
    chain_id::Int64,
    next_id::Int64,
    inf_time::Real
) -> Individual

Create a new person and draw their population characteristics. They start as not infected; only the infection decision later in the generation marks them infected. Interventions set up their information on the person later in the generation, so creating a contact (in make_contact! or a model's contacts_of) involves no intervention.

EpiBranch.logsumexp Function

julia

logsumexp(x) -> Any

julia

logsumexp(x)

Numerically stable log-sum-exp. Any iterable is accepted (generators are collected first to allow two-pass computation).

EpiBranch.required_fields Function

julia

required_fields(_::AbstractTrigger) -> Vector{Symbol}

julia

required_fields(trigger) -> Vector{Symbol}

The population characteristics a trigger reads from Individual.state, such as :onset_time set by clinical_presentation. Default: none. Scheduled adds these to the fields its wrapped intervention already requires, so ModelSpec catches a model missing them the same way it catches an intervention missing a field of its own.

julia

required_fields(
    _::AbstractClinicalTransition
) -> Vector{Symbol}

julia

required_fields(transition::AbstractClinicalTransition) -> Vector{Symbol}

Case information a step needs before it can run, for example :onset_time. Each must be set by the model's attributes or recorded by another step in its progression; if neither provides one, the simulation stops with an error when it starts. None by default.

julia

required_fields(_::AbstractIntervention) -> Vector{Symbol}

The information an intervention needs recorded on each person, such as :onset_time for isolation, checked when the outbreak starts. Default: none.

EpiBranch._validate_required_fields Function

julia

_validate_required_fields(individual, items)

Check that a person has all the information the interventions or natural-history steps need (their required_fields), and raise an error naming the missing one and how to provide it.

EpiBranch._column_order Function

julia

_column_order(ks) -> Vector{Symbol}

Sensible column ordering for the linelist DataFrame: core simulation columns first, then date columns (alphabetical), then the rest (alphabetical).

EpiBranch._chain_length_ll_negbin Function

julia

_chain_length_ll_negbin(
    data,
    offspring::Distributions.NegativeBinomial;
    n_initial
) -> Any

Analytical chain length likelihood for NegBin offspring.

EpiBranch._borel_logpdf Function

julia

_borel_logpdf(μ, x::Integer) -> Any
_borel_logpdf(μ, x::Integer, s::Integer) -> Any

Log-PDF of the Borel distribution. Accepts any numeric type for μ (AD-compatible). With s > 1, this is the Borel-Tanner generalisation for the total chain size starting from s independent index cases: P(X = x | s, μ) = (s/x) * (xμ)^(x-s) * exp(-xμ) / (x-s)!.

The single-seed case (s = 1) is a port of the Borel chain-size log-likelihood (dborel / .pois_size_ll) in epichains (Azam, Funk & Finger, https://github.com/epiverse-trace/epichains, MIT).

EpiBranch._gammaborel_logpdf Function

julia

_gammaborel_logpdf(k, R, x::Integer) -> Any
_gammaborel_logpdf(k, R, x::Integer, s::Integer) -> Any

Log-PDF of the GammaBorel distribution. Accepts any numeric type for k, R (AD-compatible). With s > 1, this is the multi-seed generalisation: the PGF factors as T(z)^s and Lagrange inversion gives P(X = x | s) = (s/x) * C(kx + x - s - 1, x - s) * k^(kx) * (k+R)^(s - kx - x) * R^(x - s).

The single-seed case (s = 1) is a port of .nbinom_size_ll in epichains (Azam, Funk & Finger, https://github.com/epiverse-trace/epichains, MIT).

EpiBranch._chain_size_loglik Function

julia

_chain_size_loglik(
    dist,
    data::ChainSizes;
    prob_concluded
) -> Any

julia

_chain_size_loglik(dist, data::ChainSizes; prob_concluded = nothing)

Sum over clusters of the chain-size log-likelihood under dist: each cluster's log P(X = x_i | seeds_i), or with prob_concluded the mixture π_i · P(X = x_i) + (1 − π_i) · P(X ≥ x_i) for clusters that may still grow.

EpiBranch._empirical_ll Function

julia

_empirical_ll(
    observed,
    simulated;
    min_val,
    censored,
    cap,
    prob_concluded
) -> Any

Empirical log-likelihood with Laplace smoothing and right-censoring.

When censored is provided, simulated values flagged as censored contribute to P(size >= cap) rather than P(size = cap). Observed values at or above cap are evaluated as P(size >= cap).

prob_concluded, when given (one value per observation, in [0, 1]), mixes each observation's point mass with its right tail under the same simulated empirical law, exactly as the analytical _chain_size_loglik mixture does: π_i · P(X = x_i) + (1 − π_i) · P(X ≥ x_i).

API reference ​

Transmission models ​

Population characteristics ​

Clinical progression (natural history) ​

Interventions ​

Isolation ​

Contact tracing ​

Vaccination ​

Timing and limited resources ​

Transmission routes ​

Running simulations ​

Recording contacts ​

Information about each case ​

Output ​

Closed-form results ​

Reproduction number, extinction and superspreading ​

Chain-size distributions ​

Other chain-size distributions ​

A different offspring distribution for the index case ​

Household reproduction number ​

Inference ​

Data types ​

Likelihood and fitting ​

Distributions for Bayesian fitting ​

Incidence time-series likelihood ​

Observation models ​

Reproduction number varying between chains ​

Pairwise survival likelihood ​

Household data ​

Network data ​

Distribution helpers ​

Default placeholders ​

For extension authors ​

Writing a transmission model ​

Stopping rules ​

Writing an intervention ​

Isolation and contact tracing rules ​

Vaccination ​

Interventions with limited resources ​

Start triggers ​

Writing a transition ​

Population characteristics and recorded contacts ​

Writing an observation model ​

Pairwise likelihood building blocks ​

Internal functions ​

API reference

Transmission models

Population characteristics

Clinical progression (natural history)

Interventions

Isolation

Contact tracing

Vaccination

Timing and limited resources

Transmission routes

Running simulations

Recording contacts

Information about each case

Output

Closed-form results

Reproduction number, extinction and superspreading

Chain-size distributions

Other chain-size distributions

A different offspring distribution for the index case

Household reproduction number

Inference

Data types

Likelihood and fitting

Distributions for Bayesian fitting

Incidence time-series likelihood

Observation models

Reproduction number varying between chains

Pairwise survival likelihood

Household data

Network data

Distribution helpers

Default placeholders

For extension authors

Writing a transmission model

Stopping rules

Writing an intervention

Isolation and contact tracing rules

Vaccination

Interventions with limited resources

Start triggers

Writing a transition

Population characteristics and recorded contacts

Writing an observation model

Pairwise likelihood building blocks

Internal functions