MARSS: Fit a MARSS Model via Maximum-Likelihood Estimation

In MARSS: Multivariate Autoregressive State-Space Modeling

Fit a MARSS Model via Maximum-Likelihood Estimation

Description

This is the main function for fitting multivariate autoregressive state-space (MARSS) models with linear constraints. Scroll down to the bottom to see some short examples. To open a guide to show you how to get started quickly, type RShowDoc("Quick_Start",package="MARSS"). To open the MARSS User Guide from the command line, type RShowDoc("UserGuide",package="MARSS"). To get an overview of the package and all its main functions and how to get output (parameter estimates, fitted values, residuals, Kalmin filter or smoother output, or plots), go to MARSS-package. If MARSS() is throwing errors or warnings that you don't understand, try the Troubleshooting section of the user guide or type MARSSinfo() at the command line.

The default MARSS model form is "marxss", which is Multivariate Auto-Regressive(1) eXogenous inputs State-Space model:

\mathbf{x}_{t} = \mathbf{B}_t \mathbf{x}_{t-1} + \mathbf{u}_t + \mathbf{C}_t \mathbf{c}_t + \mathbf{G}_t \mathbf{w}_t, \textrm{ where } \mathbf{W}_t \sim \textrm{MVN}(0,\mathbf{Q}_t)

\mathbf{y}_t = \mathbf{Z}_t \mathbf{x}_t + \mathbf{a}_t + \mathbf{D}_t \mathbf{d}_t + \mathbf{H}_t \mathbf{v}_t, \textrm{ where } \mathbf{V}_t \sim \textrm{MVN}(0,\mathbf{R}_t)

\mathbf{X}_1 \sim \textrm{MVN}(\mathbf{x0}, \mathbf{V0}) \textrm{ or } \mathbf{X}_0 \sim \textrm{MVN}(\mathbf{x0}, \mathbf{V0})

The parameters are everything except \mathbf{x}, \mathbf{y}, \mathbf{v}, \mathbf{w}, \mathbf{c} and \mathbf{d}. \mathbf{y} are data (missing values allowed). \mathbf{c} and \mathbf{d} are inputs (no missing values allowed). All parameters (except \mathbf{x0} and \mathbf{V0}) can be time-varying but by default, all are time-constant (and the MARSS equation is generally written without the t subscripts on the parameter matrices). All parameters can be zero, including the variance matrices.

The parameter matrices can have fixed values and linear constraints. This is an example of a 3x3 matrix with linear constraints. All matrix elements can be written as a linear function of a, b, and c:

\left[\begin{array}{c c c} a+2b & 1 & a\\ 1+3a+b & 0 & b \\ 0 & -2 & c\end{array}\right]

Values such as a b or a^2 or log(a) are not linear constraints.

Usage

MARSS(y, 
    model = NULL, 
    inits = NULL, 
    miss.value = as.numeric(NA), 
    method = c("kem", "BFGS", "TMB", "BFGS_TMB", "nlminb_TMB"), 
    form = c("marxss", "dfa", "marss"), 
    fit = TRUE, 
    silent = FALSE, 
    control = NULL, 
    fun.kf = c("MARSSkfas", "MARSSkfss"), 
    ...) 

Arguments

The default settings for the optional arguments are set in MARSSsettings.R and are given below in the details section. For form specific defaults see the form help file (e.g. MARSS.marxss or MARSS.dfa).

