Public Documentation

Module for defining epidemiological models.

struct DirectInfections{S<:Distributions.Sampleable} <: AbstractTuringEpiModel

Model unobserved/latent infections as a transformation on a sampled latent process.

Mathematical specification

If $Z_t$ is a realisation of the latent model, then the unobserved/latent infections are given by

\[I_t = g(\hat{I}_0 + Z_t).\]

where $g$ is a transformation function and the unconstrained initial infections $\hat{I}_0$ are sampled from a prior distribution.

DirectInfections are constructed by passing an EpiData object data and an initialisation_prior for the prior distribution of $\hat{I}_0$. The default initialisation_prior is Normal().

Constructors

• DirectInfections(; data, initialisation_prior)

Example usage with generate_latent_infs

generate_latent_infs can be used to construct a Turing model for the latent infections conditional on the sample path of a latent process. In this example, we generate a sample of a white noise latent process.

First, we construct a DirectInfections struct with an EpiData object, an initialisation prior and a transformation function.

using Distributions, Turing, EpiAware
gen_int = [0.2, 0.3, 0.5]
g = exp

# Create an EpiData object
data = EpiData(gen_int, g)

# Create a DirectInfections model
direct_inf_model = DirectInfections(data = data, initialisation_prior = Normal())

Then, we can use generate_latent_infs to construct a Turing model for the unobserved infection generation model set by the type of direct_inf_model.

# Construct a Turing model
Z_t = randn(100)
latent_inf = generate_latent_infs(direct_inf_model, Z_t)

Now we can use the Turing PPL API to sample underlying parameters and generate the unobserved infections.

# Sample from the unobserved infections model

#Sample random parameters from prior
θ = rand(latent_inf)
#Get unobserved infections as a generated quantities from the model
I_t = generated_quantities(latent_inf, θ)

Fields

• data::EpiData: Epidata object.

• initialisation_prior::Distributions.Sampleable: Prior distribution for the initialisation of the infections. Default is Normal().

struct EpiData{T<:Real, F<:Function}

The EpiData struct represents epidemiological data used in infectious disease modeling.

Constructors

• EpiData(gen_int, transformation::Function). Constructs an EpiData object with discrete

generation interval gen_int and transformation function transformation.

• EpiData(;gen_distribution::ContinuousDistribution, D_gen, Δd = 1.0, transformation::Function = exp).

Constructs an EpiData object with double interval censoring discretisation of the continuous next generation interval distribution gen_distribution with additional right truncation at D_gen. Δd sets the interval width (default = 1.0). transformation sets the transformation function

Examples

Construction direct from discrete generation interval and transformation function:

using EpiAware
gen_int = [0.2, 0.3, 0.5]
g = exp
data = EpiData(gen_int, g)

Construction from continuous distribution for generation interval.

using Distributions

gen_distribution = Uniform(0.0, 10.0)

data = EpiData(;gen_distribution
    D_gen = 10.0)

Fields

• gen_int::Vector{T} where T<:Real: Discrete generation interval.

• len_gen_int::Integer: Length of the discrete generation interval.

• transformation::Function: Transformation function defining constrained and unconstrained domain bijections.

struct ExpGrowthRate{S<:Distributions.Sampleable} <: AbstractTuringEpiModel

Model unobserved/latent infections as due to time-varying exponential growth rate $r_t$ which is generated by a latent process.

Mathematical specification

If $Z_t$ is a realisation of the latent model, then the unobserved/latent infections are given by

\[I_t = g(\hat{I}_0) \exp(Z_t).\]

where $g$ is a transformation function and the unconstrained initial infections $\hat{I}_0$ are sampled from a prior distribution.

ExpGrowthRate are constructed by passing an EpiData object data and an initialisation_prior for the prior distribution of $\hat{I}_0$. The default initialisation_prior is Normal().

Constructor

• ExpGrowthRate(; data, initialisation_prior).

Example usage with generate_latent_infs

generate_latent_infs can be used to construct a Turing model for the latent infections conditional on the sample path of a latent process. In this example, we generate a sample of a white noise latent process.

First, we construct an ExpGrowthRate struct with an EpiData object, an initialisation prior and a transformation function.

using Distributions, Turing, EpiAware
gen_int = [0.2, 0.3, 0.5]
g = exp

