phylolm {phylolm}

Phylogenetic Linear Model

Description

Fits a phylogenetic linear regression model. The likelihood is calculated with an algorithm that is linear in the number of tips in the tree.

Usage

phylolm(formula, data = list(), phy, model = c("BM", "OUrandomRoot",
       "OUfixedRoot", "lambda", "kappa", "delta", "EB", "trend"),
       lower.bound = NULL, upper.bound = NULL,
       starting.value = NULL, measurement_error = FALSE,
       boot=0,full.matrix = TRUE, ...)

Arguments

formula

a model formula.

data

a data frame containing variables in the model. If not found in data, the variables are taken from current environment.

phy

a phylogenetic tree of type phylo with branch lengths.

model

a model for the covariance (see Details).

lower.bound

optional lower bound for the optimization of the phylogenetic model parameter.

upper.bound

optional upper bound for the optimization of the phylogenetic model parameter.

starting.value

optional starting value for the optimization of the phylogenetic model parameter.

measurement_error

a logical value indicating whether there is measurement error sigma2_error (see Details).

boot

number of independent bootstrap replicates, 0 means no bootstrap.

full.matrix

if TRUE, the full matrix of bootstrap estimates (coefficients and covariance parameters) will be returned.

...

further arguments to be passed to the function optim.

Details

This function uses an algorithm that is linear in the number of tips in the tree to calculate the likelihood. Possible phylogenetic models for the error term are the Brownian motion model (BM), the Ornstein-Uhlenbeck model with an ancestral state to be estimated at the root (OUfixedRoot), the Ornstein-Uhlenbeck model with the ancestral state at the root having the stationary distribution (OUrandomRoot), Pagel's \lambda model (lambda), Pagel's \kappa model (kappa), Pagel's \delta model (delta), the early burst model (EB), and the Brownian motion model with a trend (trend).

Using measurement error means that the covariance matrix is taken to be \sigma^2*V + \sigma^2_{error}*I where V is the phylogenetic covariance matrix from the chosen model, I is the identity matrix, and \sigma^2_{error} is the variance of the measurement error (which could include environmental variability, sampling error on the species mean, etc.).

By default, the bounds on the phylogenetic parameters are [10^{-7}/T,50/T] for \alpha, [10^{-7},1] for \lambda, [10^{-6},1] for \kappa, [10^{-5},3] for \delta and [-3/T,0] for rate, where T is the mean root-to-tip distance. [10^{-16}, 10^{16}] for the ratio sigma2_error/sigma2 (if measurement errors is used).

Bootstrapping can be parallelized using the future package on any future compatible back-end. For example, run library(future); plan(multiprocess)), after which bootstrapping will automatically occur in parallel. See plan for options.

Value

coefficients

the named vector of coefficients.

sigma2

the maximum likelihood estimate of the variance rate \sigma^2.

sigma2_error

the maximum likelihood estimate of the variance of the measurement errors.

optpar

the optimized value of the phylogenetic correlation parameter (alpha, lambda, kappa, delta or rate).

logLik

the maximum of the log likelihood.

p

the number of all parameters of the model.

aic

AIC value of the model.

vcov

covariance matrix for the regression coefficients, given the phylogenetic correlation parameter (if any).

fitted.values

fitted values

residuals

raw residuals

y

response

X

design matrix

n

number of observations (tips in the tree)

d

number of dependent variables

mean.tip.height

mean root-to-tip distance, which can help choose good starting values for the correlation parameter.

formula

the model formula

call

the original call to the function

model

the phylogenetic model for the covariance

bootmean

(boot > 0 only) bootstrap means of the parameters estimated.

bootsd

(boot > 0 only) bootstrap standard deviations of the estimated parameters.

bootconfint95

(boot > 0 only) bootstrap 95% confidence interval.

bootmeansdLog

(boot > 0 only) bootstrap mean and standard deviation of the logs of the estimated covariance parameters.

bootnumFailed

(boot > 0 only) number of independent bootstrap replicates for which phylolm failed.

bootstrap

(boot > 0 and full.matrix = TRUE only) matrix of all bootstrap estimates.

r.squared

The r^2 for the model.

adj.r.squared

The adjusted r^2 for the model.

Note

The tip labels in the tree are matched to the data names (row names in the data frame). If no data frame is provided through the argument data, taxon labels in the tree are matched to taxon labels in the response variable based on the row names of the response vector, and the taxa are assumed to come in the same order for all variables in the model.

For the delta model, the tree is rescaled back to its original height after each node's distance from the root is raised to the power delta. This is to provide a stable estimate of the variance parameter \sigma^2. For non-ultrametric trees, the tree is rescaled to maintain the longest distance from the root to its original value.

The trend model can only be used with non-ultrametric trees. For this model, one predictor variable is added to the model whose values are the distances from the root to every tip of the tree. The estimate of the coefficent for this variable forms the trend value.

Pagel's \lambda model and measurement error cannot be used together: the parameters \lambda, \sigma^2 and \sigma^2_{error} are not distinguishable (identifiable) from each other.

Author(s)

Lam Si Tung Ho

References

Ho, L. S. T. and Ane, C. 2014. "A linear-time algorithm for Gaussian and non-Gaussian trait evolution models". Systematic Biology 63(3):397-408.

Butler, M. A. and King, A. A. 2004. "Phylogenetic comparative analysis: A modeling approach for adaptive evolution". The American Naturalist 164:683-695.

Hansen, T. F. 1997. "Stabilizing selection and the comparative analysis of adaptation". Evolution 51:1341-1351.

Harmon, L. J. et al. 2010. "Early bursts of body size and shape evolution are rare in comparative data". Evolution 64:2385-2396.

Ho, L. S. T. and Ane, C. 2013. "Asymptotic theory with hierarchical autocorrelation: Ornstein-Uhlenbeck tree models". The Annals of Statistics 41(2):957-981.

Pagel, M. 1997. "Inferring evolutionary processes from phylogenies". Zoologica Scripta 26:331-348.

Pagel, M. 1999. "Inferring the historical patterns of biological evolution". Nature 401:877-884.

See Also

corBrownian, corMartins, corPagel, fitContinuous, pgls.

Examples

set.seed(123456)
tre = rcoal(60)
taxa = sort(tre$tip.label)
b0=0; b1=1;
x <- rTrait(n=1, phy=tre,model="BM",
            parameters=list(ancestral.state=0,sigma2=10))
y <- b0 + b1*x + 
     rTrait(n=1,phy=tre,model="lambda",parameters=list(
              ancestral.state=0,sigma2=1,lambda=0.5))
dat = data.frame(trait=y[taxa],pred=x[taxa])
fit = phylolm(trait~pred,data=dat,phy=tre,model="lambda")
summary(fit)

# adding measurement errors and bootstrap
z <- y + rnorm(60,0,1)
dat = data.frame(trait=z[taxa],pred=x[taxa])
fit = phylolm(trait~pred,data=dat,phy=tre,model="BM",measurement_error=TRUE,boot=100)
summary(fit)