This function fits an Ornstein–Uhlenbeck model giving a phylogenetic tree, and a continuous character. The user specifies the node(s) where the optimum changes. The parameters are estimated by maximum likelihood; their standard-errors are computed assuming normality of these estimates.
a numeric vector giving the values of a continuous character.
an object of class
a vector giving the number(s) of the node(s) where the
parameter ‘theta’ (the trait optimum) is assumed to change. The
node(s) can be specified with their labels if
the value of alpha to be used when fitting the model. By default, this parameter is estimated (see details).
The Ornstein–Uhlenbeck (OU) process can be seen as a generalization of the Brownian motion process. In the latter, characters are assumed to evolve randomly under a random walk, that is change is equally likely in any direction. In the OU model, change is more likely towards the direction of an optimum (denoted theta) with a strength controlled by a parameter denoted alpha.
The present function fits a model where the optimum parameter
theta, is allowed to vary throughout the tree. This is
specified with the argument
node: theta changes
after each node whose number is given there. Note that the optimum
changes only for the lineages which are descendants of this
Hansen (1997) recommends to not estimate alpha together
with the other parameters. The present function allows this by giving
a numeric value to the argument
alpha. By default, this
parameter is estimated, but this seems to yield very large
standard-errors, thus validating Hansen's recommendation. In practice,
a “poor man estimation” of alpha can be done by
repeating the function call with different values of
selecting the one that minimizes the deviance (see Hansen 1997 for an
x has names, its values are matched to the tip labels of
phy, otherwise its values are taken to be in the same order
than the tip labels of
The user must be careful here since the function requires that both
series of names perfectly match, so this operation may fail if there
is a typing or syntax error. If both series of names do not match, the
values in the
x are taken to be in the same order than the tip
phy, and a warning message is issued.
an object of class
"compar.ou" which is list with the following
the deviance (= -2 * loglik).
a data frame with the maximum likelihood estimates and their standard-errors.
the function call.
The inversion of the variance-covariance matrix in the likelihood
function appeared as somehow problematic. The present implementation
uses a Cholevski decomposition with the function
chol2inv instead of the usual function
Hansen, T. F. (1997) Stabilizing selection and the comparative analysis of adaptation. Evolution, 51, 1341–1351.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
data(bird.orders) ### This is likely to give you estimates close to 0, 1, and 0 ### for alpha, sigma^2, and theta, respectively: compar.ou(x <- rnorm(23), bird.orders) ### Much better with a fixed alpha: compar.ou(x, bird.orders, alpha = 0.1) ### Let us 'mimick' the effect of different optima ### for the two clades of birds... x <- c(rnorm(5, 0), rnorm(18, 5)) ### ... the model with two optima: compar.ou(x, bird.orders, node = 25, alpha = .1) ### ... and the model with a single optimum: compar.ou(x, bird.orders, node = NULL, alpha = .1) ### => Compare both models with the difference in deviances ## which follows a chi^2 with df = 1.
$deviance  61.58694 $para estimate stderr alpha 0.42272352 86.7549317 sigma2 0.85198597 0.1776928 theta1 -0.05165616 0.1360976 $call compar.ou(x = x <- rnorm(23), phy = bird.orders) attr(,"class")  "compar.ou" $deviance  61.6933 $para estimate stderr sigma2 0.8592617 0.1792100 theta1 -0.0601906 0.1494551 $call compar.ou(x = x, phy = bird.orders, alpha = 0.1) attr(,"class")  "compar.ou" $deviance  69.2047 $para estimate stderr sigma2 1.1911272 0.2484164 theta1 0.6231958 0.3756685 theta2 5.0216417 0.2000396 $call compar.ou(x = x, phy = bird.orders, node = 25, alpha = 0.1) attr(,"class")  "compar.ou" $deviance  96.70752 $para estimate stderr sigma2 3.938013 0.8212964 theta1 4.041715 0.3199533 $call compar.ou(x = x, phy = bird.orders, node = NULL, alpha = 0.1) attr(,"class")  "compar.ou"
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.