3 \alias{extract.popsize}
6 \title{Reversible Jump MCMC to Infer Demographic History}
8 mcmc.popsize(tree,nstep, thinning=1, burn.in=0,progress.bar=TRUE,
9 method.prior.changepoints=c("hierarchical", "fixed.lambda"), max.nodes=30,
10 lambda=0.5, gamma.shape=0.5, gamma.scale=2,
11 method.prior.heights=c("skyline", "constant", "custom"),
14 extract.popsize(mcmc.out, credible.interval=0.95, time.points=200, thinning=1, burn.in=0)
15 \method{plot}{popsize}(x, show.median=TRUE, show.years=FALSE, subst.rate, present.year, ...)
16 \method{lines}{popsize}(x, show.median=TRUE,show.years=FALSE, subst.rate, present.year, ...)
20 \item{tree}{Either an ultrametric tree (i.e. an object of class \code{"phylo"}),
21 or coalescent intervals (i.e. an object of class \code{"coalescentIntervals"}). }
22 \item{nstep}{Number of MCMC steps, i.e. length of the Markov chain (suggested value: 10,000-50,000).}
23 \item{thinning}{Thinning factor (suggest value: 10-100).}
24 \item{burn.in}{Number of steps dropped from the chain to allow for a burn-in phase (suggest value: 1000).}
26 \item{progress.bar}{Show progress bar during the MCMC run.}
28 \item{method.prior.changepoints}{If \code{hierarchical}is chosen (the default) then the smoothing parameter lambda is drawn from
29 a gamma distribution with some specified shape and scale parameters.
30 Alternatively, for \code{fixed.lambda} the value of lambda is a given constant.
33 \item{max.nodes}{Upper limit for the number of internal nodes of the approximating spline (default: 30).}
34 \item{lambda}{Smoothing parameter. For \code{method="fixed.lambda"} the specifed value of lambda determines
35 the mean of the prior distribution for the number of internal nodes of the approximating
36 spline for the demographic function (suggested value: 0.1-1.0).}
37 \item{gamma.shape}{Shape parameter of the gamma function from which \code{lambda} is drawn for
38 \code{method="hierarchical"}.}
39 \item{gamma.scale}{Scale parameter of the gamma function from which \code{lambda} is drawn for
40 \code{method="hierarchical"}.}
41 \item{method.prior.heights}{Determines the prior for the heights of the change points.
42 If \code{custom} is chosen then two functions describing the mean and variance
43 of the heigths in depence of time have to be specified (via \code{prior.height.mean}
44 and \code{prior.height.var} options). Alternatively, two built-in priors are available:
45 \code{constant} assumes constant population size and variance determined by Felsenstein
46 (1992), and \code{skyline} assumes a skyline plot (see Opgen-Rhein et al. 2004 for
48 \item{prior.height.mean}{Function describing the mean of the prior distribution for the heights
49 (only used if \code{method.prior.heights = custom}).}
51 \item{prior.height.var}{Function describing the variance of the prior distribution for the heights
52 (only used if \code{method.prior.heights = custom}).}
53 \item{mcmc.out}{Output from \code{mcmc.popsize} - this is needed as input for \code{extract.popsize}.}
54 \item{credible.interval}{Probability mass of the confidence band (default: 0.95).}
56 \item{time.points}{Number of discrete time points in the table output by \code{extract.popsize}.}
58 \item{x}{Table with population size versus time, as computed by \code{extract.popsize}. }
60 \item{show.median}{Plot median rather than mean as point estimate for demographic function (default: TRUE).}
62 \item{show.years}{Option that determines whether the time is plotted in units of
63 of substitutions (default) or in years (requires specification of substution rate
64 and year of present).}
65 \item{subst.rate}{Substitution rate (see option show.years).}
66 \item{present.year}{Present year (see option show.years).}
67 \item{\dots}{Further arguments to be passed on to \code{plot}.}
70 These functions implement a reversible jump MCMC framework to infer the demographic history,
71 as well as corresponding confidence bands,
72 from a genealogical tree. The computed demographic history is a continous
73 and smooth function in time.
74 \code{mcmc.popsize} runs the actual MCMC chain and outputs information about the
75 sampling steps, \code{extract.popsize} generates from this MCMC
76 output a table of population size in time, and \code{plot.popsize} and \code{lines.popsize}
77 provide utility functions to plot the corresponding demographic functions.
81 Please refer to Opgen-Rhein et al. (2005) for methodological details, and the help page of
82 \code{\link{skyline}} for information on a related approach.
86 \author{Rainer Opgen-Rhein (\url{http://www.stat.uni-muenchen.de/~opgen/}) and
87 Korbinian Strimmer (\url{http://www.stat.uni-muenchen.de/~strimmer/}).
88 Parts of the rjMCMC sampling procedure are adapted from R code by Karl Browman
89 (\url{http://www.biostat.jhsph.edu/~kbroman/})}
92 \code{\link{skyline}} and \code{\link{skylineplot}}. }
94 Opgen-Rhein, R., Fahrmeir, L. and Strimmer, K. 2005. Inference of
95 demographic history from genealogical trees using reversible jump
96 Markov chain Monte Carlo. \emph{BMC Evolutionary Biology}, \bold{5},
101 data("hivtree.newick") # example tree in NH format
102 tree.hiv <- read.tree(text = hivtree.newick) # load tree
105 mcmc.out <- mcmc.popsize(tree.hiv, nstep=100, thinning=1, burn.in=0,progress.bar=FALSE) # toy run
106 #mcmc.out <- mcmc.popsize(tree.hiv, nstep=10000, thinning=5, burn.in=500) # remove comments!!
108 # make list of population size versus time
109 popsize <- extract.popsize(mcmc.out)
111 # plot and compare with skyline plot
112 sk <- skyline(tree.hiv)
113 plot(sk, lwd=1, lty=3, show.years=TRUE, subst.rate=0.0023, present.year = 1997)
114 lines(popsize, show.years=TRUE, subst.rate=0.0023, present.year = 1997)