y	A n x T matrix of n time series over T time steps. Only y is required for the function. A ts object (univariate or multivariate) can be used and this will be converted to a matrix with time in the columns.
inits	A list with the same form as the list outputted by coef(fit) that specifies initial values for the parameters. See also MARSS.marxss.
model	Model specification using a list of parameter matrix text shortcuts or matrices. See Details and MARSS.marxss for the default form. Or better yet open the Quick Start Guide RShowDoc("Quick_Start",package="MARSS").
miss.value	Deprecated. Denote missing values by NAs in your data.
method	Estimation method. MARSS provides an EM algorithm (method="kem") (see MARSSkem) and the BFGS algorithm (method="BFGS") (see MARSSoptim).
form	The equation form used in the MARSS() call. The default is "marxss". See MARSS.marxss or MARSS.dfa.
fit	TRUE/FALSE Whether to fit the model to the data. If FALSE, a marssMLE object with only the model is returned.
silent	Setting to TRUE(1) suppresses printing of full error messages, warnings, progress bars and convergence information. Setting to FALSE(0) produces error output. Setting silent=2 will produce more verbose error messages and progress information.
fun.kf	What Kalman filter function to use. MARSS has two: MARSSkfas() which is based on the Kalman filter in the KFAS package based on Koopman and Durbin and MARSSkfss() which is a native R implementation of the Kalman filter and smoother in Shumway and Stoffer. The KFAS filter is much faster. MARSSkfas() modifies the input and output in order to output the lag-one covariance smoother needed for the EM algorithm (per page 321 in Shumway and Stoffer (2000).
control	Estimation options for the maximization algorithm. The typically used control options for method="kem" are below but see marssMLE for the full list of control options. Note many of these are not allowed if method="BFGS"; see MARSSoptim for the allowed control options for this method. minit The minimum number of iterations to do in the maximization routine (if needed by method). If method="kem", this is an easy way to up the iterations and see how your estimates are converging. (positive integer) maxit Maximum number of iterations to be used in the maximization routine (if needed by method) (positive integer). min.iter.conv.test Minimum iterations to run before testing convergence via the slope of the log parameter versus log iterations. conv.test.deltaT=9 Number of iterations to use for the testing convergence via the slope of the log parameter versus log iterations. conv.test.slope.tol The slope of the log parameter versus log iteration to use as the cut-off for convergence. The default is 0.5 which is a bit high. For final analyses, this should be set lower. If you want to only use abstol as your convergence test, then to something very large, for example conv.test.slope.tol=1000. Type MARSSinfo(11) to see some comments on when you might want to do this. abstol The logLik.(iter-1)-logLik.(iter) convergence tolerance for the maximization routine. To meet convergence both the abstol and slope tests must be passed. allow.degen Whether to try setting \mathbf{Q} or \mathbf{R} elements to zero if they appear to be going to zero. trace An integer specifying the level of information recorded and error-checking run during the algorithms. trace=0, specifies basic error-checking and brief error-messages; trace>0 will print full error messages. In addition if trace>0, the Kalman filter output will be added to the outputted marssMLE object. Additional information recorded depends on the method of maximization. For the EM algorithm, a record of each parameter estimate for each EM iteration will be added. See optim for trace output details for the BFGS method. trace=-1 will turn off most internal error-checking and most error messages. The internal error checks are time expensive so this can speed up model fitting. This is particularly useful for bootstrapping and simulation studies. It is also useful if you get an error saying that MARSS() stops in MARSSkfss() due to a chol() call. MARSSkfss() uses matrix inversions and for some models these are unstable (high condition value). MARSSkfss() is used for error-checks and does not need to be called normally. safe Setting safe=TRUE runs the Kalman smoother after each parameter update rather than running the smoother only once after updated all parameters. The latter is faster but is not a strictly correct EM algorithm. In most cases, safe=FALSE (default) will not change the fits. If this setting does cause problems, you will know because you will see an error regarding the log-likelihood dropping and it will direct you to set safe=TRUE.
...	Optional arguments passed to function specified by form.

Details

The model argument specifies the structure of your model. There is a one-to-one correspondence between how you would write your model in matrix form on the whiteboard and how you specify the model for MARSS(). Many different types of multivariate time-series models can be converted to the MARSS form. See the User Guide and Quick Start Guide for examples.

The MARSS package has two forms for standard users: marxss and dfa.

MARSS.marxss

This is the default form. This is a MARSS model with (optional) inputs \mathbf{c}_t or \mathbf{d}_t. Most users will want this help page.

MARSS.dfa

This is a model form to allow easier specification of models for Dynamic Factor Analysis. The \mathbf{Z} parameters has a specific form and the \mathbf{Q} is set at i.i.d (diagonal) with variance of 1.

Those looking to modify or understand the base code, should look at MARSS.marss and MARSS.vectorized. These describe the forms used by the base functions. The EM algorithm uses the MARSS model written in vectorized form. This form is what allows linear constraints.

The likelihood surface for MARSS models can be multimodal or with strong ridges. It is recommended that for final analyses the estimates are checked by using a Monte Carlo initial conditions search; see the chapter on initial conditions searches in the User Guide. This requires more computation time, but reduces the chance of the algorithm terminating at a local maximum and not reaching the true MLEs. Also it is wise to check the EM results against the BFGS results (if possible) if there are strong ridges in the likelihood. Such ridges seems to slow down the EM algorithm considerably and can cause the algorithm to report convergence far from the maximum-likelihood values. EM steps up the likelihood and the convergence test is based on the rate of change of the log-likelihood in each step. Once on a strong ridge, the steps can slow dramatically. You can force the algorithm to keep working by setting minit. BFGS seems less hindered by the ridges but can be prodigiously slow for some multivariate problems. BFGS tends to work better if you give it good initial conditions (see Examples below for how to do this).

If you are working with models with time-varying parameters, it is important to notice the time-index for the parameters in the process equation (the \mathbf{x} equation). In some formulations (e.g. in KFAS), the process equation is \mathbf{x}_t=\mathbf{B}_{t-1}\mathbf{x}_{t-1}+\mathbf{w}_{t-1} so \mathbf{B}_{t-1} goes with \mathbf{x}_t not \mathbf{B}_t. Thus one needs to be careful to line up the time indices when passing in time-varying parameters to MARSS(). See the User Guide for examples.

Value

An object of class marssMLE. The structure of this object is discussed below, but if you want to know how to get specific output (like residuals, coefficients, smoothed states, confidence intervals, etc), see print.marssMLE(), tidy.marssMLE(), MARSSresiduals() and plot.marssMLE().

The outputted marssMLE object has the following components:

model	MARSS model specification. It is a marssMODEL object in the form specified by the user in the MARSS() call. This is used by print functions so that the user sees the expected form.
marss	The marssMODEL object in marss form. This form is needed for all the internal algorithms, thus is a required part of a marssMLE object.
call	All the information passed in in the MARSS() call.
start	List with specifying initial values that were used for each parameter matrix.
control	A list of estimation options, as specified by arguments control.
method	Estimation method.

If fit=TRUE, the following are also added to the marssMLE object. If fit=FALSE, a marssMLE object ready for fitting via the specified method is returned.

par	A list of estimated parameter values in marss form. Use print(), tidy() or coef() for outputing the model estimates in the MARSS() call (e.g. the default "marxss" form).
states	The expected value of \mathbf{X} conditioned on all the data, i.e. smoothed states.
states.se	The standard errors of the expected value of \mathbf{X}.
ytT	The expected value of \mathbf{Y} conditioned on all the data. Note this is just y for those y that are not missing.
ytT.se	The standard errors of the expected value of \mathbf{Y}. Note this is 0 for any non-missing y.
numIter	Number of iterations required for convergence.
convergence	Convergence status. 0 means converged successfully, 3 means all parameters were fixed (so model did not need to be fit) and -1 means call was made with fit=FALSE and parameters were not fixed (thus no $par element and Kalman filter/smoother cannot be run). Anything else is a warning or error. 2 means the marssMLE object has an error; the object is returned so you can debug it. The other numbers are errors during fitting. The error code depends on the fitting method. See MARSSkem and MARSSoptim.
logLik	Log-likelihood.
AIC	Akaike's Information Criterion.
AICc	Sample size corrected AIC.

If control$trace is set to 1 or greater, the following are also added to the marssMLE object.

kf	A list containing Kalman filter/smoother output from MARSSkf(). This is not normally added to a marssMLE object since it is verbose, but can be added using MARSSkf().
Ey	A list containing output from MARSShatyt. This isn't normally added to a marssMLE object since it is verbose, but can be computed using MARSShatyt().

Author(s)

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

References

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 Type RShowDoc("UserGuide",package="MARSS") to open a copy.

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]

