mcmc.popsize {ape}  R Documentation 
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.
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, xlab = NULL, ylab = "Effective population size", log = "y", ...) ## S3 method for class 'popsize' lines(x, show.median=TRUE,show.years=FALSE, subst.rate, present.year, ...)
tree 
Either an ultrametric tree (i.e. an object of class 
nstep 
Number of MCMC steps, i.e. length of the Markov chain (suggested value: 10,00050,000). 
thinning 
Thinning factor (suggest value: 10100). 
burn.in 
Number of steps dropped from the chain to allow for a burnin phase (suggest value: 1000). 
progress.bar 
Show progress bar during the MCMC run. 
method.prior.changepoints 
If 
max.nodes 
Upper limit for the number of internal nodes of the approximating spline (default: 30). 
lambda 
Smoothing parameter. For 
gamma.shape 
Shape parameter of the gamma function from which 
gamma.scale 
Scale parameter of the gamma function from which 
method.prior.heights 
Determines the prior for the heights of the change points.
If 
prior.height.mean 
Function describing the mean of the prior distribution for the heights
(only used if 
prior.height.var 
Function describing the variance of the prior distribution for the heights
(only used if 
mcmc.out 
Output from 
credible.interval 
Probability mass of the confidence band (default: 0.95). 
time.points 
Number of discrete time points in the table output by 
x 
Table with population size versus time, as computed by 
show.median 
Plot median rather than mean as point estimate for demographic function (default: TRUE). 
show.years 
Option that determines whether the time is plotted in units of of substitutions (default) or in years (requires specification of substution rate and year of present). 
subst.rate 
Substitution rate (see option show.years). 
present.year 
Present year (see option show.years). 
xlab 
label on the xaxis (depends on the value of

ylab 
label on the yaxis. 
log 
logtransformation of axes; by default, the yaxis is logtransformed. 
... 
Further arguments to be passed on to 
Please refer to OpgenRhein et al. (2005) for methodological details, and the help page of
skyline
for information on a related approach.
Rainer OpgenRhein and Korbinian Strimmer. Parts of the rjMCMC sampling procedure are adapted from R code by Karl Broman.
OpgenRhein, 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.
skyline
and skylineplot
.
# 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)