| fit.multivariate.OU.user.defined {evoTS} | R Documentation |
Fit user-defined multivariate Ornstein-Uhlenbeck models to multivariate evolutionary sequence (time-series) data.
Description
Function to find maximum likelihood solutions to a multivariate Ornstein-Uhlenbeck model fitted using user-defined A and R matrices.
Usage
fit.multivariate.OU.user.defined(
yy,
A.user = NULL,
R.user = NULL,
method = "Nelder-Mead",
hess = FALSE,
pool = TRUE,
trace = FALSE,
iterations = NULL,
iter.sd = NULL,
user.init.diag.A = NULL,
user.init.upper.diag.A = NULL,
user.init.lower.diag.A = NULL,
user.init.diag.R = NULL,
user.init.off.diag.R = NULL,
user.init.theta = NULL,
user.init.anc = NULL
)
Arguments
yy |
a multivariate evoTS object. |
A.user |
the pull matrix. A user-defined A matrix. |
R.user |
the drift matrix. A user-defined R matrix. |
method |
optimization method, passed to function optim. Default is "Nelder-Mead". |
hess |
logical, indicating whether to calculate standard errors from the Hessian matrix. |
pool |
indicating whether to pool variances across samples |
trace |
logical, indicating whether information on the progress of the optimization is printed. |
iterations |
the number of times the optimization method is run from different starting points. Default is NULL, meaning the optimization is run once. |
iter.sd |
defines the standard deviation of the Gaussian distribution from which starting values for the optimization routine is run. Default is 1. |
user.init.diag.A |
starting values for the optimization routine of the diagonal elements of the A matrix. Default is NULL. |
user.init.upper.diag.A |
starting values for the optimization routine of the upper diagonal elements of the A matrix. Default is NULL. |
user.init.lower.diag.A |
starting values for the optimization routine of the lower diagonal elements of the A matrix. Default is NULL. |
user.init.diag.R |
starting values for the optimization routine of the diagonal elements of the R matrix. Default is NULL. |
user.init.off.diag.R |
starting values for the optimization routine of the off-diagonal elements of the R matrix. Default is NULL. |
user.init.theta |
starting values for the optimization routine of the optima. Default is NULL. |
user.init.anc |
starting values for the optimization routine of the ancestral values. Default is NULL. |
Details
This function provides users the flexibility to define their own A and R matrices. The possibility to define any A matrices enable detailed investigation of specific evolutionary hypotheses. The parameters to be estimated in the matrices are indicated by the value 1. All other entries in the matrix must be 0.
The function searches - using an optimization routine - for the maximum-likelihood solution for the chosen multivariate Ornstein-Uhlenbeck model. The argument 'method' is passed to the 'optim' function and is included for the convenience of users to better control the optimization routine. Note that the the default method (Nelder-Mead) seems to work for most evolutionary sequences. The method L-BFGS-B allows box-constraints on some parameters (e.g. non-negative variance parameters) and is faster than Nelder-Mead, but is less stable than the default method (Nelder-Mead).
Initial estimates to start the optimization come from maximum-likelihood estimates of the univariate Ornstein-Uhlenbeck model (from the paleoTS package) fitted to each time-series separately.
It is good practice to repeat any numerical optimization procedure from different starting points. This is especially important for complex models as the log-likelihood surface might contain more than one peak. The number of iterations is controlled by the argument 'iterations'. The function will report the model parameters from the iteration with the highest log-likelihood.
There is no guarantee that the likelihood can be computed with the initial parameters provided by the function. The starting values for fitting the multivariate OU model are based on maximum likelihood parameter estimates for the univariate OU model fitted to each trait separately, which seems to provide sensible (and working) initial parameter estimates for almost all tested data sets. However, the provided initial parameters may fail depending on the nature of the data. If an error message is returned saying "function cannot be evaluated at initial parameters", the user can try to start the optimization procedure from other initial parameter values using "user.init.diag.A", "user.init.upper.diag.A", "user.init.lower.diag.A", "user.init.diag.R", "user.init.off.diag.R", "user.init.theta", and "user.init.anc." It is usually the initial guess of the off-diagonal elements of the A and R matrices that prevents the optimization routine to work. It is therefore recommended to only try to change these initial values before experimenting with different starting values for the diagonal of the A and R matrices.
Value
First part of the output reports the log-likelihood of the model and its AICc score. The second part of the output is the maximum log-likelihood model parameters (ancestral.values, optima, A, and R). The half-life is also provided, which is the The last part of the output gives information about the number of parameters in the model (K), number of samples in the data (n) and number of times the optimization routine was run (iter).
Note
The models have been implemented to be compatible with the joint parameterization routine in the package paleoTS. The optimization is therefore fit using the actual sample values, with the autocorrelation among samples accounted for in the log-likelihood function. The joint distribution of sample means is multivariate normal, with means and variance-covariances determined by evolutionary parameters and sampling errors.
Author(s)
Kjetil Lysne Voje
References
Reitan, T., Schweder, T. & Henderiks, J. Phenotypic evolution studied by layered stochastic differential equations. Ann Appl Statistics 6, 1531–1551 (2012).
Bartoszek, K., Pienaar, J., Mostad, P., Andersson, S. & Hansen, T. F. A phylogenetic comparative method for studying multivariate adaptation. J Theor Biol 314, 204–215 (2012).
Clavel, J., Escarguel, G. & Merceron, G. mvmorph: an r package for fitting multivariate evolutionary models to morphometric data. Methods Ecol Evol 6, 1311–1319 (2015).
Examples
## Generate a evoTS object by simulating a multivariate dataset
x <- sim.multi.OU(15)
## Define an A matrix that is lower diagonal.
A <- matrix(c(1,0,1,1), nrow=2, byrow=TRUE)
## Define a diagonal R matrix.
R <- matrix(c(1,0,0,1), nrow=2, byrow=TRUE)
## Fit the multivariate Ornstein-Uhlenbeck model to the data. This example will run for a long time.
fit.multivariate.OU.user.defined(x, A.user=A, R.user=R, trace=TRUE)