Fits the parameters of a BiSSeCo model by maximum likelihood

This method estimates a pair of relative fitness coefficients from pathogen phylogenies representing differences in transmissibility and differences in within-host fitness. Relative fitness is paramaterised by 1) the parameter "alpha" which represents the fold-change in within-host replicative fitness and 2) the parameter "omega" which represents the fold-change in between-host replicative fitness. See details for more information.

Usage

fitbisseco(
  tr,
  isvariant,
  Tg,
  mu,
  Net = NULL,
  theta0 = log(c(alpha = 15, omega = 0.95, yscale = 1)),
  augment_likelihood = TRUE,
  optim_parms = list(),
  mlesky_parms = list(),
  res = 200,
  ...
)

Arguments

tr: an ape::phylo representing a time-scaled phylogeny. These can be computed with the treedater package
isvariant: vector of type boolean with length equal to ape::Ntip(tr). Each element is TRUE if the corresponding element in tr$tip.label is a variant type
Tg: Generation time, i.e. the mean time elapsed between generations. Units of this variable should match those used in tr$edge.length and 1/mu, e.g. days or years
mu: Molecular clock rate of evolution
Net: NULL or a matrix with two columns giving the effective population size. If NULL, the effective population size will be computed with the mlesky package. The first column should be time since some point in the past, and the second column should be an estimate of the effective population size. Time zero should correspond to the root of the tree
theta0: Initial guess of (log) parameter values: alpha, omega, and Ne scale. Can be a named vector or in the aforementioned order.
augment_likelihood: if TRUE (default), will combine the coalescent likelihood with a binomial likelihood of sampling variants or ancestral types under the assumption of random sampling and mutation-selection balance
optim_parms: optional list of parameters passed to optim when fitting the model
mlesky_parms: optional list of parameters passed to mlskygrid if estimating Ne(t)
res: Integer number time steps used in coalescent likelihood
...: Additional parameters are passed to colik

Value

A fitted BiSSeCo model with coef and summary methods

Details

Pathogen phylogenies are assumed to be reconstructed from population-based random samples of pathogen genomes and at most one sequence per host. Phylogenies should be time-scaled, and an estimate of the molecular clock rate of evolution should be provided, such as estimated with the treedater package. If not provided, the effective population size over time is estimated with the mlesky package, and this estimate is also provided in the returned BiSSeCo fit.

Examples


# load tree simulated with TiPS
tr <- readRDS(system.file("extdata/tips_params_1_rep_2.rds", package = "musseco"))

# set to TRUE tips that are from variant
isvariant <- grepl( tr$tip.label, pattern = 'V_' )
isvariant <- setNames(isvariant, tr$tip.label)

#this will take a few minutes to run
#the values of mu and Tg was based on true values used to simulate the trees
#see vignette for more details.
if (FALSE) { # \dontrun{
fb_au <- fitbisseco( tr,
                     isvariant,
                     Tg = 10.2,
                     mu = 0.0016,
                     Net = NULL,
                     theta0 = log(c(2,.75,1/2)),
                     augment_likelihood = TRUE,
                     optim_parms = list(),
                     mlesky_parms = list(tau = NULL,
                                         tau_lower = 0.1,
                                         tau_upper = 1e7,
                                         ncpu = 1, 
                                         model = 1 ) )
} # }
# you can load the results from here
load(system.file("extdata/fit.rda", package = "musseco"))
fb_au
#> Binary state speciation and extinction coalescent fit
#> 
#>               2.5% 97.5%
#> alpha  0.910 0.570 1.453
#> omega  0.918 0.874 0.964
#> yscale 0.844 0.769 0.926
#> 
#> Likelihood: -5331.4195