Parameter estimation for MARSS models using optim
Description
Parameter estimation for MARSS models using R's optim
function. This allows access to R's quasiNewton algorithms available via the optim
function. The MARSSoptim
function is called when MARSS
is called with method="BFGS"
. This is a base function in the MARSSpackage
. neglogLik
is a helper function for MARSSoptim
that returns the negative loglikelihood given a vector of the estimated parameters and a marssMLE
object. When possible, the Kalman filter and smoother functions from the KFAS R package are used.
Usage
1 2  MARSSoptim(MLEobj)
neglogLik(x, MLEobj)

Arguments
MLEobj 
An object of class 
x 
An vector of the estimated parameters as output by 
Details
Objects of class marssMLE
may be built from scratch but are easier to construct using MARSS
with MARSS(..., fit=FALSE, method="BFGS")
.
Options for optim
are passed in using MLEobj$control
. See optim
for a list of that function's control options. If lower
and upper
for optim
need to be passed in, they should be passed in as part of control
as control$lower
and control$upper
. Additional control
arguments affect printing and initial conditions.
MLEobj$control$MCInit
If TRUE, Monte Carlo initialization will be performed by
MARSSmcinit
.MLEobj$control$numInits
Number of random initial value draws to be used with
MARSSmcinit
. Ignored ifcontrol$MCInit=FALSE
.MLEobj$control$numInitSteps
Maximum number of EM iterations for each random initial value draw to be used with
MARSSmcinit
. Ignored ifcontrol$MCInit=FALSE
.MLEobj$control$boundsInits
Length 6 list. Each component is a length 2 vector of bounds on the uniform distributions from which initial values will be drawn to be used with
MARSSmcinit()
. Ignored ifcontrol$MCInit=FALSE
. See Examples.MLEobj$control$kf.x0
The initial condition is at $t=0$ if kf.x0="x00". The initial condition is at $t=1$ if kf.x0="x10".
MLEobj$marss$diffuse
If diffuse=TRUE, a diffuse initial condition is used. MLEobj$par$V0 is then the scaling function for the diffuse part of the prior. Thus the prior is V0*kappa where kappaâ€“>Inf. Note that setting a diffuse prior does not change the correlation structure within the prior. If diffuse=FALSE, a nondiffuse prior is used and MLEobj$par$V0 is the nondiffuse prior variance on the initial states. The the prior is V0.
MLEobj$control$silent
Suppresses printing of progress bars, error messages, warnings and convergence information.
Value
The marssMLE
object which was passed in, with additional components:
method 
String "BFGS". 
kf 
Kalman filter output. 
iter.record 
If 
numIter 
Number of iterations needed for convergence. 
convergence 
Did estimation converge successfully?

logLik 
Loglikelihood. 
states 
State estimates from the Kalman filter. 
states.se 
Confidence intervals based on state standard errors, see caption of Fig 6.3 (p. 337) Shumway & Stoffer. 
errors 
Any error messages. 
Discussion
The function only returns parameter estimates. To compute CIs, use MARSSparamCIs
but if you use parametric or nonparametric bootstrapping with this function, it will use the EM algorithm to compute the bootstrap parameter estimates! The quasiNewton estimates are too fragile for the bootstrap routine since one often needs to search to find a set of initial conditions that work (i.e. don't lead to numerical errors).
Estimates from MARSSoptim
(which come from optim
) should be checked against estimates from the EM algorithm. If the quasiNewton algorithm works, it will tend to find parameters with higher likelihood faster than the EM algorithm. However, the MARSS likelihood surface can be multimodal with sharp peaks at degenerate solutions where a Q or R diagonal element equals 0. The quasiNewton algorithm sometimes gets stuck on these peaks even when they are not the maximum. Neither an initial conditions search nor starting near the known maximum (or from the parameters estimates after the EM algorithm) will necessarily solve this problem. Thus it is wise to check against EM estimates to ensure that the BFGS estimates are close to the MLE estimates (and visaversa, it's wise to rerun with method="BFGS" after using method="kem"). Conversely, there is a strong flat ridge in your likelihood, the EM algorithm can report early convergence while the BFGS may continue much further along the ridge and find very different parameter values. Of course a likelihood surface with strong flat ridges makes the MLEs less informative...
Note this is mainly a problem if the time series are short or very gappy. If the time series are long, then the likelihood surface should be nice with a single interior peak. In this case, the quasiNewton algorithm works well but it can still be sensitive (and slow) if not started with a good initial condition. Thus starting it with the estimates from the EM algorithm is often desirable.
One should be aware that the prior set on the variance of the initial states at t=0 or t=1 can have catastrophic effects on one's estimates if the presumed prior covariance structure conflicts with the structure implied by the MARSS model. For example, if you use a diagonal variancecovariance matrix for the prior but the model implies a matrix with nonzero covariances, your MLE estimates can be strongly influenced by the prior variancecovariance matrix. Setting a diffuse prior does not help because the diffuse prior still has the correlation structure specified by V0. One way to detect priors effects is to compare the BFGS estimates to the EM estimates. Persistent differences typically signify a problem with the correlation structure in the prior conflicting with the implied correlation structure in the MARSS model. If this is the case, using V0=0 and estimating the x0 elements (with control$kf.x0="x10") can often help.
Author(s)
Eli Holmes, NOAA, Seattle, USA.
eli(dot)holmes(at)noaa(dot)gov
See Also
MARSS
MARSSkem
marssMLE
optim
Examples
1 2 3 4 5 6 7 8 9 10  dat = t(harborSealWA)
dat = dat[2:4,] #remove the year row
#fit a model with EM and then use that fit as the start for BFGS
#fit a model with 1 hidden state where obs errors are iid
#R="diagonal and equal" is the default so not specified
#Q is fixed
kemfit = MARSS(dat, model=list(Z=matrix(1,3,1),Q=matrix(.01)))
bfgsfit = MARSS(dat, model=list(Z=matrix(1,3,1),Q=matrix(.01)),
inits=coef(kemfit,form="marss"), method="BFGS")
