EM Algorithm function for MARSS models

MARSSkem() performs maximum-likelihood estimation, using an EM algorithm for constrained and unconstrained MARSS models. Users would not call this function directly normally. The function MARSS() calls MARSSkem(). However users might want to use MARSSkem() directly if they need to avoid some of the error-checking overhead associated with the MARSS() function.

Usage

MARSSkem(MLEobj)

Arguments

MLEobj: An object of class marssMLE.

Details

Objects of class marssMLE may be built from scratch but are easier to construct using MARSS() with MARSS(..., fit=FALSE).

Options for MARSSkem() may be set using MLEobj$control. The commonly used elements of control are as follows (see marssMLE):

minit: Minimum number of EM iterations. You can use this to force the algorithm to do a certain number of iterations. This is helpful if your solution is not converging.
maxit: Maximum number of EM iterations.
min.iter.conv.test: The minimum number of iterations before the log-log convergence test will be computed. If maxit is set less than this, then convergence will not be computed (and the algorithm will just run for maxit iterations).
kf.x0: Whether to set the prior at $t=0$ ("x00") or at $t=1$ ("x10"). The default is "x00".
conv.test.deltaT: The number of iterations to use in the log-log convergence test. This defaults to 9.
abstol: Tolerance for log-likelihood change for the delta logLik convergence test. If log-likelihood changes less than this amount relative to the previous iteration, the EM algorithm exits. This is normally (default) set to NULL and the log-log convergence test is used instead.
allow.degen: Whether to try setting $\mathbf{Q}$ or $\mathbf{R}$ elements to zero if they appear to be going to zero.
trace: A positive integer. If not 0, a record will be created of each variable over all EM iterations and detailed warning messages (if appropriate) will be printed.
safe: If TRUE, MARSSkem will rerun MARSSkf after each individual parameter update rather than only after all parameters are updated. The latter is slower and unnecessary for many models, but in some cases, the safer and slower algorithm is needed because the ML parameter matrices have high condition numbers.
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 "kem".

kf

Kalman filter output.

iter.record

If MLEobj$control$trace = TRUE, a list with par = a record of each estimated parameter over all EM iterations and logLik = a record of the log likelihood at each iteration.

numIter

Number of iterations needed for convergence.

convergence

Did estimation converge successfully?

convergence=0: Converged in both the abstol test and the log-log plot test.
convergence=1: Some of the parameter estimates did not converge (based on the log-log plot test AND abstol tests) before MLEobj$control$maxit was reached. This is not an error per se.
convergence=3: No convergence diagnostics were computed because all parameters were fixed thus no fitting required.
convergence=-1: No convergence diagnostics were computed because the MLE object was not fit (called with fit=FALSE). This isn't a convergence error just information. There is not par element so no functions can be run with the object.
convergence=2: No convergence diagnostics were computed because the MLE object had problems and was not fit. This isn't a convergence error just information.
convergence=10: Abstol convergence only. Some of the parameter estimates did not converge (based on the log-log plot test) before MLEobj$control$maxit was reached. However MLEobj$control$abstol was reached.
convergence=11: Log-log convergence only. Some of the parameter estimates did not converge (based on the abstol test) before MLEobj$control$maxit was reached. However the log-log convergence test was passed.
convergence=12: Abstol convergence only. Log-log convergence test was not computed because MLEobj$control$maxit was set to less than control$min.iter.conv.test.
convergence=13: Lack of convergence info. Parameter estimates did not converge based on the abstol test before MLEobj$control$maxit was reached. No log-log information since control$min.iter.conv.test is less than MLEobj$control$maxit so no log-log plot test could be done.
convergence=42: MLEobj$control$abstol was reached but the log-log plot test returned NAs. This is an odd error and you should set control$trace=TRUE and look at the outputted $iter.record to see what is wrong.
convergence=52: The EM algorithm was abandoned due to numerical errors. Usually this means one of the variances either went to zero or to all elements being equal. This is not an error per se. Most likely it means that your model is not very good for your data (too inflexible or too many parameters). Try setting control$trace=1 to view a detailed error report.
convergence=53: The algorithm was abandoned due to numerical errors in the likelihood calculation from MARSSkf.
convergence=62: The algorithm was abandoned due to errors in the log-log convergence test. You should not get this error (it is included for debugging purposes to catch improper arguments passed into the log-log convergence test).
convergence=63: The algorithm was run for control$maxit iterations, control$abstol not reached, and the log-log convergence test returned errors. You should not get this error (it is included for debugging purposes to catch improper arguments passed into the log-log convergence test).
convergence=72: Other convergence errors. This is included for debugging purposes to catch misc. errors.

logLik

Log-likelihood.

states

State estimates from the Kalman smoother.

states.se

Confidence intervals based on state standard errors, see caption of Fig 6.3 (p. 337) in Shumway & Stoffer (2006).

errors

Any error messages.

Discussion

