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

```
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, ...)
```

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`

.

the named vector of coefficients.

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

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

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

the maximum of the log likelihood.

the number of all parameters of the model.

AIC value of the model.

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

fitted values

raw residuals

response

design matrix

number of observations (tips in the tree)

number of dependent variables

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

the model formula

the original call to the function

the phylogenetic model for the covariance

(`boot > 0`

only) bootstrap means of the parameters estimated.

(`boot > 0`

only) bootstrap standard deviations of the estimated parameters.

(`boot > 0`

only) bootstrap 95% confidence interval.

(`boot > 0`

only) bootstrap mean and standard deviation of the logs of the estimated covariance parameters.

(`boot > 0`

only) number of independent bootstrap replicates for which
`phylolm`

failed.

(`boot > 0`

and `full.matrix`

= `TRUE`

only) matrix of all bootstrap estimates.

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.

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.

# NOT RUN { 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) # }