cladeMode: Analyse nested modes of evolution
In PuttickMacroevolution/cladeMode: nested modes of evolution

Description Usage Arguments Value Examples

Analyses various nested modes of continous trait evolution applied to the whole phylogeny (BM, OU, EB) or nested within subclades of the phylogeny ("nested EB", "nested EB rate", "rate Shift", "nested OU")

1
2
3

cladeMode(phy, y, node = F, minSize = 5, mode = c("BM", "OU", "EB",
  "nestedEB", "nestedEBRate", "rateShift", "nestedOU"), Ncores = 1,
  contrasts = FALSE, returnPhy = FALSE, measureError = NULL)

`phy`	dated tree in APE format. Although non-ultrametric phylogenies can be used the method has only been tested on ultrametric trees
`y`	named vector of trait data
`node`	number of node at which to estimate nested evolutionary mode. If FALSE (default) the model tests model fit at all nodes using an iterative procedure. The nodes at which the model is applied it determined by the minSize argument
`minSize`	the minimum size of node at which to select a nested rate. The number represents the minimum size of the clade for which all models will be nested.
`mode`	of evolution for the model. One of "BM", "EB", "OU", "nestedEB", "nestedEBRate", "rateShift", "nestedOU"). "BM" fits a simple model of Brownian Motion; "EB" fits the early burst model (Harmon et al. 2010); "OU" fits the Orstein-Uhlenbeck model (Hansen). The 'nested' models allow shifts to occur in subclades of the phylogeny with an ancestral BM model. 'nested EB' allows for a EB model to a subclade; 'nested EB rate' fits a model a similar model but with the addition of a scalar that allows for an increase on the ancestral BM rate; "nested OU" fits a nested model of the OU process which optimises the alpha parameter; and the "rateShift" model allows for a scalar-based increase or decrease rate in a nested subclade (O'Meara et al. 2006; Thomas et al. 2006)
`contrasts`	logical (default = FALSE) as to whether to estimate parameters with the default likelihood search or independent contrasts
`returnPhy`	logical (default = FALSE) returns re-scaled phylogeny according the best-fitting model
`measureError`	a named vector of measurement error around the trait data y. If NULL (default) the model assumes no measurement error. Warning: this is currently not fully tested and only applies to all models (except BM) when 'contrasts=FALSE'
`ncores`	number of seperate cores to simultaneously estimate rate shifts - defaults to 1

results a list containing the model parameters, log likelihood, aicc, optionally the rescaled tree

library(geiger)
set.seed(30)
# simulate tree with 20 species
simTree <- rcoal(20)
simData <- cladeModeSim(n=1, phy=simTree, mode="BM")

##### comparison between geiger and cladeMode parameter estimates using early burst model

#cladeMode

unlist(cladeMode(phy=simTree, y=simData[,1], mode="EB", cont=T))
#cladeMode paramater estimates 

#	lnl		root.state	beta	a	aic	aicc	k 
#	-13.12825138	-0.38104635	1.14123829	-0.06945831	32.25650276	33.75650276	3.00000000 

#fitContinuous
unlist(fitContinuous(simTree, simData[,1], model="EB")[[4]])
#	a	sigsq	z0	lnL	method	k	aic	aicc
#	-0.06945	1.14123	-0.3810463	-13.1282513	subplex	3	32.2565	33.7565

################### Comparison of models            
# example use of each model in the function cladeMode          
# bm  
bmRes <- cladeMode(phy=simTree, y=simData[,1], mode="BM", cont=T)
# early burst
ebRes <-cladeMode(phy=simTree, y=simData[,1], mode="EB", cont=T)
# nested early burst with concurrent rate increase
nestedEBRateRes <-cladeMode(phy=simTree, y=simData[,1], mode="nestedEBRate", cont=T)
# nested early burst 
nestedEBRes <-cladeMode(phy=simTree, y=simData[,1], mode="nestedEB", cont=T)
# nested OU
nestedOURes <-cladeMode(phy=simTree, y=simData[,1], mode="nestedOU", cont=T)
# nested rate shift
nestedRateShiftRes <-cladeMode(phy=simTree, y=simData[,1], mode="rateShift", cont=T)
#comparison of models to summarise the best relative fit
aicW <- aiccWeight(c(bmRes$aicc, ebRes$aicc, nestedEBRateRes$aicc, nestedEBRes$aicc, nestedOURes$aicc, nestedRateShiftRes$aicc))
rownames(aicW) <- c("bm", "eb", "nestedEBRate", "nestedEB", "nestedOU", "nestedRateShift")
aicW
# nested rate shift has highest, but erroneus support (data were simulated under BM)
# Error correction. We will simulate 100 datasets (1000 were used in the ms, 100 used here for brevity) under BM and identify the error rate for the nested rate shift model
hunderedSim <- cladeModeSim(n=100, simTree, mode="BM")
# get aicc of nested rate models
nestedRateAICcSim <- sapply(1:100, function(x) cladeMode(phy=simTree, y=hunderedSim[,x], mode="rateShift", cont=T)$aicc)
# get aicc of bm models
bmAICcSim <- sapply(1:100, function(x) cladeMode(phy=simTree, y=hunderedSim[,x], mode="BM", cont=T)$aicc)
#calculate the 95th quantile under which nested rate is wrongly support
cutOffAICc <- as.numeric(quantile(nestedRateAICcSim - bmAICcSim, 0.95))
# add this number to nested rate aicc score and bm receives true model support
aicW <- aiccWeight(c(bmRes$aicc, nestedRateShiftRes$aicc + cutOffAICc))
aicW