To ensure that the global maximum-likelihood values are found, it is recommended that you test the fit under different initial parameter values, particularly if the model is not a good fit to the data. This requires more computation time, but reduces the chance of the algorithm terminating at a local maximum and not reaching the true MLEs. For many models and for draft analyses, this is unnecessary, but answers should be checked using an initial conditions search before reporting final values. See the chapter on initial conditions in the User Guide for a discussion on how to do this.

MARSSkem() calls a Kalman filter/smoother MARSSkf() for hidden state estimation. The algorithm allows two options for the initial state conditions: fixed but unknown or a prior. In the first case, x0 (whether at t=0 or t=1) is treated as fixed but unknown (estimated); in this case, fixed$V0=0 and x0 is estimated. This is the default behavior. In the second case, the initial conditions are specified with a prior and V0!=0. In the later case, x0 or V0 may be estimated. MARSS will allow you to try to estimate both, but many researchers have noted that this is not robust so you should fix one or the other.

If you get errors, you can type MARSSinfo() for help. Fitting problems often mean that the solution involves an ill-conditioned matrix. For example, your $\mathbf{Q}$ or $\mathbf{R}$ matrix is going to a value in which all elements have the same value, for example zero. If for example, you tried to fit a model with a fixed $\mathbf{R}$ matrix with high values on the diagonal and the variance in that $\mathbf{R}$ matrix (diagonal terms) was much higher than what is actually in the data, then you might drive $\mathbf{Q}$ to zero. Also if you try to fit a structurally inadequate model, then it is not unusual that $\mathbf{Q}$ will be driven to zero. For example, if you fit a model with 1 hidden state trajectory to data that clearly have 2 quite different hidden state trajectories, you might have this problem. Comparing the likelihood of this model to a model with more structural flexibility should reveal that the structurally inflexible model is inadequate (much lower likelihood).

Convergence testing is done via a combination of two tests. The first test (abstol test) is the test that the change in the absolute value of the log-likelihood from one iteration to another is less than some tolerance value (abstol). The second test (log-log test) is that the slope of a plot of the log of the parameter value or log-likelihood versus the log of the iteration number is less than some tolerance. Both of these must be met to generate the Success! parameters converged output. If you want to circumvent one of these tests, then set the tolerance for the unwanted test to be high. That will guarantee that that test is met before the convergence test you want to use is met. The tolerance for the abstol test is set by control$abstol and the tolerance for the log-log test is set by control$conv.test.slope.tol. Anything over 1 is huge for both of these.

References

R. H. Shumway and D. S. Stoffer (2006). Chapter 6 in Time series analysis and its applications. Springer-Verlag, New York.

Ghahramani, Z. and Hinton, G. E. (1996) Parameter estimation for linear dynamical systems. Technical Report CRG-TR-96-2, University of Toronto, Dept. of Computer Science.

Harvey, A. C. (1989) Chapter 5 in Forecasting, structural time series models and the Kalman filter. Cambridge University Press, Cambridge, UK.

The MARSS User Guide: Holmes, E. E., E. J. Ward, and M. D. Scheuerell (2012) Analysis of multivariate time-series using the MARSS package. NOAA Fisheries, Northwest Fisheries Science Center, 2725 Montlake Blvd E., Seattle, WA 98112 Go to User Guide to open the most recent version.

Holmes, E. E. (2012). Derivation of the EM algorithm for constrained and unconstrained multivariate autoregressive state-space (MARSS) models. Technical Report. arXiv:1302.3919 [stat.ME] EMDerivation has the most recent version.

Author

Eli Holmes and Eric Ward, NOAA, Seattle, USA.

Examples

dat <- t(harborSeal)
dat <- dat[2:4, ]
# you can use MARSS to construct a proper marssMLE object.
fit <- MARSS(dat, model = list(Q = "diagonal and equal", U = "equal"), fit = FALSE)
#> 
#> Model form is marxss. Model Structure is
#> m: 3 state process(es) named X.CoastalEstuaries X.OlympicPeninsula X.StraitJuanDeFuca
#> n: 3 observation time series named CoastalEstuaries OlympicPeninsula StraitJuanDeFuca
#> 
#> Z :  design matrix with rows: X.CoastalEstuaries X.OlympicPeninsula X.StraitJuanDeFuca (3 x 3) 
#> A :  fixed and zero (3 x 1) 
#> R :  diagonal and equal (3 x 3) 
#> B :  identity (3 x 3) 
#> U :  all equal (3 x 1) 
#> Q :  diagonal and equal (3 x 3) 
#> x0 :  unconstrained (3 x 1) 
#> V0 :  fixed and zero (3 x 3) 
#> D :  fixed and zero (3 x 1) 
#> C :  fixed and zero (3 x 1) 
#> d :  fixed and zero (1 x 1) 
#> c :  fixed and zero (1 x 1) 
#> G :  identity (3 x 3) 
#> H :  identity (3 x 3) 
#> L :  identity (3 x 3) 
# Pass this marssMLE object to MARSSkem to do the fit.
kemfit <- MARSSkem(fit)
#> Success! abstol and log-log tests passed at 59 iterations.
#> Alert: conv.test.slope.tol is 0.5.
#> Test with smaller values (<0.1) to ensure convergence.