# Multivariate continuous trait evolution for a stationary time series (Ornstein-Uhlenbeck model) This function allows the fitting of a multivariate Ornstein-Uhlenbeck (OU) model to a time series. Species measurement errors or dispersions can also be included in the model. ```r mvOUTS(times, data, error = NULL, param = list(sigma = NULL, alpha = NULL, vcv = "randomRoot", decomp = c("cholesky","spherical","eigen","qr", "diagonal","upper","lower")), method = c("rpf", "inverse", "pseudoinverse", "univarpf"), scale.height = FALSE, optimization = c("L-BFGS-B", "Nelder-Mead", "subplex"), control = list(maxit = 20000), precalc = NULL, diagnostic = TRUE, echo = TRUE) ``` ## Arguments - `times`: Time series - vector of sample ages. - `data`: Matrix or data frame with species in rows and continuous traits in columns. NA values are allowed. - `error`: Matrix or data frame with species in rows and continuous trait sampling variance (squared standard errors) in columns. - `param`: List of arguments to be passed to the function. See details below. - `method`: Choose between "rpf", "inverse", "pseudoinverse", or "univarpf" for computing the log-likelihood during the fitting process. See details below. - `scale.height`: Whether the time series should be scaled to unit length or not. - `optimization`: Methods used by the optimization routines (see ?optim and ?subplex for details). The "fixed" method returns the log-likelihood function only. - `control`: Max. bound for the number of iteration of the optimizer; other options can be fixed in the list (see ?optim or ?subplex). - `precalc`: Optional. Precalculation of fixed parameters. See ?mvmorph.Precalc for details. - `diagnostic`: Whether the convergence diagnostics should be returned or not. - `echo`: Whether the results must be returned or not. ## Details The mvOUTS function fits a multivariate model of trait evolution on a time series according to an Ornstein-Uhlenbeck process. The user can include measurement errors to the analyzed dataset. The "method" argument allows the user to try different algorithms for computing the log-likelihood. The **"rpf"**, **"univarpf"** (for univariate analysis) methods use fast GLS algorithms based on factorization for avoiding the computation of the inverse of the variance-covariance matrix and its determinant for the log-likelihood estimation. The **"inverse"** approach uses the **"stable"** standard explicit computation of the inverse and determinant of the matrix and is therefore slower. The **"pseudoinverse"** method uses a generalized inverse that is safer for matrix near singularity but highly time consuming. See ?mvLL for details. Arguments in the "param" **list** are: "sigma" or "alpha" - Starting values for the likelihood search can be specified through the "alpha" and "sigma" arguments in the param list. It is also possible to test for the significance of the off-diagonal sigma (scatter) and alpha (drift) matrix in the full model by making comparison with a constrained model (using sigma="constraint", or alpha="constraint") in the "param" argument list. 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))). "decomp" - You can further constrain the alpha matrix by specifying the decomposition of the matrix through the "decomp" argument in the "param" list. Indeed, the multivariate Ornstein-Uhlenbeck model is described by the spectral decomposition of the alpha matrix. Thus it is possible to parameterize the alpha matrix to be decomposable using various parameterizations (e.g., on its eigenvalues with different biological interpretations; Sy et al. 1997, Bartoszek et al. 2012). For a symmetric matrix parameterization the user can choose the **"cholesky"**, **"eigen"**, or **"spherical"** option. For general square (non-symmetric) matrices the **"svd"**, **"qr"** and **"schur"** parameterizations can be used. The **"schur"** parameterization constrains the eigenvalues of the alpha matrix to be real numbers. The **"svd+"**, **"qr+"** or **"eigen+"** options forces the eigenvalues to be positives by taking their logarithm. It is also possible to specify **"diagonal"** which is similar to the use of the "constraint" argument for the "alpha" argument, or to use **"equal"** and **"equaldiagonal"**. Finally, one can specify that the alpha matrix is **"upper"** or **"lower"** triangular (i.e., one process affect the other unilateraly). Details can be found in the package vignette: browseVignettes("mvMORPH"). "decompSigma" - The sigma matrix is parameterized by various methods to ensure its positive definiteness (Pinheiro and Bates, 1996). These methods can be accessed through the "decompSigma" argument and are the **"cholesky"**, **"eigen+"**, and **"spherical"** parameterization. The sigma matrix can also be forced to be diagonal using **"diagonal"** or **"equaldiagonal"** and forced to have the same variances using **"equal"**. Details can be found in the package vignette: browseVignettes("mvMORPH"). "vcv" - It is possible to specify in the "param" list what kind of variance-covariance matrix to use with the "vcv" argument, depending on how the root is treated. The **vcv="randomRoot"** option assumes that the value at the root is a random variable with the stationary distribution of the process. The **vcv="fixedRoot"** option assumes that the root is a fixed parameter. "root" - If root=TRUE, the ancestral state and the optimum (stationary mean) are estimated, otherwise (root=FALSE) the ancestral (initial) state and the optimum (long-term expectation) are assumed to be the same. Note: for the "decomp" and "decompSigma arguments, an user-defined matrix with integer values taken as indices of the parameters to be estimated can be provided. See ?mvBM and ?mvRWTS. ## Returns - **LogLik**: The log-likelihood of the optimal model. - **AIC**: Akaike Information Criterion for the optimal model. - **AICc**: Sample size-corrected AIC. - **theta**: Estimated ancestral states. - **alpha**: Matrix of estimated alpha values (strength of selection, drift matrix). - **sigma**: Evolutionary rate matrix (scatter). - **convergence**: Convergence status of the optimizing function; "0" indicates convergence. (See ?optim for details). - **hess.values**: Reliability of the likelihood estimates calculated through the eigen-decomposition of the hessian matrix. "0" means that a reliable estimate has been reached. See details on ?mvOU. - **param**: List of model fit parameters (optimization, method, model, number of parameters...). - **llik**: The log-likelihood function evaluated in the model fit "$llik(par, root.mle=TRUE)". ## 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. Clavel J., Escarguel G., Merceron G. 2015. mvMORPH: an R package for fitting multivariate evolutionary models to morphometric data. Methods Ecol. Evol. 6(11):1311-1319. Hunt G., Bell M.A., Travis M.P. 2008. Evolution toward a new adaptive optimum: phenotypic evolution in a fossil stickleback lineage. Evolution 62(3):700-710. Pinheiro J.C., Bates D.M. 1996. Unconstrained parameterizations for variance-covariance matrices. Stat. Comput. 6:289-296. Sy J.P., Taylor J.M.G., Cumberland W.G. 1997. A stochastic model for the analysis of bivariate longitudinal AIDS data. Biometrics. 53:542-555. ## Author(s) Julien Clavel ## See Also `mvMORPH` `halflife` `stationary` `mvOU` `mvRWTS` `mvBM` `mvEB` `mvSHIFT` `mvSIM` `LRT` `optim` ## Examples ```r # Simulate the time series set.seed(14) timeseries <- 0:49 # Parameters with general alpha matrix on two competitive species (or two traits) # asymetric (drift) matrix with intervention from the lowest layer alpha <- matrix(c(0.15,0,0.1,0.1),2,2) # scatter matrix sigma <- matrix(c(0.01,0.005,0.005,0.01),2) # ancestral states and long term optimum expectation theta <- matrix(c(0,1,0,.5),2) # columns=traits # Simulate the data traits <- mvSIM(timeseries, model="OUTS", param=list(theta=theta, alpha=alpha, sigma=sigma)) # Plot the time series matplot(traits,type="o",pch=1, xlab="Time (relative)") fit1 <- mvOUTS(timeseries, traits, param=list(decomp="qr")) fit2 <- mvOUTS(timeseries, traits, param=list(decomp="eigen")) fit3 <- mvOUTS(timeseries, traits, param=list(decomp="diagonal")) results <- list(fit1,fit2,fit3) aicw(results) # Simulate under the MLE traits2 <- simulate(fit1,tree=timeseries) matplot(traits2, type="o", pch=1, xlab="Time (relative)") mvOUTS(timeseries, traits2, param=list(decomp="eigen")) mvOUTS(timeseries, traits2, param=list(decomp="diagonal")) mvOUTS(timeseries, traits2, param=list(decomp="upper")) mvOUTS(timeseries, traits2, param=list(decomp="lower")) # try user defined constraints set.seed(100) ts <- 49 timeseries <- 1:ts sigma <- matrix(c(0.01,0.005,0.003,0.005,0.01,0.003,0.003,0.003,0.01),3) # upper triangular matrix with effect of trait 2 on trait 1. alpha <- matrix(c(0.4,0,0,-0.5,0.3,0,0,0,0.2),3,3) theta <- matrix(c(0,0,0,1,0.5,0.5),byrow=TRUE, ncol=3); root=TRUE data <- mvSIM(timeseries, model="OUTS", param=list(alpha=alpha, sigma=sigma, theta=theta, root=root, names_traits=c("sp 1", "sp 2", "sp 3"))) # plot matplot(data, type="o", pch=1, xlab="Time (relative)") legend("bottomright", inset=.05, legend=colnames(data), pch=19, col=c(1,2,3), horiz=TRUE) # define an user constrained drift matrix indice <- matrix(NA,3,3) diag(indice) <- c(1,2,3) indice[1,2] <- 4 # fit the model fit_1 <- mvOUTS(timeseries, data, param=list(vcv="fixedRoot", decomp=indice)) fit_2 <- mvOUTS(timeseries, data, param=list(vcv="fixedRoot", decomp="diagonal")) LRT(fit_1, fit_2) ```