Simulate a coalescent genealogy based on a demographic process.

A tree is simulated based on a demographic process and user-supplied parameters and initial conditions. The times and types of samples must also be supplied.

Usage

sim.co.tree(
  theta,
  demographic.process.model,
  x0,
  t0,
  sampleTimes,
  sampleStates = NULL,
  res = 1000,
  integrationMethod = "lsoda",
  finiteSizeCorrections = TRUE,
  substitutionRates = NULL,
  sequenceLength = 0
)

Arguments

theta

A named vector of parameter values required by the model

demographic.process.model

An object of class demographic.process. See build.demographic.process

x0

A named vector of initial conditions required by the model. This includes demes and any other dynamic variables.

t0

The time of origin of the process

sampleTimes

A numeric vector providing the times of samples. If named, taxon labels will be based on names of the corresponding sample times

sampleStates

For models with more than one deme, a matrix must be supplied describing the probability of sampling from each deme. Each row corresponds to a sample in the same order as sampleTimes. Each column corresponds to probability of sampling each deme. Column names should be defined which correspond to deme names in the model. Rows should sum to one.

res

Integer number of time steps to use when simulating model.

integrationMethod

If simulating an ODE (ordinary differential equation) model, this provides the integration routine corresponding to options in deSolve.

finiteSizeCorrections

to do

substitutionRates

to do

sequenceLength

to do

Value

An object of class DatedTree, which subclasses ape::phylo. The $heights attribute provides the time of each node in the tree. The $maxHeight attribute provides the time of the root.

Author

Erik Volz

Examples

# A simple exponential growth model with birth rates beta, and death rates gamma: 
# I is the number of infected individuals
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)
# plot the tree 
plot(tre)