Holmes, E. E., E. J. Ward and K. Wills. (2012) MARSS: Multivariate autoregressive state-space models for analyzing time-series data. R Journal 4: 11-19.

See Also

marssMLE, MARSSkem(), MARSSoptim(), MARSSkf(), MARSS-package, print.marssMLE(), plot.marssMLE(), print.marssMODEL(), MARSS.marxss(), MARSS.dfa(), fitted(), residuals(), MARSSresiduals(), predict(), tsSmooth(), tidy(), coef()

Examples

dat <- t(harborSealWA)
dat <- dat[2:4, ] # remove the year row
# fit a model with 1 hidden state and 3 observation time series
kemfit <- MARSS(dat, model = list(
  Z = matrix(1, 3, 1),
  R = "diagonal and equal"
))
kemfit$model # This gives a description of the model
print(kemfit$model) # same as kemfit$model
summary(kemfit$model) # This shows the model structure

# add CIs to a marssMLE object
# default uses an estimated Hessian matrix
kem.with.hess.CIs <- MARSSparamCIs(kemfit)
kem.with.hess.CIs

# fit a model with 3 hidden states (default)
kemfit <- MARSS(dat, silent = TRUE) # suppress printing
kemfit

# Fit the above model with BFGS using a short EM fit as initial conditions
kemfit <- MARSS(dat, control=list(minit=5, maxit=5))
bffit <- MARSS(dat, method="BFGS", inits=kemfit)

