mecca: running a MECCA analysis

In geiger: Analysis of Evolutionary Diversification

running a MECCA analysis

Description

Runs MECCA's hybrid ABC-MCMC algorithm to jointly estimate diversification rates and trait evolution from incompletely sampled comparative data. Many of the arguments taken by this function are the same as those in calibrateMecca().

Usage

mecca(phy, richness, cladeMean, cladeVariance, model = c("BM", "Trend", "twoRate"),
   prior.list = list(priorSigma = c(-4.961845, 4.247066), priorMean = c(-10, 10)),
   start = start, Ngens = 10000, printFreq = 100, sigmaPriorType = "uniform",
   rootPriorType = "uniform", SigmaBounds = c(-4.961845, 4.247066), hotclade = NULL,
   divPropWidth = 0.1, scale = 1, divSampleFreq = 0, BoxCox = TRUE, outputName = "mecca")

Arguments

phy	time-calibrated phylogenetic tree of class 'phylo'
richness	named vector of species richnesses corresponding to tips in the tree (if a tip does not represent a higher level taxon, its richness should be 1)
cladeMean	named vector of trait means (all tips in the tree must be represented)
cladeVariance	named vector of trait variances (all tips in the tree must be represented; if only one taxon is present, use 0 for the variance)
model	model of trait evolution to be used – options currently implemented are "BM" = Brownian Motion, "Trend" = Brownian moton with a trend, and "twoRate" = two Brownian rate model (see hotclade below).
prior.list	a list containing prior distribution parameters (default values used if this argument is empty)
start	ouput of startingpt.mecca
Ngens	number of generations to run MECCA
printFreq	frequency of printing the acceptance rate
sigmaPriorType	type of prior distribution on the Brownian rate parameter (currently either "uniform" or "normal")
rootPriorType	type of prior distribution on the root state rate parameter (currently "uniform" is available)
SigmaBounds	bounds for sigma (default values correspond to a wide range taken from Harmon et al. 2010)
hotclade	if a two-rate model is to be fit, this specifies which clade takes the second rate – two names should be specified in a vector; either the two tip names that span the clade of interest, or the name of a terminal/internal edge and NULL if only one branch takes the second rate
divPropWidth	proposal width for the diversification MCMC (default value of 0.1 seems to work well)
scale	a numeric value by which the proposal width for trait evolution parameters will be multiplied (a value of 2 seems to work well, but this should be adjusted for each individual dataset)
divSampleFreq	whether new trees are simulated at every step – the default (0) is yes; if a non-zero value is given, this will determine the frequency (every n steps) with which new tip trees are simulated
BoxCox	whether summaries are BOX-COX standardized – default (1) is yes and is recommended; this should always be consistent with the calibration step
outputName	name stem for output file names

Details

The output files produced are formatted to be used with the C++ Program ABCtoolbox (Wegmann et al. 2010), which produces adjusted posterior distributions and can perform model selection without likelihoods.

Value

MECCA does not store any output in memory. Instead, five output files are generated to the current working directory. This files are fully compatible with ABCtoolbox (Wegmann et al. 2011). The first file (outputname_bdSimFile.tx) will output the posterior sample for diversification parameters. The second file (outputname_bmSimFile.txt) ouputs the sampled trait evolution parameters and their associated raw summary statistics while outputname_ObsFile.txt gives the observed summaries.For ABC toolbox though, it will often be more efficient to use pls-transformed versions of the observed and simulated summary statistics. These are available in outputname_distObs.txt, and outputname_distSimFile.txt.

Note

The numbers printed to the screen during the run give the current generation, acceptance rate for diversification parameters, acceptance rate for trait evolutionary rate parameters and acceptance rate for root state parameters, respectively

Author(s)

Graham Slater, Luke Harmon, Daniel Wegmann

References

Slater GJ, LJ Harmon, D Wegmann, P Joyce, LJ Revell, and ME Alfaro. 2012. Fitting models of continuous trait evolution to incompletely sampled comparative data using approximate Bayesian computation. Evolution 66:752-762.

Examples

## Not run: 
data(carnivores)
phy <- carnivores$phy
data <- carnivores$dat
richness <- data[,1]
names(richness) <- rownames(data)

priors <- list(priorSigma = c(-4.5, 4.5), priorMean = c(-5, 2))

## CALIBRATION (far too short for a real analysis)
Cal <- calibrate.mecca(phy, richness, model = "BM", prior.list = priors, Ncalibrations = 1000)

params <- Cal$trait[, c(1,2)] ## extract the calibration BM parameters
stats <- Cal$trait[, -c(1,2)] ## extract the calibration summary stats

## now we run pls, determining combinations of summaries that explain variation in our parameters
## For BM, 2 components is sufficient. For more complex models, more componenets will be required.
require(pls)
myPlsr<-pls::plsr(as.matrix(params) ~ as.matrix(stats), scale=F, ncomp = 2)

plot(RMSEP(myPlsr)) ## Look at Root Mean Square error plots

summary(myPlsr) ## take a look at 

plsdat <- myPlsr$loadings

## extract means and variances from the carnivore data ##

cladeMean<-data[,2]
names(cladeMean)<-rownames(data)
cladeVariance<-data[,3]
names(cladeVariance)<-rownames(data)

## STARTING POINT
## And now we can compute starting values for the ABC-MCMC
start <- startingpt.mecca(Cal, phy, cladeMean, cladeVariance,
tolerance = 0.05, plsdat, BoxCox = TRUE)


## MECCA (far too short for a real analysis)
mecca(phy, richness, cladeMean, cladeVariance, model = "BM", prior.list = priors, start = start,
   Ngens = 1000, printFreq = 100, sigmaPriorType = "uniform", rootPriorType = "uniform",
   SigmaBounds = c(-4.5, 4.5), divPropWidth = 0.1, scale = 2, divSampleFreq = 0, BoxCox = TRUE,
   outputName ="MeccaBM.txt")

## PASTE UNCOMMENTED FOLLOWING LINE TO DROP FILES CREATED BY MECCA
 # unlink(dir(pattern=paste(r)),recursive=TRUE)


## End(Not run)

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