arbutus | R Documentation |
Arbutus was developed to evaluate the adequacy (or absolute goodness-of-fit) of phylogenetic models of continuous character evolution. The basic principle underlying the approach is that if the generating model is a Brownian motion process, the contrasts (sensu Felsenstein, 1985) will be independent and identically distributed (I.I.D.). We can evaluate this condition by calculating a set of test statistics the contrasts on our observed data, then simulating datasets under a Brownian motion process along the phylogeny and calculating the test statistics on the contrasts of each simulated data set. We can then compare our observed test statistics to the simulated distribution of test statistics.
The I.I.D. property of the contrasts will hold only for the case of Brownian motion. However, the approach can be extended to any arbitrarily complex model of continuous character evolution (as long as the model is based on a multivariate normal distribution) by first rescaling the phylogeny based on the fitted parameter estimates, creating what we refer to as a 'unit.tree'. Once the rescaling is done, the contrasts will be I.I.D. if the model that was fit is a adequate one.
The basic functionality of arbutus can be employed using the function arbutus
.
This function wraps several other functions which can be used separately:
make_unit_tree
: rescales phylogeny based on fitted parameter values
calculate_pic_stat
: calculates test statistics on observed data
simulate_char_unit
: simulates datasets under Brownian motion
calculate_pic_stat
: calculates test statistics on simulated data
compare_pic_stat
: compares observed to simulated statistics
Arbutus is designed to interact with a number of other packages for fitting trait models. The architecture is flexible such that additional packages can be incorporated relatively easily.
Use test statistics to assess model adequacy for phylogenetic models of continuous character evolution. This function is a simply a wrapper which does all the steps involved in evaluating model adequacy using the approach outlined in Pennell, FitzJohn, Cornwell and Harmon (in review).
arbutus(x, nsim = 1000, stats = NULL, ...)
x |
a fitted model object or a |
nsim |
the number of datasets to simulate. This is passed as an argument to the
function |
stats |
a named list of test statistics to calculate on observed and simulated
datasets. See |
... |
additional arguments to be passed to |
arbutus
was developed to evaluate the adequacy (or absolute goodness-of-fit)
of phylogenetic models of continuous character evolution. The basic principle
underlying the approach is that if the generating model is a Brownian motion
process, the contrasts (sensu Felsenstein, 1985) will be independent and
identically distributed (I.I.D.). We can evaluate this condition by calculating
a set of test statistics on the contrasts of our observed data, then simulating datasets under
a Brownian motion process along the phylogeny and calculating the test statistics
on the contrasts of each simulated data set. We can then compare our observed test statistics to the
simulated distribution of summary statistics.
The I.I.D. property of the contrasts will hold only for the case of Brownian motion. However, the approach can be extended to any arbitrarily complex model of continuous character evolution (as long as the model is based on a multivariate normal distribution) by first rescaling the phylogeny based on the fitted parameter estimates, creating what we refer to as a 'unit.tree'. Once the rescaling is done, the contrasts will be I.I.D. if the model that was fit is adequate.
P-values represent the two-tailed probability that the observed summary statistic (or distribution of summary statistics) came from the same distribution as the simulated summary statistics. Low p-values provide evidence that the model is inadequate.
The function arbutus
wraps several other functions:
make_unit_tree
: rescales phylogeny based on fitted parameter values
calculate_pic_stat
: calculates test statistics on observed data
simulate_char_unit
: simulates datasets under Brownian motion with rate 1
calculate_pic_stat
: calculates test statistics on simulated data
compare_pic_stat
: compares observed to simulated test statistics
arbutus
interacts with objects produced by fitting evolutionary models using a variety of packages.
Currently supported objects are as follows:
a gfit
object returned from fitting a model of continuous character evolution using
fitContinuous
in the geiger
package.
a fit.mle
object returned from fitting a model of continuous character evolution
using find.mle
in the diversitree
package.
a mcmcsamples
object returned from fitting a model of continuous character evolution
using MCMC methods in the diversitree
package.
make_unit_tree
will apply the same
trait dataset to a set of unit trees based on sampled parameters. By default this will create a
unit tree for every sample in the mcmc chain. To modify this, additional arguments can be use.
burnin
specifies how many samples to remove from the beginning of the chain.
thin
specifies the thinning interval (e.g. if thin=5
, the function will create a unit
tree from every fifth parameter set sampled.
sample
specifies how many samples to draw from the MCMC run.
a gls
object returned from fitting a phylogenetic generalized least squares model
of character correlation using gls
in the nlme
package.
a pgls
object returned from fitting a phylogenetic generalized least squares model
of character correlation using pgls
in the caper
package.
a phylolm
object returned from fitting a phylogenetic generalized linear model of
character correlation using phylolm
in the phylolm
package. As the phylogeny is not
returned with the phylolm
object, a phy
argument must also be specified.
a phylo
object. If a phylo
object is supplied, the tree is assumed to have been
rescaled previously. A second argument data
must also be provided included the trait
data as a named vector with names equal to the tip.labels of the phylogeny.
a multiPhylo object
. If a multiPhylo
object is supplied, the tree is assumed to have been
rescaled previously. A second argument data
must also be provided included the trait
data as a named vector with names equal to the tip.labels of the phylogenies. Note that this
function will append the same data set to every tree in the multiPhylo
object.
a OUwie
object returned from fitting a model of continuous character evolution using
OUwie
in the OUwie
package.
Matthew W. Pennell, Richard G. FitzJohn, William K. Cornwell and Luke J. Harmon (in review.). Model adequacy and the macroevolution of angiosperm functional traits.
Joseph Felsenstein (1985). Phylogenies and the comparative method. The American Naturalist 125:1-15.
make_unit_tree
, calculate_pic_stat
, simulate_char_unit
, compare_pic_stat
## finch data data(finch) phy <- finch$phy data <- finch$data[,"wingL"] ## Not run: require(geiger) ## fit Brownian motion model ## using geiger's fitContinuous function fit.bm <- fitContinuous(phy=phy, dat=data, model="BM", control=list(niter=10), ncores=1) ## check adequacy of BM model ## get p-values for default summary statistics modelad.bm <- arbutus(fit.bm, nsim=10) modelad.bm ## fit Ornstein-Uhlenbeck model ## again, using geiger's fitContinuous function fit.ou <- suppressWarnings(fitContinuous(phy=phy, dat=data, model="OU", control=list(niter=10), ncores=1)) ## check adequacy of OU model modelad.ou <- arbutus(fit.ou, nsim=10) require(diversitree) ## fit Brownian motion model using ML ## using diversitree's find.mle function bmlik <- make.bm(phy, data) fit.bm.dt <- find.mle(bmlik, 0.1) ## this creates a 'fit.mle' object which can be used ## in 'arbutus' modelad.bm.dt <- arbutus(fit.bm.dt) ## fit Brownian motion model using MCMC mcmc.bm.dt <- mcmc(bmlik, x.init=1, nsteps=1000, w=1) ## construct a unit tree object from mcmcsamples ## removing 100 samples as burnin ## and sampling 10 parameter sets modelad.bm.dt.mcmc <- arbutus(mcmc.bm.dt, burnin=100, samples=10) require(nlme) ## Use pgls to look for a correlation between two traits t1 <- data t2 <- td$data[,"tarsusL"] dd <- cbind.data.frame(t1, t2) ## fit gls model with corPagel correlation structure fit.gls <- gls(t1~t2, data=dd, correlation=corPagel(phy=phy, value=1)) ## this creates a 'gls' object which can be used ## in 'arbutus' modelad.gls <- arbutus(fit.gls) ## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.