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

In geiger: Analysis of Evolutionary Diversification

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

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)

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

Value

an object of class rjmcmc or rjmcmcmc

Author(s)

JM Eastman

See Also

rjmcmc.bm

Examples

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)

geiger documentation built on April 4, 2023, 1:07 a.m.