congruent_hbds_model: Generate a congruent homogenous-birth-death-sampling model.

Description

Given a homogenous birth-death-sampling (HBDS) model (or abstract congruence class), obtain the congruent model (or member of the congruence class) with a specific speciation rate \(\lambda\), or extinction rate \(\mu\), or sampling rate \(\psi\), or effective reproduction ratio \(R_e\) or removal rate \(\mu+\psi\) (aka. "become uninfectious"" rate). All input and output time-profiles are specified as piecewise polynomial curves (splines), defined on a discrete grid of ages. This function allows exploration of a model's congruence class, by obtaining various members of the congruence class depending on the specified \(\lambda\), \(\mu\), \(\psi\), \(R_e\) or removal rate. For more details on HBDS models and congruence classes see the documentation of simulate_deterministic_hbds.

Usage

congruent_hbds_model(age_grid,
                     PSR,
                     PDR,
                     lambda_psi,
                     lambda             = NULL,
                     mu                 = NULL,
                     psi                = NULL,
                     Reff               = NULL,
                     removal_rate       = NULL,
                     lambda0            = NULL,
                     CSA_ages           = NULL,
                     CSA_pulled_probs   = NULL,
                     CSA_PSRs           = NULL,
                     splines_degree     = 1,
                     ODE_relative_dt    = 0.001,
                     ODE_relative_dy    = 1e-4)

Value

A named list with the following elements:

success: Logical, indicating whether the calculation was successful. If FALSE, then the returned list includes an additional `error' element (character) providing a description of the error; all other return variables may be undefined.
valid: Logical, indicating whether the returned model appears to be biologically valid (for example, does not have negative \(\lambda\), \(\mu\) or \(\psi\)). In principle, a congruence class may include biologically invalid models, which might be returned depending on the input to congruent_hbds_model. Note that only biologically valid models can be subsequently simulated using simulate_deterministic_hbds.
ages: Numeric vector of size NG, specifying the discrete ages (time before present) on which all returned time-curves are specified. Will always be equal to age_grid.
lambda: Numeric vector of size NG, listing the speciation rates \(\lambda\) of the returned model at the ages given in ages[].
mu: Numeric vector of size NG, listing the extinction rates \(\mu\) of the returned model at the ages given in ages[].
psi: Numeric vector of size NG, listing the (Poissonian) sampling rates \(\psi\) of the returned model at the ages given in ages[].
lambda_psi: Numeric vector of size NG, listing the product \(\lambda\psi\) at the ages given in ages[].
Reff: Numeric vector of size NG, listing the effective reproduction ratio \(R_e\) of the returned model at the ages given in ages[].
removal_rate: Numeric vector of size NG, listing the removal rate (\(\mu+\psi\), aka. "become uninfectious" rate) of the returned model at the ages given in ages[].
Pmissing: Numeric vector of size NG, listing the probability that a past lineage extant during ages[] will be missing from a tree generated by the model.
CSA_probs: Numeric vector of the same size as CSA_ages, listing the sampling probabilities at each of the CSAs.
CSA_Pmissings: Numeric vector of the same size as CSA_ages, listing the probability that a past lineage extant during each of CSA_ages[] will be missing from a tree generated by the model.

Arguments

age_grid: Numeric vector, listing discrete ages (time before present) on which the various model variables (e.g., \(\lambda\), \(\mu\) etc) are specified. Listed ages must be strictly increasing, and must include the present-day (i.e. age 0).
PSR: Numeric vector, of the same size as age_grid, specifying the pulled speciation rate (PSR) (in units 1/time) at the ages listed in age_grid. The PSR is assumed to vary polynomially between grid points (see argument splines_degree). Can also be a single number, in which case PSR is assumed to be time-independent.
PDR: Numeric vector, of the same size as age_grid, specifying the pulled diversification rate (PDR) (in units 1/time) at the ages listed in age_grid. The PDR is assumed to vary polynomially between grid points (see argument splines_degree). The PDR of a HBDS model is defined as \(PDR=\lambda-\mu-\psi+(1/\lambda)d\lambda/dt\) (where \(t\) denotes age). Can also be a single number, in which case PDR is assumed to be time-independent.
lambda_psi: Numeric vector, of the same size as age_grid, specifying the product of speciation rate and sampling rate (\(\lambda\psi\), in units 1/time^2) at the ages listed in age_grid. \(\lambda\psi\) is assumed to vary polynomially between grid points (see argument splines_degree). Can also be a single number, in which case \(\lambda\psi\) is assumed to be time-independent.
lambda: Numeric vector, of the same size as age_grid, specifying the speciation rate (\(\lambda\), in units 1/time) at the ages listed in age_grid. The speciation rate is assumed to vary polynomially between grid points (see argument splines_degree). Can also be a single number, in which case \(\lambda\) is assumed to be time-independent. By providing \(\lambda\), one can select a specific model from the congruence class. Note that exactly one of lambda, mu, psi, Reff or removal_rate must be provided.
mu: Numeric vector, of the same size as age_grid, specifying the extinction rate (\(\mu\), in units 1/time) at the ages listed in age_grid. The extinction rate is assumed to vary polynomially between grid points (see argument splines_degree). Can also be a single number, in which case \(\mu\) is assumed to be time-independent. In an epidemiological context, \(\mu\) typically corresponds to the recovery rate plus the death rate of infected individuals. By providing \(\mu\) (together with lambda0, see below), one can select a specific model from the congruence class. Note that exactly one of lambda, mu, psi, Reff or removal_rate must be provided.
psi: Numeric vector, of the same size as age_grid, specifying the (Poissonian) sampling rate (\(\psi\), in units 1/time) at the ages listed in age_grid. The sampling rate is assumed to vary polynomially between grid points (see argument splines_degree). Can also be a single number, in which case \(\psi\) is assumed to be time-independent. By providing \(\psi\), one can select a specific model from the congruence class. Note that exactly one of lambda, mu, psi, Reff or removal_rate must be provided.
Reff: Numeric vector, of the same size as age_grid, specifying the effective reproduction ratio (\(R_e\), unitless) at the ages listed in age_grid. The \(R_e\) is assumed to vary polynomially between grid points (see argument splines_degree). Can also be a single number, in which case \(R_e\) is assumed to be time-independent. By providing \(R_e\) (together with lambda0, see below), one can select a specific model from the congruence class. Note that exactly one of lambda, mu, psi, Reff or removal_rate must be provided.
removal_rate: Numeric vector, of the same size as age_grid, specifying the removal rate (\(\mu+\psi\), in units 1/time) at the ages listed in age_grid. IN an epidemiological context this is also known as "become uninfectious" rate. The removal rate is assumed to vary polynomially between grid points (see argument splines_degree). Can also be a single number, in which case the removal rate is assumed to be time-independent. By providing \(\mu+\psi\) (together with lambda0, see below), one can select a specific model from the congruence class. Note that exactly one of lambda, mu, psi, Reff or removal_rate must be provided.
lambda0: Numeric, specifying the speciation rate at the present-day (i.e., at age 0). Must be provided if and only if one of mu, Reff or removal_rate is provided.
CSA_ages: Optional numeric vector, listing the ages of concentrated sampling attempts, in ascending order. Concentrated sampling is assumed to occur in addition to any continuous (Poissonian) sampling specified by psi.
CSA_pulled_probs: Optional numeric vector of the same size as CSA_ages, listing pulled sampling probabilities at each concentrated sampling attempt (CSA). Note that in contrast to the sampling rates psi, the CSA_pulled_probs are interpreted as probabilities and must thus be between 0 and 1. CSA_pulled_probs must be provided if and only if CSA_ages is provided.
CSA_PSRs: Optional numeric vector of the same size as CSA_ages, specifying the pulled sampling rate (PSR) during each concentrated sampling attempt. While in principle the PSR is already provided by the argument PSR, the PSR may be non-continuous at CSAs, which makes a representation as piecewise polynomial function difficult; hence, you must explicitly provide the correct PSR at each CSA. CSA_PSRs must be provided if and only if CSA_ages is provided.
splines_degree: Integer, either 0,1,2 or 3, specifying the polynomial degree of the provided time-dependent variables between grid points in age_grid. For example, if splines_degree==1, then the provided PDR, PSR and so on are interpreted as piecewise-linear curves; if splines_degree==2 they are interpreted as quadratic splines; if splines_degree==3 they are interpreted as cubic splines. The splines_degree influences the analytical properties of the curve, e.g. splines_degree==1 guarantees a continuous curve, splines_degree==2 guarantees a continuous curve and continuous derivative, and so on.
ODE_relative_dt: Positive unitless number, specifying the default relative time step for internally used ordinary differential equation solvers. Typical values are 0.01-0.001.
ODE_relative_dy: Positive unitless number, specifying the relative difference between subsequent simulated and interpolated values, in internally used ODE solvers. Typical values are 1e-2 to 1e-5. A smaller ODE_relative_dy increases interpolation accuracy, but also increases memory requirements and adds runtime.

Author

Stilianos Louca

Details

The PDR, PSR and the product \(\lambda\psi\) are variables that are invariant across the entire congruence class of an HBDS model, i.e. any two congruent models have the same PSR, PDR and product \(\lambda\psi\). Reciprocally, any HBDS congruence class is fully determined by its PDR, PSR and \(\lambda\psi\). This function thus allows "collapsing" a congruence class down to a single member (a specific HBDS model) by specifying one or more additional variables over time (such as \(\lambda\), or \(\psi\), or \(\mu\) and \(\lambda_0\)). Alternatively, this function may be used to obtain alternative models that are congruent to some reference model, for example to explore the robustness of specific epidemiological quantities of interest. The function returns a specific HBDS model in terms of the time profiles of various variables (such as \(\lambda\), \(\mu\) and \(\psi\)).

In the current implementation it is assumed that any sampled lineages are immediately removed from the pool, that is, this function cannot accommodate models with a non-zero retention probability upon sampling. This is a common assumption in molecular epidemiology. Note that in this function age always refers to time before present, i.e., present day age is 0, and age increases towards the root.

References

T. Stadler, D. Kuehnert, S. Bonhoeffer, A. J. Drummond (2013). Birth-death skyline plot reveals temporal changes of epidemic spread in HIV and hepatitis C virus (HCV). PNAS. 110:228-233.

A. MacPherson, S. Louca, A. McLaughlin, J. B. Joy, M. W. Pennell (in review as of 2020). A general birth-death-sampling model for epidemiology and macroevolution. DOI:10.1101/2020.10.10.334383

Examples

Run this code

# define an HBDS model with exponentially decreasing speciation/extinction rates
# and constant Poissonian sampling rate psi
oldest_age= 10
age_grid  = seq(from=0,to=oldest_age,by=0.1) # choose a sufficiently fine age grid
lambda    = 1*exp(0.01*age_grid) # define lambda on the age grid
mu        = 0.2*lambda # assume similarly shaped but smaller mu

# simulate deterministic HBD model
# scale LTT such that it is 100 at age 1
sim = simulate_deterministic_hbds(age_grid  = age_grid,
                                  lambda    = lambda,
                                  mu        = mu,
                                  psi       = 0.1,
                                  age0      = 1,
                                  LTT0      = 100)
                                         
# calculate a congruent HBDS model with an alternative sampling rate
# use the previously simulated variables to define the congruence class
new_psi = 0.1*exp(-0.01*sim$ages) # consider a psi decreasing with age
congruent = congruent_hbds_model(age_grid   = sim$ages,
                                 PSR        = sim$PSR,
                                 PDR        = sim$PDR,
                                 lambda_psi = sim$lambda_psi,
                                 psi        = new_psi)
											 
# compare the deterministic LTT of the two models
# to confirm that the models are indeed congruent
if(!congruent$valid){
    cat("WARNING: Congruent model is not biologically valid\n")
}else{
    # simulate the congruent model to get the LTT
    csim = simulate_deterministic_hbds(age_grid = congruent$ages,
                                       lambda   = congruent$lambda,
                                       mu       = congruent$mu,
                                       psi      = congruent$psi,
                                       age0     = 1,
                                       LTT0     = 100)

    # plot deterministic LTT of the original and congruent model
    plot( x = sim$ages, y = sim$LTT, type='l',
          main='dLTT', xlab='age', ylab='lineages', 
          xlim=c(oldest_age,0), col='red')
    lines(x= csim$ages, y=csim$LTT, 
          type='p', pch=21, col='blue')
}

Run the code above in your browser using DataLab