# Create an EpiData object
data = EpiData(gen_int, g)

# Create an ExpGrowthRate model
exp_growth_model = ExpGrowthRate(data = data, initialisation_prior = Normal())

Then, we can use generate_latent_infs to construct a Turing model for the unobserved infection generation model set by the type of direct_inf_model.

# Construct a Turing model
Z_t = randn(100) * 0.05
latent_inf = generate_latent_infs(exp_growth_model, Z_t)

Now we can use the Turing PPL API to sample underlying parameters and generate the unobserved infections.

# Sample from the unobserved infections model

#Sample random parameters from prior
θ = rand(latent_inf)
#Get unobserved infections as a generated quantities from the model
I_t = generated_quantities(latent_inf, θ)

Fields

• data::EpiData

• initialisation_prior::Distributions.Sampleable

struct ODEProcess{P<:AbstractTuringLatentModel, S, F<:Function, D<:Union{Dict, NamedTuple}} <: AbstractTuringEpiModel

A structure representing an infection process modeled by an Ordinary Differential Equation (ODE). At a high level, an ODEProcess struct object combines:

• An AbstractTuringParamModel which defines the ODE model in terms of OrdinaryDiffEq types, the parameters of the ODE model and a method to generate the parameters.
• A technique for solving and interpreting the ODE model using the SciML ecosystem. This includes the solver used in the ODE solution, keyword arguments to send to the solver and a function to map the ODESolution solution object to latent infections.

Constructors

• ODEProcess(prob::ODEProblem; ts, solver, sol2infs): Create an ODEProcess

object with the ODE problem prob, time points ts, solver solver, and function sol2infs.

Predefined ODE models

Two basic ODE models are provided in the EpiAware package: SIRParams and SEIRParams. In both cases these are defined in terms of the proportions of the population in each compartment of the SIR and SEIR models respectively.

SIR model

\[\begin{aligned} \frac{dS}{dt} &= -\beta SI \\ \frac{dI}{dt} &= \beta SI - \gamma I \\ \frac{dR}{dt} &= \gamma I \end{aligned}\]

Where S is the proportion of the population that is susceptible, I is the proportion of the population that is infected and R is the proportion of the population that is recovered. The parameters are the infectiousness β and the recovery rate γ.

using EpiAware, OrdinaryDiffEq, Distributions

# Create an instance of SIRParams
sirparams = SIRParams(
    tspan = (0.0, 100.0),
    infectiousness = LogNormal(log(0.3), 0.05),
    recovery_rate = LogNormal(log(0.1), 0.05),
    initial_prop_infected = Beta(1, 99)
)
nothing

SEIR model

\[\begin{aligned} \frac{dS}{dt} &= -\beta SI \\ \frac{dE}{dt} &= \beta SI - \alpha E \\ \frac{dI}{dt} &= \alpha E - \gamma I \\ \frac{dR}{dt} &= \gamma I \end{aligned}\]

Where S is the proportion of the population that is susceptible, E is the proportion of the population that is exposed, I is the proportion of the population that is infected and R is the proportion of the population that is recovered. The parameters are the infectiousness β, the incubation rate α and the recovery rate γ.

using EpiAware, OrdinaryDiffEq, Distributions, Random
Random.seed!(1234)

# Create an instance of SIRParams
seirparams = SEIRParams(
    tspan = (0.0, 100.0),
    infectiousness = LogNormal(log(0.3), 0.05),
    incubation_rate = LogNormal(log(0.1), 0.05),
    recovery_rate = LogNormal(log(0.1), 0.05),
    initial_prop_infected = Beta(1, 99)
)
nothing

Usage example with ODEProcess and predefined SIR model

In this example we define an ODEProcess object using the predefined SIRParams model from above. We then generate latent infections using the generate_latent_infs function, and refit the model using a Turing model.

We assume that the latent infections are observed with a Poisson likelihood around their ODE model prediction. The population size is N = 1000, which we put into the sol2infs function, which maps the ODE solution to the number of infections. Recall that the EpiAware default SIR implementation assumes the model is in density/proportion form. Also, note that since the sol2infs function is a link function that maps the ODE solution to the expected number of infections we also apply the LogExpFunctions.softplus function to ensure that the expected number of infections is non-negative. Note that the softplus function is a smooth approximation to the ReLU function x -> max(0, x). The utility of this approach is that small negative output from the ODE solver (e.g. ~ -1e-10) will be mapped to small positive values, without needing to use strict positivity constraints in the model.

