man/mcmc.popsize.Rd

   1 \name{mcmc.popsize}
   2 \alias{mcmc.popsize}
   3 \alias{extract.popsize}
   4 \alias{plot.popsize}
   5 \alias{lines.popsize}
   6 \title{Reversible Jump MCMC to Infer Demographic History}
   7 \usage{
   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"),
  12     prior.height.mean,
  13     prior.height.var)
  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, ...)
  17
  18 }
  19 \arguments{
  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).}
  25
  26   \item{progress.bar}{Show progress bar during the MCMC run.}
  27
  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.
  31   }
  32
  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
  47           more details).}
  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}).}
  50
  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).}
  55
  56  \item{time.points}{Number of discrete time points in the table output by \code{extract.popsize}.}
  57
  58  \item{x}{Table with population size versus time, as computed by \code{extract.popsize}. }
  59
  60  \item{show.median}{Plot median rather than mean as point estimate for demographic function (default: TRUE).}
  61
  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}.}
  68 }
  69 \description{
  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.
  78 }
  79
  80 \details{
  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.
  83 }
  84
  85
  86 \author{Rainer Opgen-Rhein and
  87         Korbinian Strimmer (\url{http://strimmerlab.org}).
  88         Parts of the rjMCMC sampling procedure are adapted from R code by Karl Browman
  89          (\url{http://www.biostat.wisc.edu/~kbroman/})}
  90
  91 \seealso{
  92 \code{\link{skyline}} and \code{\link{skylineplot}}. }
  93 \references{
  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},
  97   6.
  98 }
  99 \examples{
 100 # get tree
 101 data("hivtree.newick") # example tree in NH format
 102 tree.hiv <- read.tree(text = hivtree.newick) # load tree
 103
 104 # run mcmc chain
 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!!
 107
 108 # make list of population size versus time
 109 popsize  <- extract.popsize(mcmc.out)
 110
 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)
 115 }
 116 \keyword{manip}