Fit multivariate Unbiased Random Walk models to multivariate evolutionary sequence (time-series) data.

Function to find maximum likelihood solutions to a multivariate Unbiased Random Walk model.

Usage

fit.multivariate.URW(
  yy,
  R = "symmetric",
  r = "fixed",
  method = "L-BFGS-B",
  hess = FALSE,
  pool = TRUE,
  trace = FALSE,
  iterations = NULL,
  iter.sd = NULL
)

Arguments

yy

a multivariate evoTS object.

R

the drift matrix. The options are "diagonal" and "symmetric"

r

parameter describing the exponential increase/decrease in rate across time. The options are "fixed", "accel" and "decel".

method

optimization method, passed to function optim. Default is "L-BFGS-B".

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.

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, R). 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).

Details

The function allows the users to test six variants of multivariate Unbiased Random Walk models. There are two options for the structure of the R matrix. A "diagonal" R matrix means the stochastic changes in the traits are assumed to be uncorrelated. A "symmetric" R matrix means the stochastic changes in the traits are assumed to be correlated, i.e. that they are non-independent.

There are three options for the 'r' parameter. The "fixed" option means there is no change in the rate of change across time (r = 0). Setting r to "fixed" therefore fits a regular multivariate Unbiased Random Walk. The "decel" and "accel" options make the rate of change (the R matrix) decay (r < 0) and increase (r > 0) exponentially through time, respectively.

The function searches - using an optimization routine - for the maximum-likelihood solution for the chosen multivariate Unbiased Random Walk model. The argument 'method' is passed to the 'optim' function and is included for the convenience of users to better control the optimization routine. The the default method (L-BFGS-B) seems to work for most evolutionary sequences.

Initial estimates to start the optimization come from maximum-likelihood estimates of the univariate Unbiased Random Walk 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.

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.

References

Revell, L. J. & Harmon, L. Testing quantitative genetic hypotheses about the evolutionary rate matrix for continuous characters. Evolutionary Ecology Research 10, 311–331 (2008).

Voje, K. L. Testing eco‐evolutionary predictions using fossil data: Phyletic evolution following ecological opportunity. Evolution 74, 188–200 (2020).

Author

Kjetil Lysne Voje

Examples

## Generate an evoTS object by simulating a multivariate dataset
x <- sim.multi.URW(30)

## Fit a multivariate Unbiased Random Walk model to the data, allowing for correlated changes.
fit.multivariate.URW(x, R = "symmetric", r = "fixed")
#> $converge
#> [1] "Model converged successfully"
#> 
#> $modelName
#> [1] "Multivariate Random walk (R matrix with non-zero off-diagonal elements)"
#> 
#> $logL
#> [1] 22.2172
#> 
#> $AICc
#> [1] -31.93439
#> 
#> $ancestral.values
#> [1]  0.04037286 -0.24930581
#> 
#> $SE.anc
#> [1] NA
#> 
#> $R
#>           [,1]      [,2]
#> [1,] 0.7015753 0.2259568
#> [2,] 0.2259568 0.7795191
#> 
#> $SE.R
#> [1] NA
#> 
#> $method
#> [1] "Joint"
#> 
#> $K
#> [1] 5
#> 
#> $n
#> [1] 30
#> 
#> $iter
#> [1] NA
#> 
#> attr(,"class")
#> [1] "paleoTSfit"