opt.joint.GRW {paleoTS} | R Documentation |
Functions to find maximum likelihood solutions to general random walk (opt.joint.GRW
), unbiased random walk (opt.joint.URW
), stasis (opt.joint.Stasis
) and OU models (opt.joint.OU
).
opt.joint.GRW(x, pool = TRUE, cl = list(fnscale = -1), meth = "L-BFGS-B", hess = FALSE) opt.joint.URW(x, pool = TRUE, cl = list(fnscale = -1), meth = "L-BFGS-B", hess = FALSE) opt.joint.Stasis(x, pool = TRUE, cl = list(fnscale = -1), meth = "L-BFGS-B", hess = FALSE) opt.joint.OU(x, pool = TRUE, cl = list(fnscale = -1), meth = "L-BFGS-B", hess = FALSE)
x |
a paleoTS object |
pool |
logical indicating whether to pool variances across samples |
cl |
control list, passed to function optim |
meth |
optimization method, passed to function optim |
hess |
logical, indicating whether to calculate standard errors from the Hessian matrix |
These functions numerically search a log-likelihood surface for its optimum–they are a convenient wrapper to optim
.
Arguments meth
, cl
, and hess
are passed to optim
; see the help for that function for details.
These are included to allow sophisticated users greater control over the optimization; the defaults seem to work well for most,
but not all sequences. For meth="L-BFGS-B"
, some parameters are constrained to be non-negative, which is useful parameters
which cannot truly be negative, such as vstep
(random walk) and omega
(stasis model).
Initial estimates to start the optimization come in part from analytical solutions based on assuming equal sampling error across
samples and evenly spaced samples in time (functions mle.GRW
, mle.URW
and mle.Stasis
).
par |
parameter estimates |
value |
the log-likelihood of the optimal solution |
counts |
returned by optim |
convergence |
returned by optim |
message |
returned by optim |
p0 |
initial guess for parameter values at start of optimization |
K |
number of parameters in the model |
n |
the number of observations, equal to the number of samples |
AIC |
Akaike information criterion |
AICc |
bias-corrected Akaike information criterion |
BIC |
Bayes (or Schwarz) information criterion |
se |
standard errors for parameter estimates, computed from the curvature of the log-likelihood surface (only if hess = TRUE ) |
... |
other output from call to optim |
Measures of model fit (log-likelihoods, AIC scores, etc) are not comparable between the two parameterizations.
These optimizations are performed using a parameterization of the GRW, URW and Stasis models that differs from the functions opt.GRW
,
opt.URW
, and opt.Stasis
. Those functions models are fit from the differences between adjacent samples,
removing the autocorrelation in the time-series. In the joint parameterization implemented in functions opt.alt.GRW
, opt.alt.URW
,
opt.alt.Stasis
and opt.alt.OU
, models are fit using the actual sample values, with the autocorrelation among samples accounted for in
the log-likelihood function. For each model, the joint distribution of sample means is multivariate normal, with means
and variance-covariances determined by evolutionary parameters and sampling errors. Note that the Orstein-Uhlenbeck model is at present only implemented
using the joint parameterization.
For details on this parameterization of the models, see the paper Hunt, et al (2008). In general,
the two different parameterizations yield similar results as to the relative support for different models. In my experience, the two approaches tend to differ appreciably
only with relatively long sequences that have small differences between consecutive samples. In these circumstances, the alternative parameterization is better able to
distinguish true evolutionary patterns from sampling error, and is less prone to falsely favoring the Stasis model.
The general random walk (GRW
) and unbiased random walk (URW
) models require an additional parameter that specifies the phenotype at the start of
the sequence. This parameter is called anc
in the vector of returned parameter estimates, and in general is rather similar
to the trait value of the first sample in the sequence. This extra parameter is not necessary for the stasis model because the initial value does not
figure into the likelihood calculations (which are comletely determined by the position of optimum and the variance around this optimum).
Gene Hunt
Hunt, G. 2006. Fitting and comparing models of phyletic evolution: random walks and beyond. Paleobiology32:578–601.
Hunt, G., M. Bell & M. Travis. 2008. Evolution towards a new adaptive optimum: phenotypic evolution in a fossil stickleback lineage. Evolution 62:700-710.
x<- sim.GRW(ns=30, ms=1, vs=1) plot(x) # easier to use fit3models.joint() m.urw<- opt.joint.URW(x) m.grw<- opt.joint.GRW(x) m.sta<- opt.joint.Stasis(x) cat(m.urw$AICc, m.grw$AICc, m.sta$AICc, "\n") # print AICc scores