First, we define the ODEProcess object which combines the SIR model with the sol2infs link function and the solver options.

using Turing, LogExpFunctions
N = 1000.0

sir_process = ODEProcess(
    params = sirparams,
    sol2infs = sol -> softplus.(N .* sol[2, :]),
    solver_options = Dict(:verbose => false, :saveat => 1.0)
)
nothing

Second, we define a PoissionError observation model for linking the the number of infections.

pois_obs = PoissonError()
nothing

Next, we create a Turing model for the full generative process: this solves the ODE model for the latent infections and then samples the observed infections from a Poisson distribution with this as the average.

NB: The nothing argument is a dummy latent process, e.g. a log-Rt time series, that is not used in the SIR model, but might be used in other models.

@model function fit_ode_model(data)
    @submodel I_t = generate_latent_infs(sir_process, nothing)
    @submodel y_t = generate_observations(pois_obs, data, I_t)

    return y_t
end
nothing

We can generate some test data from the model by passing missing as the argument to the model. This tells Turing that there is no data to condition on, so it will sample from the prior parameters and then generate infections. In this case, we do it in a way where we cache the sampled parameters as θ for later use.

# Sampled parameters
gen_mdl = fit_ode_model(missing)
θ = rand(gen_mdl)
test_data = (gen_mdl | θ)()
nothing

Now, we can refit the model but this time we condition on the test data. We suppress the output of the sampling process to keep the output clean, but you can remove the @suppress macro.

using Suppressor
inference_mdl = fit_ode_model(test_data)
chn = Suppressor.@suppress sample(inference_mdl, NUTS(), 2_000)
summarize(chn)
nothing

We can compare the summarized chain to the sampled parameters in θ to see that the model is fitting the data well and recovering a credible interval containing the true parameters.

Custom ODE models

To define a custom ODE model, you need to define:

• Some CustomModel <: AbstractTuringLatentModel struct that contains the ODE problem as a field called prob, as well as sufficient fields to define or sample the parameters of the ODE model.
• A method for EpiAwareBase.generate_latent(params::CustomModel, Z_t) that generates the initial condition and parameters of the ODE model, potentially conditional on a sample from a latent process Z_t. This method must return a Tuple (u0, p) where u0 is the initial condition and p is the parameters.

Here is an example of a simple custom ODE model for specified exponential growth:

using EpiAware, Turing, OrdinaryDiffEq
# Define a simple exponential growth model for testing
function expgrowth(du, u, p, t)
    du[1] = p[1] * u[1]
end

r = log(2) / 7 # Growth rate corresponding to 7 day doubling time

# Define the ODE problem using SciML
prob = ODEProblem(expgrowth, [1.0], (0.0, 10.0), [r])

# Define the custom parameters struct
struct CustomModel <: AbstractTuringLatentModel
    prob::ODEProblem
    r::Float64
    u0::Float64
end
custom_ode = CustomModel(prob, r, 1.0)

# Define the custom generate_latent function
@model function EpiAwareBase.generate_latent(params::CustomModel, n)
    return ([params.u0], [params.r])
end
nothing

This model is not random! But we can still use it to generate latent infections.

# Define the ODEProcess
expgrowth_model = ODEProcess(
    params = custom_ode,
    sol2infs = sol -> sol[1, :]
)
infs = generate_latent_infs(expgrowth_model, nothing)()
nothing

Fields

• params::AbstractTuringLatentModel: The ODE problem and parameters, where P is a subtype of AbstractTuringLatentModel.

• solver::Any: The solver used for the ODE problem. Default is AutoVern7(Rodas5()), which is an auto switching solver aimed at medium/low tolerances.

• sol2infs::Function: A function that maps the solution object of the ODE to infection counts.

• solver_options::Union{Dict, NamedTuple}: The extra solver options for the ODE problem. Can be either a Dict or a NamedTuple containing the solver options.

struct Renewal{E, S<:Distributions.Sampleable, A} <: AbstractTuringRenewal

