mvMORPH (version 1.0.2)

mvOU: Multivariate Ornstein-Uhlenbeck model of continuous traits evolution

Description

This function allows the fitting of a multivariate OU model while allowing a given tree branch to be subdivided by multiple selective regimes using SIMMAP-like mapping of ancestral states. Species measurement error or dispersion can also be included in the model.

Usage

mvOU(tree, data, error = NULL, sigma = NULL, alpha = NULL, 
 model = c("OUM", "OU1"), scale.height=FALSE, diagnostic = TRUE, 
 method = c("L-BFGS-B","Nelder-Mead","subplex"), pseudoinverse=FALSE, 
 echo = TRUE, control=list(maxit=20000))

Arguments

tree

Phylogenetic tree with mapped ancestral states in SIMMAP format. (See make.simmap function from phytools package). A "phylo" object can be used with model "OU1".

data

Matrix or data frame with species in rows and continuous traits in columns.

error

Matrix or data frame with species in rows and continuous traits standard error (squared) in columns.

sigma

Starting values of the sigma matrix before optimization (optional). Note that if you use "constraint" instead of NULL, sigma covariances matrix off-diagonal are omitted. You can also use a list argument to provide starting values for each traits (see deta

alpha

Starting values of the alpha matrix before optimization (optional). Note that if you use "constraint" instead of NULL, alpha covariances matrix off-diagonal are omitted. You can also use a list argument to provide starting values for each traits (see deta

model

Choose between "OUM" for a multiple selective regime model, or "OU1" for a unique selective regime for the whole tree.

scale.height

Whether the tree should be scaled to length 1.

diagnostic

Whether the diagnostics of convergence should be returned.

method

Methods used by the optimization function. (see ?optim and ?subplex for details).

pseudoinverse

Whether Moore-Penrose pseudoinverse should be used in calculation (slower)

echo

Should the results be returned

control

Max. bound for the number of iteration of the optimizer; other options can be fixed on the list (see ?optim or ?subplex).

Value

• LogLik
  Log-Likelihood of the optimized model.
• AIC
  Akaike Information Criterion calculated for the best model.
• AICc
  Akaike Information Criterion corrected for small sample size.
• theta
  Matrix of estimated theta values for each trait and selective regimes.
• alpha.mat
  Matrix of estimated alpha values (strength of selection) for the studied traits (diagonal).
• sigma.mat
  Matrix of estimated sigma values (drift) for the studied traits (diagonal).
• alpha
  alpha values estimated by the optimizing function.
• sigma
  sigma values estimated by the optimizing function.
• alpha.se
  standard errors of the alpha estimated values.
• sigma.se
  standard errors of the sigma estimated values.
• convergence
  Convergence status of the optimizing function. See ?optim for more details.
• hessian
  Hessian matrix (see details).
• hess.value
  If value is 0, this mean that the eigen-values of the hessian matrix are all positives. Value of 1 means that the optimizing function may have converged to a non-reliable estimation.

Details

The mvOU function fit a multivariate model of evolution according to a Ornstein-Uhlenbeck process. This function use a modified version of the OUCH package (Butler & King, 2004) and extends it by allowing the user to incorporate measurement error and use SIMMAP-like mapping of ancestral states. SIMMAP mapping allows one to assign parts of branchs to different selective regimes, and allows testing for changes in trait variance which are not synchronous with the species divergence event. Mapping of ancestral states can be done using the "make.simmap", "make.era.map" or "paintSubTree" functions from the "phytools" package. It is possible to test for the significance of the off-diagonal sigma and alpha matrix in the full model by making comparison with a constrained model (using sigma="constraint", or alpha="constraint"). You can also provide starting values for the constrained model. For instance, for two traits use sigma=list("constraint", c(0.5,0.5)) (or alpha=list("constraint", c(0.5,0.5))). The hessian is the matrix of second order partial derivatives of the likelihood function with respect to the maximum likelihood parameter values. This matrix provide a measure of the steepness of the likelihood surface in the region of the optimum. The eigen-decomposition of the hessian matrix returned by the optimizing function allow to assess the reliability of the fit of the model (even if the optimizer have converged). When the optimization function does not converge on a stable result, the user may consider increasing the "maxit" argument in the "control" option, or try a simpler model with less parameters to estimate. Changing the starting values ("alpha" and "sigma" options) as well as the optimizing function ("method" options) may help sometimes (e.g., alpha=runif(3) for two traits analysis with random starting values (i.e., the lower triangular alpha matrix)). Note that Bartoszek et al.,(2012) proposed the mvSLOUCH package dedicated to multivariate Ornstein-Uhlenbeck processes, which allows fitting models with a randomly evolving predictor variables. They also provide a detailled mathematical description of the multivariate OU model.

References

Bartoszek K., Pienaar J., Mostad P., Andersson S., Hansen T.F. 2012. A phylogenetic comparative method for studying multivariate adaptation. J. Theor. Biol. 314:204-215. Beaulieu J.M., Jhwueng D.-C., Boettiger C., O'Meara B.C. 2012. Modeling stabilizing selection: Expanding the Ornstein-Uhlenbeck model of adaptive evolution. Evolution. 66:2369-2389. Butler M.A., King A.A. 2004. Phylogenetic comparative analysis: a modeling approach for adaptive evolution. Am. Nat. 164:683-695.

See Also

mvMORPH mvBM mvEB mvSHIFT optim make.simmap make.era.map paintSubTree

Examples

## Toy Exemple
  set.seed(123)
  # Generating a random tree
  tree<-pbtree(n=25)

  # Setting the regime states of tip species
  sta<-as.vector(c(rep("Forest",10),rep("Savannah",15))); names(sta)<-tree$tip.label

  # Making the simmap tree with mapped states
  tree<-make.simmap(tree,sta , model="ER", nsim=1)
  col<-c("blue","orange"); names(col)<-c("Forest","Savannah")

  # Plot of the phylogeny for illustration
  plotSimmap(tree,col,fsize=0.6,node.numbers=FALSE,lwd=3, pts=FALSE)

  # 2 Random traits evolving along the phylogeny
  data<-data.frame(head.size=rTraitCont(tree), mouth.size=rTraitCont(tree))

  # Names of the species
  rownames(data)<-tree$tip.label

## Run the analysis!
  # OU model with unique optima
  mvOU(tree,data, model="OU1", scale.height=TRUE)
  
  # Analysis with multiple optima!
  mvOU(tree,data, scale.height=TRUE)