geiger (version 2.0.11)

load.rjmcmc: posterior samples from single or multiple MCMC runs

Description

generating a pooled and thinned sample of posterior estimates from rjmcmc

Usage

load.rjmcmc(x, phy = NULL, burnin = NULL, thin = NULL, ...)

Value

an object of class rjmcmc or rjmcmcmc

Arguments

x

a vector of directory names (character strings) in which samples generated by geiger are found

phy

a phylogenetic tree of class 'phylo' for which to compile results; if NULL, the tree stored within each generated rda file is used

burnin

the proportion of each chain to be removed as burnin

thin

an integer specifying the thinning interval for chains

...

arguments to be passed to other functions (has no effect in the present context)

Author

JM Eastman

Details

This function provides a means to compile results from at least a single MCMC run, and is especially useful in pooling of multiple independent runs. In cases where MCMC chains were sampled using a set of (different) trees, the argument phy can be used to summarize compiled results against a single summary tree. Branchwise samples are stored with unique edge identifiers (character strings) that are more reliable than the use of numeric node-identifiers (see digest for the function used to generate these unique edge labels).

If it is desired to run analyses for the same trait data across a set of trees (see Examples), it is strongly recommended that the set of trees be sent to rjmcmc.bm as a multiPhylo object (see read.tree).

See Also

rjmcmc.bm

Examples

# \donttest{
sal=get(data(caudata))
a<-sim<-sal$phy
bl=c(386,387,388,183,184,185,186)
mod=match(bl, sim$edge[,2])
sim$edge.length[mod]=sim$edge.length[mod]*64
dat=rTraitCont(sim)
while(1){
	b=a
	b$tip.label[183:186]=sample(b$tip.label[183:186])
	if(!all(a$tip.label==b$tip.label)) break()
}

trees=list(a=a,b=b, c=ladderize(a, right=TRUE), d=ladderize(a, right=FALSE))
class(trees)="multiPhylo"
rjmcmc.bm(trees, dat, ngen=1e3, type="rbm")

res=load.rjmcmc(paste("relaxedBM", names(trees), sep="."), phy=trees$d, burnin=0.25)
plot(res, par="shifts", show.tip=FALSE, edge.width=2.5)
# }