hansen {ouch} | R Documentation |

The function `hansen`

fits an Ornstein-Uhlenbeck model to data.
The fitting is done using `optim`

or `subplex`

.

hansen( data, tree, regimes, sqrt.alpha, sigma, fit = TRUE, method = c("Nelder-Mead", "subplex", "BFGS", "L-BFGS-B"), hessian = FALSE, ... )

`data` |
Phenotypic data for extant species, i.e., species at the terminal twigs of the phylogenetic tree.
This can either be a single named numeric vector, a list of |

`tree` |
A phylogenetic tree, specified as an |

`regimes` |
A vector of codes, one for each node in the tree, specifying the selective regimes hypothesized to have been operative.
Corresponding to each node, enter the code of the regime hypothesized for the branch segment terminating in that node.
For the root node, because it has no branch segment terminating on it, the regime specification is irrelevant.
If there are |

`sqrt.alpha, sigma` |
These are used to initialize the optimization algorithm.
The selection strength matrix |

`fit` |
If |

`method` |
The method to be used by the optimization algorithm.
See |

`hessian` |
If |

`...` |
Additional arguments will be passed as |

The Hansen model for the evolution of a multivariate trait *X* along a lineage can be written as a stochastic differential equation (Ito diffusion)

*dX = alpha (theta(t)-X(t)) dt + sigma dB(t),*

where *t* is time along the lineage,
*theta(t)* is the optimum trait value, *B(t)* is a standard Wiener process (Brownian motion),
and *alpha* and *sigma* are matrices
quantifying, respectively, the strength of selection and random drift.
Without loss of generality, one can assume *sigma* is lower-triangular.
This is because only the infinitesimal variance-covariance matrix
*sigma^2 = sigma%*%transpose(sigma)*
is identifiable, and for any admissible variance-covariance matrix, we can choose *sigma* to be lower-triangular.
Moreover, if we view the basic model as describing evolution on a fitness landscape, then *alpha* will be symmetric.
If we further restrict ourselves to the case of stabilizing selection, *alpha* will be positive definite as well.
We make these assumptions and therefore can assume that the matrix *alpha* has a lower-triangular square root.

The `hansen`

code uses unconstrained numerical optimization to maximize the likelihood.
To do this, it parameterizes the *alpha* and *sigma^2* matrices in a special way:
each matrix is parameterized by `nchar*(nchar+1)/2`

parameters, where `nchar`

is the number of quantitative characters.
Specifically, the parameters initialized by the `sqrt.alpha`

argument of `hansen`

are used
to fill the nonzero entries of a lower-triangular matrix (in column-major order),
which is then multiplied by its transpose to give the selection-strength matrix.
The parameters specified in `sigma`

fill the nonzero entries in the lower triangular *sigma* matrix.
When `hansen`

is executed, the numerical optimizer maximizes the likelihood over these parameters.

`hansen`

returns an object of class `hansentree`

.

Aaron A. King

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

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

Cressler, C. E., Butler, M. A., and King, A. A. 2015. Detecting adaptive evolution in phylogenetic comparative analysis using the Ornstein-Uhlenbeck model. Systematic Biology, 64:953–968.

`stats::optim`

, `subplex::subplex`

, `bimac`

, `anolis.ssd`

Other phylogenetic comparative models:
`brown()`

,
`ouch-package`

,
`ouchtree`

,
`paint()`

## Analysis of sexual size dimorphism data tree <- with(anolis.ssd,ouchtree(node,ancestor,time/max(time),species)) plot(tree,node.names=TRUE) h1 <- brown(anolis.ssd['log.SSD'],tree) h1 plot(h1) h2 <- hansen(anolis.ssd['log.SSD'],tree,anolis.ssd['OU.1'],sqrt.alpha=1,sigma=1) h2 plot(h2) h3 <- hansen(anolis.ssd['log.SSD'],tree,anolis.ssd['OU.7'],sqrt.alpha=1,sigma=1) h3 plot(h3) ### Darwin's finches library(geiger) data(geospiza) ### check the correspondence between data and tree tips: print(nc <- with(geospiza,name.check(geospiza.tree,geospiza.data))) ### looks like one of the terminal twigs has no data associated ### drop that tip: tree <- with(geospiza,drop.tip(geospiza.tree,nc$tree_not_data)) dat <- as.data.frame(geospiza$dat) ### make an ouchtree out of the phy-format tree ot <- ape2ouch(tree) ### merge data with tree info otd <- as(ot,"data.frame") ### in these data, it so happens that the rownames correspond to node names ### we will exploit this correspondence in the 'merge' operation: dat$labels <- rownames(dat) otd <- merge(otd,dat,by="labels",all=TRUE) rownames(otd) <- otd$nodes print(otd) ### this data-frame now contains the data as well as the tree geometry ### now remake the ouch tree ot <- with(otd,ouchtree(nodes=nodes,ancestors=ancestors,times=times,labels=labels)) b1 <- brown(tree=ot,data=otd[c("tarsusL","beakD")]) summary(b1) ### evaluate an OU model with a single, global selective regime otd$regimes <- as.factor("global") h1 <- hansen( tree=ot, data=otd[c("tarsusL","beakD")], regimes=otd["regimes"], sqrt.alpha=c(1,0,1), sigma=c(1,0,1), maxit=10000 ) summary(h1)

[Package *ouch* version 2.17-1 Index]