ouch (version 2.20)

hansen: Ornstein-Uhlenbeck models of trait evolution

Description

The function hansen fits an Ornstein-Uhlenbeck model to data. The fitting is done using optim or subplex.

Usage

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

Value

hansen returns an object of class hansentree.

Arguments

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 nchar named vectors, or a data frame containing nchar data variables. There must be an entry per variable for every node in the tree; use NA to represent missing data. If the data are supplied as one or more named vectors, the names attributes are taken to correspond to the node names specified when the ouchtree was constructed (see ouchtree). If the data are supplied as a data-frame, the rownames serve that purpose.

tree

A phylogenetic tree, specified as an ouchtree object.

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 nchar quantitative characters, then one can specify a single set of regimes for all characters or a list of nchar regime specifications, one for each character.

sqrt.alpha, sigma

These are used to initialize the optimization algorithm. The selection strength matrix \(\alpha\) and the random drift variance-covariance matrix \(\sigma^2\) are parameterized by their matrix square roots. Specifically, these initial guesses are each packed into lower-triangular matrices (column by column). The product of this matrix with its transpose is the \(\alpha\) or \(\sigma^2\) matrix. See Details for more information.

fit

If fit=TRUE, then the likelihood will be maximized. If fit=FALSE, the likelihood will be evaluated at the specified values of sqrt.alpha and sigma; the optima theta will be returned as well.

method

The method to be used by the optimization algorithm. See subplex::subplex and stats::optim for information on the available options.

hessian

If hessian=TRUE, then the Hessian matrix will be computed by optim.

...

Additional arguments will be passed as control options to optim or subplex. See stats::optim() and subplex::subplex() for information on the available options.

Author

Aaron A. King

Details

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\sigma^T\) 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.

References

T.F. Hansen. 1997. Stabilizing selection and the comparative analysis of adaptation. Evolution, 51:1341--1351. tools:::Rd_expr_doi("10.1111/j.1558-5646.1997.tb01457.x").

Butler, M.A. and A.A. King. 2004. Phylogenetic comparative analysis: a modeling approach for adaptive evolution. American Naturalist 164:683--695. tools:::Rd_expr_doi("10.1086/426002").

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. tools:::Rd_expr_doi("10.1093/sysbio/syv043").

See Also

stats::optim, subplex::subplex, bimac, anolis.ssd

Other phylogenetic comparative models: brown(), ouch-package, ouchtree, paint()

Examples

## Analysis of sexual size dimorphism data
 ## Save time for CRAN
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.
 ## Save time for CRAN
### The data were taken from package 'geiger' due to the latter being orphaned.
if (requireNamespace("ape")) {

  data(geospiza)
  plot(geospiza$phy)
  print(geospiza$dat)
  
### make an ouchtree out of the phy-format tree
  ot <- ape2ouch(geospiza$phy)

### merge data with tree info
  otd <- as(ot,"data.frame")
  otd <- merge(otd,geospiza$dat,by.x="labels",by.y="row.names",all=TRUE)
### row-names are used by 'hansen'
  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))
  plot(ot)

  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)
  plot(h1)

}