ape (version 3.1-3)

mcmc.popsize: Reversible Jump MCMC to Infer Demographic History

Description

These functions implement a reversible jump MCMC framework to infer the demographic history, as well as corresponding confidence bands, from a genealogical tree. The computed demographic history is a continous and smooth function in time. mcmc.popsize runs the actual MCMC chain and outputs information about the sampling steps, extract.popsize generates from this MCMC output a table of population size in time, and plot.popsize and lines.popsize provide utility functions to plot the corresponding demographic functions.

Usage

mcmc.popsize(tree,nstep, thinning=1, burn.in=0,progress.bar=TRUE,
    method.prior.changepoints=c("hierarchical", "fixed.lambda"), max.nodes=30,
   lambda=0.5, gamma.shape=0.5, gamma.scale=2,
    method.prior.heights=c("skyline", "constant", "custom"),
    prior.height.mean,
    prior.height.var)
extract.popsize(mcmc.out, credible.interval=0.95, time.points=200, thinning=1, burn.in=0)
## S3 method for class 'popsize':
plot(x, show.median=TRUE, show.years=FALSE, subst.rate, present.year, ...)
## S3 method for class 'popsize':
lines(x, show.median=TRUE,show.years=FALSE, subst.rate, present.year, ...)

Details

Please refer to Opgen-Rhein et al. (2005) for methodological details, and the help page of skyline for information on a related approach.

References

Opgen-Rhein, R., Fahrmeir, L. and Strimmer, K. 2005. Inference of demographic history from genealogical trees using reversible jump Markov chain Monte Carlo. BMC Evolutionary Biology, 5, 6.

See Also

skyline and skylineplot.

Examples

# get tree
data("hivtree.newick") # example tree in NH format
tree.hiv <- read.tree(text = hivtree.newick) # load tree

# run mcmc chain
mcmc.out <- mcmc.popsize(tree.hiv, nstep=100, thinning=1, burn.in=0,progress.bar=FALSE) # toy run
#mcmc.out <- mcmc.popsize(tree.hiv, nstep=10000, thinning=5, burn.in=500) # remove comments!!

# make list of population size versus time
popsize  <- extract.popsize(mcmc.out)

# plot and compare with skyline plot
sk <- skyline(tree.hiv)
plot(sk, lwd=1, lty=3, show.years=TRUE, subst.rate=0.0023, present.year = 1997)
lines(popsize, show.years=TRUE, subst.rate=0.0023, present.year = 1997)