coalescentMCMC (version 0.4-4)

coalescentMCMC: Run MCMC for Coalescent Trees

Description

These are the main function of the package to run a Markov chain Monte Carlo (MCMC) to generate a set of trees which is returned with their likelihoods, the coalescent likelihoods and the respective parameter(s).

The logLik method returns the average log-likelihood of the coalescent model. AIC, BIC, and anova use this average log-likelihood.

Usage

coalescentMCMC(x, ntrees = 3000, model = "constant", tree0 = NULL,
               printevery = 100, degree = 1, nknots = 0,
               knot.times = NULL, moves = 1:6)
# S3 method for coalescentMCMC
logLik(object, ...)
# S3 method for coalescentMCMC
AIC(object, ..., k = 2)
# S3 method for coalescentMCMC
BIC(object, ...)
# S3 method for coalescentMCMC
anova(object, ...)

Arguments

x

a set of DNA sequences, typically an object of class "DNAbin" or "phyDat".

ntrees

the number of trees to output.

tree0

the initial tree of the chain; by default, a UPGMA tree with a JC69 distance is generated.

model

the coalescent model to be used for resampling. By default, a constant-THETA is used.

printevery

an integer specifying the frequency at which to print the numbers of trees proposed and accepted; set to 0 to cancel all printings.

degree, nknots, knot.times

parameters used if model = "splines".

moves

the tree moves used by the MCMC (see details).

…

options passed to other methods.

object

an bject of class "coalescentMCMC".

k

the coefficient used to calculate the AIC (see AIC).

Value

coalescentMCMC returns an object of class c("coalescentMCMC", "coda") with the log-likelihood and the parameters of each tree.

logLik, AIC and BIC return a numeric vector.

anova return an object of class "anova".

Details

Six tree moves are programmed and one is chosen randomly at each step of the MCMC. The steps are: (1) NeighborhoodRearrangement (Kuhner et al., 1995), (2) ScalingMove, (3) branchSwapping, (4) subtreeExchange, (5) NodeAgeMove, and (6) randomWalkThetaMu (all five from Drummond et al., 2002). In practice, it appears that in many situations moves = c(1, 3) is a good selection resulting in around 50% acceptance rate.

References

Drummond, A. J., Nicholls, G. K., Rodrigo, A. G. and Solomon, W. (2002) Estimating mutation parameters, population history and genealogy simultaneously from temporally spaced sequence data. Genetics, 161, 1307--1320.

Hastings, W. K. (1970) Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57, 97--109.

Kuhner, M. K., Yamato, J. and Felsenstein, J. (1995) Estimating effective population size and mutation rate from sequence data using Metropolis-Hastings sampling. Genetics, 140, 1421--1430.

See Also

getMCMCtrees, dcoal, treeOperators

Examples

# NOT RUN {
data(woodmouse)
out <- coalescentMCMC(woodmouse)
plot(out)
getMCMCtrees() # returns 3000 trees
# }