compar.ou: Ornstein-Uhlenbeck Model for Continuous Characters

Description Usage Arguments Details Value Note Author(s) References See Also Examples

View source: R/compar.ou.R

Description

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.

Usage

1
compar.ou(x, phy, node = NULL, alpha = NULL)

Arguments

x

a numeric vector giving the values of a continuous character.

phy

an object of class "phylo".

node

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 phy has node labels. By default there is no change (same optimum thoughout lineages).

alpha

the value of alpha to be used when fitting the model. By default, this parameter is estimated (see details).

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 node.

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 alpha, and selecting the one that minimizes the deviance (see Hansen 1997 for an example).

If 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 phy.

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 labels of phy, and a warning message is issued.

Value

an object of class "compar.ou" which is list with the following components:

deviance

the deviance (= -2 * loglik).

para

a data frame with the maximum likelihood estimates and their standard-errors.

call

the function call.

Note

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 solve.

Author(s)

Emmanuel Paradis

References

Hansen, T. F. (1997) Stabilizing selection and the comparative analysis of adaptation. Evolution, 51, 1341–1351.

See Also

ace, compar.lynch, corBrownian, corMartins, pic

Examples

 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.

Example output

$deviance
[1] 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")
[1] "compar.ou"
$deviance
[1] 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")
[1] "compar.ou"
$deviance
[1] 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")
[1] "compar.ou"
$deviance
[1] 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")
[1] "compar.ou"

ape documentation built on April 25, 2021, 9:06 a.m.