Model unobserved/latent infections as due to time-varying Renewal model with reproduction number $\mathcal{R}_t$ which is generated by a latent process.

Mathematical specification

If $Z_t$ is a realisation of the latent model, then the unobserved/latent infections are given by

\[\begin{align} \mathcal{R}_t &= g(Z_t),\\ I_t &= \mathcal{R}_t \sum_{i=1}^{n-1} I_{t-i} g_i, \qquad t \geq 1, \\ I_t &= g(\hat{I}_0) \exp(r(\mathcal{R}_1) t), \qquad t \leq 0. \end{align}\]

where $g$ is a transformation function and the unconstrained initial infections $\hat{I}_0$ are sampled from a prior distribution. The discrete generation interval is given by $g_i$.

$r(\mathcal{R}_1)$ is the exponential growth rate implied by $\mathcal{R}_1)$ using the implicit relationship between the exponential growth rate and the reproduction number.

\[\mathcal{R} \sum_{j \geq 1} g_j \exp(- r j)= 1.\]

Renewal are constructed by passing an EpiData object data and an initialisation_prior for the prior distribution of $\hat{I}_0$. The default initialisation_prior is Normal().

Constructors

• Renewal(; data, initialisation_prior). Construct a Renewal model with default update steps.
• Renewal(data; initialisation_prior). Construct a Renewal model with default update steps.
• Renewal(data, initialisation_prior, recurrent_step) Construct a Renewal model with recurrent_step update step function.

Example usage with generate_latent_infs

generate_latent_infs can be used to construct a Turing model for the latent infections conditional on the sample path of a latent process. In this example, we generate a sample of a white noise latent process.

First, we construct an Renewal struct with an EpiData object, an initialisation prior and a transformation function.

using Distributions, Turing, EpiAware
gen_int = [0.2, 0.3, 0.5]
g = exp

# Create an EpiData object
data = EpiData(gen_int, g)

# Create an Renewal model
renewal_model = Renewal(data; initialisation_prior = Normal())

Then, we can use generate_latent_infs to construct a Turing model for the unobserved infection generation model set by the type of direct_inf_model.

# Construct a Turing model
Z_t = randn(100) * 0.05
latent_inf = generate_latent_infs(renewal_model, Z_t)

Now we can use the Turing PPL API to sample underlying parameters and generate the unobserved infections.

# Sample from the unobserved infections model

#Sample random parameters from prior
θ = rand(latent_inf)
#Get unobserved infections as a generated quantities from the model
I_t = generated_quantities(latent_inf, θ)

Fields

• data::Any

• initialisation_prior::Distributions.Sampleable

• recurrent_step::Any

R_to_r(
    R₀,
    w::Array{T<:AbstractFloat, 1};
    newton_steps,
    Δd
) -> Any

This function computes an approximation to the exponential growth rate r given the reproductive ratio R₀ and the discretized generation interval w with discretized interval width Δd. This is based on the implicit solution of

\[G(r) - {1 \over R_0} = 0.\]

where

\[G(r) = \sum_{i=1}^n w_i e^{-r i}.\]

is the negative moment generating function (MGF) of the generation interval distribution.

The two step approximation is based on: 1. Direct solution of implicit equation for a small r approximation. 2. Improving the approximation using Newton's method for a fixed number of steps newton_steps.

Returns:

• The approximate value of r.

expected_Rt(
    data::EpiData,
    infections::Vector{<:Real}
) -> Any

Calculate the expected Rt values based on the given EpiData object and infections.

\[R_t = \frac{I_t}{\sum_{i=1}^{n} I_{t-i} g_i}\]

Arguments

• data::EpiData: An instance of the EpiData type containing generation interval data.
• infections::Vector{<:Real}: A vector of infection data.

Returns

• exp_Rt::Vector{Float64}: A vector of expected Rt values.

Examples

using EpiAware

data = EpiData([0.2, 0.3, 0.5], exp)
infections = [100, 200, 300, 400, 500]
expected_Rt(data, infections)

r_to_R(r, w::AbstractVector) -> Any

r_to_R(r, w)

Compute the reproductive ratio given exponential growth rate r and discretized generation interval w.

Arguments

• r: The exponential growth rate.
• w: discretized generation interval.

Returns

• The reproductive ratio.