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 popsize
plot(x, show.median=TRUE, show.years=FALSE,
subst.rate, present.year, xlab = NULL,
ylab = "Effective population size", log = "y", ...)
# S3 method for popsize
lines(x, show.median=TRUE,show.years=FALSE, subst.rate, present.year, ...)
Either an ultrametric tree (i.e. an object of class "phylo"
),
or coalescent intervals (i.e. an object of class "coalescentIntervals"
).
Number of MCMC steps, i.e. length of the Markov chain (suggested value: 10,000-50,000).
Thinning factor (suggest value: 10-100).
Number of steps dropped from the chain to allow for a burn-in phase (suggest value: 1000).
Show progress bar during the MCMC run.
If hierarchical
is chosen (the default) then the smoothing parameter lambda is drawn from
a gamma distribution with some specified shape and scale parameters.
Alternatively, for fixed.lambda
the value of lambda is a given constant.
Upper limit for the number of internal nodes of the approximating spline (default: 30).
Smoothing parameter. For method="fixed.lambda"
the specifed value of lambda determines
the mean of the prior distribution for the number of internal nodes of the approximating
spline for the demographic function (suggested value: 0.1-1.0).
Shape parameter of the gamma function from which lambda
is drawn for
method="hierarchical"
.
Scale parameter of the gamma function from which lambda
is drawn for
method="hierarchical"
.
Determines the prior for the heights of the change points.
If custom
is chosen then two functions describing the mean and variance
of the heigths in depence of time have to be specified (via prior.height.mean
and prior.height.var
options). Alternatively, two built-in priors are available:
constant
assumes constant population size and variance determined by Felsenstein
(1992), and skyline
assumes a skyline plot (see Opgen-Rhein et al. 2004 for
more details).
Function describing the mean of the prior distribution for the heights
(only used if method.prior.heights = custom
).
Function describing the variance of the prior distribution for the heights
(only used if method.prior.heights = custom
).
Output from mcmc.popsize
- this is needed as input for extract.popsize
.
Probability mass of the confidence band (default: 0.95).
Number of discrete time points in the table output by extract.popsize
.
Table with population size versus time, as computed by extract.popsize
.
Plot median rather than mean as point estimate for demographic function (default: TRUE).
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).
Substitution rate (see option show.years).
Present year (see option show.years).
label on the x-axis (depends on the value of
show.years
).
label on the y-axis.
log-transformation of axes; by default, the y-axis is log-transformed.
Further arguments to be passed on to plot
or lines
.
Rainer Opgen-Rhein and Korbinian Strimmer. Parts of the rjMCMC sampling procedure are adapted from R code by Karl Broman.
Please refer to Opgen-Rhein et al. (2005) for methodological details, and the help page of
skyline
for information on a related approach.
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.
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)
Run the code above in your browser using DataLab