# fit a model with 3 correlated hidden states
# with one variance and one  covariance
# maxit set low to speed up example, but more iters are needed for convergence
kemfit <- MARSS(dat, model = list(Q = "equalvarcov"), control = list(maxit = 50))
# use Q="unconstrained" to allow different variances and covariances

# fit a model with 3 independent hidden states
# where each observation time series is independent
# the hidden trajectories 2-3 share their U parameter
kemfit <- MARSS(dat, model = list(U = matrix(c("N", "S", "S"), 3, 1)))

# same model, but with fixed independent observation errors
# and the 3rd x processes are forced to have a U=0
# Notice how a list matrix is used to combine fixed and estimated elements
# all parameters can be specified in this way using list matrices
kemfit <- MARSS(dat, model = list(U = matrix(list("N", "N", 0), 3, 1), R = diag(0.01, 3)))

# fit a model with 2 hidden states (north and south)
# where observation time series 1-2 are north and 3 is south
# Make the hidden state process independent with same process var
# Make the observation errors different but independent
# Make the growth parameters (U) the same
# Create a Z matrix as a design matrix that assigns the "N" state to the first 2 rows of dat
# and the "S" state to the 3rd row of data
Z <- matrix(c(1, 1, 0, 0, 0, 1), 3, 2)
# You can use factor is a shortcut making the above design matrix for Z
# Z <- factor(c("N","N","S"))
# name the state vectors
colnames(Z) <- c("N", "S")
kemfit <- MARSS(dat, model = list(
  Z = Z,
  Q = "diagonal and equal", R = "diagonal and unequal", U = "equal"
))

# print the model followed by the marssMLE object
kemfit$model

## Not run: 
# simulate some new data from our fitted model
sim.data <- MARSSsimulate(kemfit, nsim = 10, tSteps = 10)

# Compute bootstrap AIC for the model; this takes a long, long time
kemfit.with.AICb <- MARSSaic(kemfit, output = "AICbp")
kemfit.with.AICb

## End(Not run)

## Not run: 
# Many more short examples can be found in the
# Quick Examples chapter in the User Guide
RShowDoc("UserGuide", package = "MARSS")

# You can find the R scripts from the chapters by
# going to the index page
RShowDoc("index", package = "MARSS")

## End(Not run)

MARSS documentation built on May 31, 2023, 9:28 p.m.