ape: mcmc.popsize – R documentation

Become an expert in R — Interactive courses, Cheat Sheets, certificates and more!

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, 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, ...)

Arguments

`tree`	Either an ultrametric tree (i.e. an object of class `"phylo"`), or coalescent intervals (i.e. an object of class `"coalescentIntervals"`).
`nstep`	Number of MCMC steps, i.e. length of the Markov chain (suggested value: 10,000-50,000).
`thinning`	Thinning factor (suggest value: 10-100).
`burn.in`	Number of steps dropped from the chain to allow for a burn-in phase (suggest value: 1000).
`progress.bar`	Show progress bar during the MCMC run.
`method.prior.changepoints`	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.
`max.nodes`	Upper limit for the number of internal nodes of the approximating spline (default: 30).
`lambda`	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).
`gamma.shape`	Shape parameter of the gamma function from which `lambda` is drawn for `method="hierarchical"`.
`gamma.scale`	Scale parameter of the gamma function from which `lambda` is drawn for `method="hierarchical"`.
`method.prior.heights`	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).
`prior.height.mean`	Function describing the mean of the prior distribution for the heights (only used if `method.prior.heights = custom`).
`prior.height.var`	Function describing the variance of the prior distribution for the heights (only used if `method.prior.heights = custom`).
`mcmc.out`	Output from `mcmc.popsize` - this is needed as input for `extract.popsize`.
`credible.interval`	Probability mass of the confidence band (default: 0.95).
`time.points`	Number of discrete time points in the table output by `extract.popsize`.
`x`	Table with population size versus time, as computed by `extract.popsize`.
`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 x-axis (depends on the value of `show.years`).
`ylab`	label on the y-axis.
`log`	log-transformation of axes; by default, the y-axis is log-transformed.
`...`	Further arguments to be passed on to `plot` or `lines`.

Details

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

Author(s)

Rainer Opgen-Rhein and Korbinian Strimmer. Parts of the rjMCMC sampling procedure are adapted from R code by Karl Broman.

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.

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)

ape

Analyses of Phylogenetics and Evolution

v5.5

GPL-2 | GPL-3

Authors

Emmanuel Paradis [aut, cre, cph] (<https://orcid.org/0000-0003-3092-2199>), Simon Blomberg [aut, cph] (<https://orcid.org/0000-0003-1062-0839>), Ben Bolker [aut, cph] (<https://orcid.org/0000-0002-2127-0443>), Joseph Brown [aut, cph] (<https://orcid.org/0000-0002-3835-8062>), Santiago Claramunt [aut, cph] (<https://orcid.org/0000-0002-8926-5974>), Julien Claude [aut, cph] (<https://orcid.org/0000-0002-9267-1228>), Hoa Sien Cuong [aut, cph], Richard Desper [aut, cph], Gilles Didier [aut, cph] (<https://orcid.org/0000-0003-0596-9112>), Benoit Durand [aut, cph], Julien Dutheil [aut, cph] (<https://orcid.org/0000-0001-7753-4121>), RJ Ewing [aut, cph], Olivier Gascuel [aut, cph], Thomas Guillerme [aut, cph] (<https://orcid.org/0000-0003-4325-1275>), Christoph Heibl [aut, cph] (<https://orcid.org/0000-0002-7655-3299>), Anthony Ives [aut, cph] (<https://orcid.org/0000-0001-9375-9523>), Bradley Jones [aut, cph] (<https://orcid.org/0000-0003-4498-1069>), Franz Krah [aut, cph] (<https://orcid.org/0000-0001-7866-7508>), Daniel Lawson [aut, cph] (<https://orcid.org/0000-0002-5311-6213>), Vincent Lefort [aut, cph], Pierre Legendre [aut, cph] (<https://orcid.org/0000-0002-3838-3305>), Jim Lemon [aut, cph], Guillaume Louvel [aut, cph] (<https://orcid.org/0000-0002-7745-0785>), Eric Marcon [aut, cph] (<https://orcid.org/0000-0002-5249-321X>), Rosemary McCloskey [aut, cph] (<https://orcid.org/0000-0002-9772-8553>), Johan Nylander [aut, cph], Rainer Opgen-Rhein [aut, cph], Andrei-Alin Popescu [aut, cph], Manuela Royer-Carenzi [aut, cph], Klaus Schliep [aut, cph] (<https://orcid.org/0000-0003-2941-0161>), Korbinian Strimmer [aut, cph] (<https://orcid.org/0000-0001-7917-2056>), Damien de Vienne [aut, cph] (<https://orcid.org/0000-0001-9532-5251>)

Initial release

2021-04-24