fitting macroevolutionary models to phylogenetic trees
fitContinuous(phy, dat, SE = 0,
model = c("BM","OU","EB","rate_trend","lambda","kappa","delta","mean_trend","white"),
bounds= list(), control = list(method = c("subplex","L-BFGS-B"),
niter = 100, FAIL = 1e+200, hessian = FALSE, CI = 0.95), ncores=NULL, ...)fitContinuous returns a list with the following four elements:
is the function used to compute the model likelihood. The returned function (lik) takes arguments that are necessary for the given model.
For instance, if estimating a Brownian-motion model with unknown standard error, the arguments (pars) to the lik function would be sigsq and SE. By default, the function evaluates the likelihood of the model by assuming the maximum likelihood root state. This behavior can be changed in the call to lik with lik(pars, root=ROOT.GIVEN) where pars includes a value for the root state (z0). See Examples for a demonstration. The tree and data are stored internally within the lik function, which permits those elements to be efficiently reused when computing the likelihood under different parameter values
is a matrix of the used bounds for the relevant parameters estimated in the model. Warnings will be issued if any parameter estimates occur at the supplied (or default) parameter bounds
is a matrix of results from optimization. Rownames of the res matrix are the optimization methods
(see optim and subplex). The columns in the res matrix are the estimated
parameter values, the estimated model likelihood, and an indication of optimization convergence. Values of convergence not
equal to zero are not to be trusted
is a list of the primary results: estimates of the parameters, the maximum-likelihood estimate (lnL) of the model, the
optimization method used to compute the MLE, the number of model parameters (k, including one parameter for the root state), the AIC (aic),
sample-size corrected AIC (aicc). The number of observations for AIC computation is taken to be the number of trait values observed.
If the Hessian is used, confidence intervals on the parameter estimates (CI) and the Hessian matrix (hessian) are also returned
a phylogenetic tree of class phylo
data vector for a single trait, with names matching tips in phy
a single value or named vector of standard errors associated with values in dat; if any elements in the vector SE are NA, SE will be estimated
model to fit to comparative data (see Details)
range to constrain parameter estimates (see Details)
settings used for optimization of the model likelihood
Number of cores. If NULL then number of cores is detected
additional arguments to be passed to the internal likelihood function bm.lik
LJ Harmon, W Challenger, and JM Eastman
This function fits various likelihood models for continuous character evolution. The function returns parameter estimates and the likelihood for univariate datasets.
The model likelihood is maximized using methods available in optim as well as subplex.
Optimization methods to be used within optim can be specified through the control object (i.e., control$method).
A number of random starting points are used in optimization and are given through the niter element within the control object
(e.g., control$niter). The FAIL value within the control object should be a large value
that will be considerably far from -lnL of the maximum model likelihood. In most cases, the default setting
for control$FAIL will be appropriate. The Hessian may be used to compute confidence intervals
(CI) for the parameter estimates if the hessian element in control is TRUE.
Beware: difficulty in finding the optimal solution is determined by an interaction between the nature and complexity of the
likelihood space (which is data- and model-dependent) and the numerical optimizer used to explore the space. There is never a
guarantee that the optimal solution is found, but using many random starting points (control$niter) and many optimization methods
(control$method) will increase these odds.
Bounds for the relevant parameters of the fitted model may be given through the bounds argument. Bounds
may be necessary (particularly under the OU model) if the likelihood surface is characterized by a long,
flat ridge which can be exceedingly difficult for optimization methods. Several bounds can be given at a time
(e.g., bounds=list(SE=c(0,0.1),alpha=c(0,1)) would constrain measurement error as well as the 'constraint'
parameter of the Ornstein-Uhlenbeck model). Default bounds under the different models are given below.
Possible models are as follows:
BM is the Brownian motion model (Felsenstein 1973), which assumes the correlation structure among trait values is proportional to the extent of shared ancestry
for pairs of species. Default bounds on the rate parameter are sigsq=c(min=exp(-500),max=exp(100)). The same bounds are applied to all other
models, which also estimate sigsq
OU is the Ornstein-Uhlenbeck model (Butler and King 2004), which fits a random walk with a central tendency with an attraction strength proportional to the parameter alpha.
The OU model is called the hansen model in ouch, although the way the parameters are fit is slightly different here. Default
bounds are alpha = c(min = exp(-500), max = exp(1))
EB is the Early-burst model (Harmon et al. 2010) and also called the ACDC model (accelerating-decelerating; Blomberg et al. 2003). Set by the a rate parameter, EB
fits a model where the rate of evolution increases or decreases exponentially through time, under the model r[t] = r[0] * exp(a * t), where r[0] is the
initial rate, a is the rate change parameter, and t is time. The maximum bound is set to -0.000001, representing a decelerating rate of evolution. The minimum bound is set to log(10^-5)/depth of the tree.
rate_trend is a diffusion model with linear trend in rates through time (toward larger or smaller rates). Used to be denominated the "trend" model, which is still accepted by fitContinuous for backward compatibility. Default bounds are slope = c(min = -100, max = 100)
lambda is one of the Pagel (1999) models that fits the extent to which the phylogeny predicts covariance among trait values for species. The model effectively transforms the tree:
values of lambda near 0 cause the phylogeny to become more star-like, and a lambda value of 1 recovers the BM model. Default
bounds are lambda = c(min = exp(-500), max = 1
kappa is a punctuational (speciational) model of trait evolution (Pagel 1999), where character divergence is related to the number of speciation events between two species. Note that if
there are speciation events that are missing from the given phylogeny (due to extinction or incomplete sampling), interpretation under the kappa model may be difficult. Considered as a tree
transformation, the model raises all branch lengths to an estimated power (kappa). Default bounds are kappa = c(min = exp(-500), max = 1)
delta is a time-dependent model of trait evolution (Pagel 1999). The delta model is similar to ACDC insofar as the delta model fits the relative contributions of
early versus late evolution in the tree to the covariance of species trait values. Where delta is greater than 1, recent evolution has been relatively fast; if delta is less
than 1, recent evolution has been comparatively slow. Intrepreted as a tree transformation, the model raises all node depths to an estimated power (delta). Default bounds are delta = c(min = exp(-500), max = 3)
mean_trend is a model of trait evolution with a directional drift or trend component (i.e., toward smaller or larger values through time). This model is sensible only for non-ultrametric trees, as the likelihood surface is entirely flat with respect to the slope of the trend if the tree is ultrametric. The model used to be denominated the "drift" model, which is still accepted by fitContinuous for backward compatibility. Default bounds are drift = c(min = -100, max = 100)
white is a white-noise (non-phylogenetic) model, which assumes data come from a single normal distribution with no covariance structure among species. The variance parameter sigsq
takes the same bounds defined under the BM model
Blomberg SP, T Garland, and AR Ives. 2003. Testing for phylogenetic signal in comparative data: behavioral traits are more labile. Evolution 57:717-745.
Butler MA and AA King, 2004. Phylogenetic comparative analysis: a modeling approach for adaptive evolution. American Naturalist 164:683-695.
Felsenstein J. 1973. Maximum likelihood estimation of evolutionary trees from continuous characters. American Journal of Human Genetics 25:471-492.
Harmon LJ et al. 2010. Early bursts of body size and shape evolution are rare in comparative data. Evolution 64:2385-2396.
Pagel M. 1999. Inferring the historical patterns of biological evolution. Nature 401:877-884