rjmcmc.bm {geiger}

Bayesian sampling of shifts in trait evolution: relaxed Brownian motion

Description

Implements reversible-jump Markov chain Monte Carlo sampling for trait evolutionary models

Usage

rjmcmc.bm(phy, dat, SE=NA, ngen = 50000, samp = 100, 
    type = c("jump-rbm", "rbm", "jump-bm", "bm"), ...)

Arguments

phy

a phylogenetic tree of class 'phylo'

dat

a named vector of continuous trait values, associated with each species in phy

SE

a named vector of standard errors for each trait value; applied to all trait values if given a single value

ngen

number of sampling generations

samp

frequency with which Markov samples are retained (e.g., samp=10 retains every tenth sample in the chain)

type

the class of model to use (see Details)

...

arguments passed to make.gbm

Details

Implemented is an MCMC sampler for a general model of Brownian motion, which in the full model (type="jump-rbm") allows relaxed local clocks and also a point process of pulses in evolutionary rate along individual branches. Restricted models include global-rate Brownian motion (type="bm"), relaxed-rates Brownian motion (type="rbm"), and models including jumps but a single rate of diffusion across the tree (type="jump-bm").

Where applicable, posterior estimates of shifts between local rates, estimates of the rates themselves, and inferred jumps (or pulses) are provided as output. Estimates are stored as an MCMC-generations-by-branches matrix (see Examples), and branches are uniquely labeled by a cryptographic function to ensure comparability amongst trees differing in topology (see digest).

Note that default settings (as the user assumes if nothing is specified in ...) provide absolutely no guarantee of the chain achieving convergence. The user is emphatically encouraged to supply informed arguments for what are the most critical aspects of this MCMC sampler (see make.gbm for more information on permissible modifications to the MCMC sampler). Finding reasonable run parameters will likely require much trial and error. Run diagnosis and inspection of chain mixing is facilitated by the R-package coda or by the Java application, Tracer (http://tree.bio.ed.ac.uk/software/tracer/).

In the Examples below, do not expect such short chains to reach stationarity!

Value

After a run has completed, acceptance rates for the primary proposal mechanisms are printed to the console, along with settings of control parameters for the run (see make.gbm).

Posterior results are written to several files within a base directory, the contents of which are as follows:

\bold{log}

is a logfile including the following for each Markov chain: the generations at which samples were retained (state), the min, max, and median rate of the diffusion process across the tree, the number of evolutionary pulses (jumps) along single branches, the variance associated with the jump process (jumpvar), the root state, and the likelihood (lnL) and prior (lnLp) of the model at sampled generations.

\bold{rda}

is a compressed R object, which stores branchwise estimates of the jump and diffusion processes. In order to be interpretable, the rda file should be processed by the function load.rjmcmc. The package coda can be used within R for run diagnostics on the processed output (see, e.g., heidel.diag and autocorr.

Author(s)

JM Eastman, LJ Harmon, AL Hipp, and JC Uyeda

References

Eastman JM, ME Alfaro, P Joyce, AL Hipp, and LJ Harmon. 2011. A novel comparative method for identifying shifts in the rate of character evolution on trees. Evolution 65:3578-3589.

See Also

load.rjmcmc

Examples

## GENERATE DATA: jump-diffusion
phy <- ladderize(sim.bdtree(n=200), right=FALSE)
r <- paste(sample(letters,9,replace=TRUE),collapse="")
defpar <- par(no.readonly=TRUE)


tmp <- ex.jumpsimulator(phy, jumps=10)
dat <- tmp$dat
hist <- tmp$hist

ex.traitgram(phy, hist, alpha=0) # plot history of trait change

## RUN ANALYSIS


## coda package is not a dependency of geiger
## but is very useful for evaluating mcmc runs
## library(coda)


rjmcmc.bm(phy,dat, prop.width=1.5, ngen=20000, samp=500, filebase=r,
          simple.start=TRUE, type="jump-bm")
outdir <- paste("jump-BM", r, sep=".")
ps <- load.rjmcmc(outdir)

dev.new()
plot(x=ps, par="jumps", burnin=0.25, legend=FALSE, show.tip=FALSE, type="fan", edge.width=2)
mm=match(phy$edge[,2],hist$descendant)
hist=hist[mm,]
edgelabels.auteur(text=NULL, pch=21, cex=hist$cex, bg=NA, col=ifelse(hist$cex>0, 1, NA), lty=2)
title("red (estimated); black (true jump size)", line=-5)
par(defpar)

dev.new()
## from the coda package
coda::autocorr.plot(ps$log, ask=dev.interactive())
plot(ps$log, ask=dev.interactive())

## GENERATE DATA: multi-rate diffusion
scl <- ex.ratesimulator(phy, min=12, show.tip=FALSE)
dat <- rTraitCont(scl)

## RUN ANALYSIS
rjmcmc.bm(phy, dat, prop.width=1.5, ngen=20000, samp=500, filebase=r, simple.start=TRUE, type="rbm")
outdir <- paste("relaxedBM", r, sep=".")
ps <- load.rjmcmc(outdir)
dev.new()
plot(x=ps, par="shifts", burnin=0.25, legend=TRUE, show.tip=FALSE, edge.width=2)

if(!interactive()) { ## clean up
    dirs <- dir(pattern="^(jump-BM|relaxedBM)")
    unlink(dirs, recursive=TRUE)
}