Compute a structured coalescent likelihood given a dated genealogy and a demographic history in FGY format

Usage

colik.pik.fgy(
  tree,
  tfgy,
  timeOfOriginBoundaryCondition = TRUE,
  maxHeight = Inf,
  forgiveAgtY = 1,
  AgtY_penalty = 1,
  returnTree = FALSE,
  step_size_res = 10,
  PL2 = FALSE
)

Arguments

tree: A DatedTree object
tfgy: class tfgy that corresponds to the time, births, migrations and number of infections. See examples for how to construct one.
timeOfOriginBoundaryCondition: If TRUE, will return -Inf if the root of the tree precedes the time of origin.
maxHeight: It will only count internode intervals in the likelihood that occur after maxHeight years before present. Useful for large trees and when you do not want to model the entire demographic history.
forgiveAgtY: If number of extant lineages exceeds simulated population size, return -Inf if this value is zero, or forgive the discrepancy if zero. If between zero and one, only forgive the discrepancy if this proportion of lineages is less than the given value.
AgtY_penalty: If number of extant lineages exceeds simulated population size, penalize likelihood with value L*AgtY_penalty where L is the cumulative coalescent rate within the given internode interval. 0<= AgtY_penalty <= Inf.
returnTree: If TRUE, a copy of the tree is also returned, which includes the inferred states of lineages and likelihood terms at each internal node.
step_size_res: Parameter for the ODE solver; it is the default number of timesteps to use when solving coalescent equations in each internode interval
PL2: Toggle likelihood approximation to be used. If TRUE, the likelihood PL2 will be used. If PL2 = FALSE, the likelihood PL1 will be used. (PL1 faster better approximation, PL2 slow/good approximation). See Volz & Siveroni 2018 for details.

Value

The coalescent (or structured coalescent) log likelihood (numeric).

Author

Erik Volz

Examples

# A simple exponential growth model with birth rates beta and death rates gamma:
dm <- build.demographic.process(births = c(I = 'parms$beta * I'),
                                deaths = c(I = 'parms$gamma * I'),
                                parameterNames = c('beta', 'gamma'),
                                rcpp = FALSE,
                                sde = FALSE)

# simulate a tree based on the model:
tre <- sim.co.tree(list(beta = 1.5, gamma = 1),
                        dm,
                        x0  = c(I = 1),
                        t0 = 0,
                        sampleTimes = seq(10, 15, length.out=50),
                        res = 1000)
                        
 #Create the tfgy using the dm function above
 tfgy <- dm(theta = list(beta = 1.5, gamma = 1),
            x0 = c(I = 1),
            t0 = 0,
            t1 = 15)
 
 # Compute a likelihood
 colik.pik.fgy(tre, tfgy)
#> [1] -207.836