Package 'MARSS'

Multivariate Autoregressive State-Space Model Estimation

Description

The MARSS package fits time-varying constrained and unconstrained multivariate autoregressive time-series models to multivariate time series data. To get started quickly, go to the Quick Start Guide (or at the command line, you can type RShowDoc("Quick_Start", package="MARSS")). To open the MARSS User Guide with many vignettes and examples, go to User Guide (or type RShowDoc("UserGuide",package="MARSS")).

The default MARSS model form is a MARXSS model: Multivariate Auto-Regressive(1) eXogenous inputs State-Space model. This model has the following form:

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

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

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

All parameters can be time-varying; the $t$ subscript is left off the parameters to remove clutter. Note, by default $\mathbf{V0}$ is a matrix of all zeros and thus $\mathbf{x}_1$ or $\mathbf{x}_0$ is treated as an estimated parameter not a diffuse prior.

The parameter matrices can have fixed values and linear constraints. This is an example of a 3x3 matrix with fixed values and linear constraints. In this example all the 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 $ln(a)$ are not allowed as those would not be linear.

The MARSS model parameters, hidden state processes ($\mathbf{x}$), and observations ($\mathbf{y}$) are matrices:

• $\mathbf{x}_t$, $\mathbf{x0}$, and $\mathbf{u}$ are m x 1

• $\mathbf{y}_t$ and $\mathbf{a}$ are n x 1 (m<=n)

• $\mathbf{B}$ and $\mathbf{V0}$ are m x m

• $\mathbf{Z}$ is n x m

• $\mathbf{Q}$ is g x g (default m x m)

• $\mathbf{G}$ is m x g (default m x m identity matrix)

• $\mathbf{R}$ is h x h (default n x n)

• $\mathbf{H}$ is n x h (default n x n identity matrix)

• $\mathbf{C}$ is m x q

• $\mathbf{D}$ is n x p

• $\mathbf{c}_t$ is q x 1

• $\mathbf{d}_t$ is p x 1

If a parameter is time-varying then the time dimension is the 3rd dimension. Thus a time-varying $\mathbf{Z}$ would be n x m x T where T is the length of the data time series.

The main fitting function is MARSS() which is used to fit a specified model to data and estimate the model parameters. MARSS() estimates the model parameters using an EM algorithm (primarily but see MARSSoptim()). Functions are provided for parameter confidence intervals and the observed Fisher Information matrix, smoothed state estimates with confidence intervals, all the Kalman filter and smoother outputs, residuals and residual diagnostics, printing and plotting, and summaries.

Details

Main MARSS functions:

MARSS()

Top-level function for specifying and fitting MARSS models.

coef()

Returns the estimated parameters in a variety of formats.

tidy()

Parameter estimates with confidence intervals

tsSmooth()

$\mathbf{x}$ and $\mathbf{y}$ estimates output as a data frame. Output can be conditioned on all the data ($T$), data up to $t-1$, or data up to $t$. From the Kalman filter and smoother output.

fitted()

Model x$\mathbf{x}$ and $\mathbf{y}$ predictions as a data frame or matrices. Another user interface for model predictions is predict.marssMLE.

residuals()

Model residuals as a data frame.

MARSSresiduals()

Model residuals as a data frame or matrices. Normal user interface to residuals is residuals.marssMLE.

predict()

Predictions and forecasts from a marssMLE object.

plot for marssMLE

A series of plots of fits and residuals diagnostics.

autoplot()

A series of plots using ggplot2 of fits and residuals diagnostics.

glance()

Brief summary of fit.

logLik()

Log-likelihood.

print()

Prints a wide variety of output from a marssMLE object.

print.marssMODEL()

Prints description of the MARSS model (marssMODEL object).

plot.marssPredict()

Plot a prediction or forecast.

toLatex.marssMODEL()

Outputs a LaTeX version of the model.

Other outputs for a fitted model:

MARSSsimulate()

Produces simulated data from a MARSS model.

MARSSkf(), MARSSkfas(), MARSSkfss()

Kalman filters and smoothers with extensive output of all the intermediate filter and smoother variances and expectations.

MARSSaic()

Calculates AICc, AICc, and various bootstrap AICs.

MARSSparamCIs()

Adds confidence intervals to a marssMLE object.

MARSShessian()

Computes an estimate of the variance-covariance matrix for the MLE parameters.

MARSSFisherI()

Returns the observed Fisher Information matrix.

Important internal MARSS functions (called by the above functions):

MARSSkem()

Estimates MARSS parameters using an EM algorithm.

MARSSoptim()

Estimates MARSS parameters using a quasi-Newton algorithm via optim.

MARSShatyt()

Calculates the expectations involving Y.

MARSSinnovationsboot()

Creates innovations bootstrapped data.

MARSS.marss()

Discusses the form in which MARSS models are stored internally.

Use help.search("internal", package="MARSS") to see the documentation of all the internal functions in the MARSS R package.

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

The MARSS Quick Start Guide: Quick Start Guide or type RShowDoc("Quick_Start",package="MARSS") to open a copy.

---

Return accuracy metrics

Description

This is a method for the generic accuracy function in the generics package. It is written to mimic the output from the accuracy function in the forecast package. See that package for details.

The measures calculated are:

• ME: Mean Error

• RMSE: Root Mean Squared Error

• MAE: Mean Absolute Error

• MPE: Mean Percentage Error

• MAPE: Mean Absolute Percentage Error

• MASE: Mean Absolute Scaled Error

• ACF1: Autocorrelation of errors at lag 1.

The MASE calculation is scaled using MAE of the training set naive forecasts which are simply $\mathbf{y}_{t-1}$.

For the training data, the metrics are shown for the one-step-ahead predictions by default (type="ytt1"). This is the prediction of $\mathbf{y}_t$ conditioned on the data up to $t-1$ (and the model estimated from all the data). With type="ytT", you can compute the metrics for the fitted ytT, which is the expected value of new data at $t$ conditioned on all the data. type does not affect test data (forecasts are past the end of the training data).

Usage

## S3 method for class 'marssPredict'
accuracy(object, x, test = NULL, type = "ytt1", verbose = FALSE, ...)
## S3 method for class 'marssMLE'
accuracy(object, x, test = NULL, type = "ytt1", verbose = FALSE, ...)

Arguments

object	A marssMLE or marssPredict object
x	A matrix or data frame with data to test against the h steps of a forecast.
test	Which time steps in training data (data model fit to) to compute accuracy for.
type	type="ytt1" is the one-step-ahead predictions. type="ytT" is the fitted ytT predictions. The former are standardly used for training data prediction metrics.
verbose	Show metrics for each time series of data.
...	Not used.

References

Hyndman, R.J. and Koehler, A.B. (2006) "Another look at measures of forecast accuracy". International Journal of Forecasting, 22(4), 679-688.

Hyndman, R.J. and Athanasopoulos, G. (2018) "Forecasting: principles and practice", 2nd ed., OTexts, Melbourne, Australia. Section 3.4 "Evaluating forecast accuracy". https://otexts.com/fpp2/accuracy.html.

Examples

dat <- t(harborSeal)
dat <- dat[c(2, 11, 12),]
train.dat <- dat[,1:12]
fit <- MARSS(train.dat, model = list(Z = factor(c("WA", "OR", "OR"))))

accuracy(fit)

# Compare to test data set
fr <- predict(fit, n.ahead=10)
test.dat <- dat[,13:22]
accuracy(fr, x=test.dat)

---

Coefficient function for MARSS MLE objects

Description

MARSS() outputs marssMLE objects. coef(object), where object is the output from a MARSS() call, will print out the estimated parameters. The default output is a list with values for each parameter, however the output can be altered using the type argument to output a vector of all the estimated values (type="vector") or a list with the full parameter matrix with the estimated and fixed elements (type="matrix"). For a summary of the parameter estimates with CIs from the estimated Hessian, use try tidy(object).

Usage

## S3 method for class 'marssMLE'
coef(object, ..., type = "list", form = NULL, what = "par")

Arguments

object	A marssMLE object.
...	Other arguments. Not used.
type	What to output. Default is "list". Options are "list" A list of only the estimated values in each matrix. Each model matrix has it's own list element. "vector" A vector of all the estimated values in each matrix. "matrix" A list of the parameter matrices each parameter with fixed values at their fixed values and the estimated values at their estimated values. Time-varying parameters, including d and c in a marxss form model, are returned as an array with time in the 3rd dimension. parameter name Returns the parameter matrix for that parameter with fixed values at their fixed values and the estimated values at their estimated values. Note, time-varying parameters, including d and c in a marxss form model, are returned as an array with time in the 3rd dimension.
form	This argument can be ignored. By default, the model form specified in the call to MARSS() is used to determine how to display the coefficients. This information is in attr(object$model,"form") . The default form is "marxss"; see MARSS.marxss(). However, the internal functions convert this to form "marss"; see MARSS.marss(). The marss form of the model is stored (in object$marss). You can look at the coefficients in marss form by passing in form="marss".
what	By default, coef() shows the parameter estimates. Other options are "par.se", "par.lowCI", "par.upCI", "par.bias", and "start".

Value

A list of the estimated parameters for each model matrix.

Author(s)

Eli Holmes, NOAA, Seattle, USA.

See Also

tidy(), print()

Examples

dat <- t(harborSeal)
dat <- dat[c(2, 11), ]
fit <- MARSS(dat)

coef(fit)
coef(fit, type = "vector")
coef(fit, type = "matrix")
# to retrieve just the Q matrix
coef(fit, type = "matrix")$Q

---

Plot Extinction Risk Metrics

Description

Generates a six-panel plot of extinction risk metrics used in Population Viability Analysis (PVA). This is a function used by one of the vignettes in the MARSS-package.

Usage

CSEGriskfigure(data, te = 100, absolutethresh = FALSE, threshold = 0.1, 
  datalogged = FALSE, silent = FALSE, return.model = FALSE, 
  CI.method = "hessian", CI.sim = 1000)

Arguments

data	A data matrix with 2 columns; time in first column and counts in second column. Note time is down rows, which is different than the base MARSS-package functions.
te	Length of forecast period (positive integer)
absolutethresh	Is extinction threshold an absolute number? (T/F)
threshold	Extinction threshold either as an absolute number, if absolutethresh=TRUE, or as a fraction of current population count, if absolutethresh=FALSE.
datalogged	Are the data already logged? (T/F)
silent	Suppress printed output? (T/F)
return.model	Return state-space model as marssMLE object? (T/F)
CI.method	Confidence interval method: "hessian", "parametrc", "innovations", or "none". See MARSSparamCIs.
CI.sim	Number of simulations for bootstrap confidence intervals (positive integer).

Details

Panel 1: Time-series plot of the data. Panel 2: CDF of extinction risk. Panel 3: PDF of time to reach threshold. Panel 4: Probability of reaching different thresholds during forecast period. Panel 5: Sample projections. Panel 6: TMU plot (uncertainty as a function of the forecast).

Value

If return.model=TRUE, an object of class marssMLE.

Author(s)

Eli Holmes, NOAA, Seattle, USA, and Steve Ellner, Cornell Univ.

References

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.

(theory behind the figure) Holmes, E. E., J. L. Sabo, S. V. Viscido, and W. F. Fagan. (2007) A statistical approach to quasi-extinction forecasting. Ecology Letters 10:1182-1198.

(CDF and PDF calculations) Dennis, B., P. L. Munholland, and J. M. Scott. (1991) Estimation of growth and extinction parameters for endangered species. Ecological Monographs 61:115-143.

(TMU figure) Ellner, S. P. and E. E. Holmes. (2008) Resolving the debate on when extinction risk is predictable. Ecology Letters 11:E1-E5.

See Also

MARSSboot, marssMLE, CSEGtmufigure

Examples

d <- harborSeal[, 1:2]
kem <- CSEGriskfigure(d, datalogged = TRUE)

---

Plot Forecast Uncertainty

Description

Plot the uncertainty in the probability of hitting a percent threshold (quasi-extinction) for a single random walk trajectory. This is the quasi-extinction probability used in Population Viability Analysis. The uncertainty is shown as a function of the forecast, where the forecast is defined in terms of the forecast length (number of time steps) and forecasted decline (percentage). This is a function used by one of the vignettes in the MARSS-package.

Usage

CSEGtmufigure(N = 20, u = -0.1, s2p = 0.01, make.legend = TRUE)

Arguments

N	Time steps between the first and last population data point (positive integer)
u	Per time-step decline (-0.1 means a 10% decline per time step; 1 means a doubling per time step.)
s2p	Process variance (Q). (a positive number)
make.legend	Add a legend to the plot? (T/F)

Details

This figure shows the region of high uncertainty in dark grey. In this region, the minimum 95 percent confidence intervals on the probability of quasi-extinction span 80 percent of the 0 to 1 probability. Green hashing indicates where the 95 percent upper bound does not exceed 5% probability of quasi-extinction. The red hashing indicates, where the 95 percent lower bound is above 95% probability of quasi-extinction. The light grey lies between these two certain/uncertain extremes. The extinction calculation is based on Dennis et al. (1991). The minimum theoretical confidence interval is based on Fieberg and Ellner (2000). This figure was developed in Ellner and Holmes (2008).

Examples using this figure are shown in the User Guide in the PVA chapter.

Author(s)

Eli Holmes, NOAA, Seattle, USA, and Steve Ellner, Cornell Univ.

References

Dennis, B., P. L. Munholland, and J. M. Scott. (1991) Estimation of growth and extinction parameters for endangered species. Ecological Monographs 61:115-143.

Fieberg, J. and Ellner, S.P. (2000) When is it meaningful to estimate an extinction probability? Ecology, 81, 2040-2047.

Ellner, S. P. and E. E. Holmes. (2008) Resolving the debate on when extinction risk is predictable. Ecology Letters 11:E1-E5.

See Also

CSEGriskfigure

Examples

CSEGtmufigure(N = 20, u = -0.1, s2p = 0.01)

---

Example Data Sets

Description

Example data sets for use in MARSS vignettes for the MARSS-package.

• plankton Plankton data sets: Lake WA plankton 32-year time series and Ives et al data from West Long Lake.

• SalmonSurvCUI Snake River spring/summer chinook survival indices.

• isleRoyal Isle Royale wolf and moose data with temperature and precipitation covariates.

• population-count-data A variety of fish, mammal and bird population count data sets.

• loggerhead Loggerhead turtle tracking (location) data from ARGOS tags.

• harborSeal Harbor seal survey data (haul out counts) from Oregon, Washington and California, USA.

---

Return fitted values for X(t) and Y(t) in a MARSS model

Description

fitted() returns the different types of fitted values for $\mathbf{x}_t$ and $\mathbf{y}_t$ in a MARSS model. The fitted values are the expected value of the right side of the MARSS equations without the error terms, thus are the model predictions of $\mathbf{y}_t$ or $\mathbf{x}_t$. fitted.marssMLE is a companion function to tsSmooth() which returns the expected value of the right side of the MARSS equations with the error terms (the Kalman filter and smoother output).

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

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

The data go from $t=1$ to $t=T$. For brevity, the parameter matrices are shown without a time subscript, however all parameters can be time-varying.

Note that the prediction of $\mathbf{x}_t$ conditioned on the data up to time $t$ is not provided since that would require the expected value of $\mathbf{X}_{t}$ conditioned on data from $t = 1$ to $t+1$, which is not output from the Kalman filter or smoother.

Usage

## S3 method for class 'marssMLE'
fitted(object, ..., 
    type = c("ytt1", "ytT", "xtT", "ytt", "xtt1"),   
    interval = c("none", "confidence", "prediction"), 
    level = 0.95, 
    output = c("data.frame", "matrix"), 
    fun.kf = c("MARSSkfas", "MARSSkfss"))

Arguments

object	A marssMLE object.
type	If type="tT", then the predictions are conditioned on all the data. If type="tt", then the predictions are conditioned on data up to time $t$. If type="tt1", the predictions are conditioned on data up to time $t-1$. The latter are also known as one-step-ahead estimates. For $\mathbf{y}$, these are also known as the innovations.
interval	If interval="confidence", then the standard error and confidence interval of the predicted value is returned. If interval="prediction", then the standard deviation and prediction interval of new data or states are returned.
level	Level for the intervals if interval is not equal to "none".
output	data frame or list of matrices
fun.kf	By default, tsSmooth() will use the Kalman filter/smoother function in object$fun.kf (either MARSSkfas() or MARSSkfss()). You can pass in fun.kf to force a particular Kalman filter/smoother function to be used.
...	Not used.

Details

In the state-space literature, the two most commonly used fitted values are "ytt1" and "ytT". The former is the expected value of $\mathbf{Y}_t$ conditioned on the data 1 to time $t-1$. These are known as the innovations and they, plus their variance, are used in the calculation of the likelihood of a MARSS model via the innovations form of the likelihood. The latter, "ytT" are the model estimates of the $\mathbf{y}$ values using all the data; this is the expected value of $\mathbf{Z}\mathbf{X}_t+\mathbf{a}+\mathbf{D}\mathbf{d}_t$ conditioned on the data 1 to $T$. The "xtT" along with "ytT" are used for computing smoothation residuals used in outlier and shock detection. See MARSSresiduals. For completeness, fitted.marssMLE will also return the other possible model predictions with different conditioning on the data (1 to $t-1$, $t$, and $T$), however only type="ytt1" (innovations) and "ytT" and "xtT" (smoothations) are regularly used. Keep in mind that the fitted "xtT" is not the smoothed estimate of $\mathbf{x}$ (unlike "ytT"). If you want the smoothed estimate of $\mathbf{x}$ (i.e. the expected value of $\mathbf{X}_t$ conditioned on all the data), you want the Kalman smoother. See tsSmooth.

Fitted versus estimated values: The fitted states at time $t$ are predictions from the estimated state at time $t-1$ conditioned on the data: expected value of $\mathbf{B}\mathbf{X}_{t-1}+\mathbf{u}+\mathbf{C}\mathbf{c}_t$ conditioned on the data. They are distinguished from the estimated states at time $t$ which would includes the conditional expected values of the error terms: $\textrm{E}[\mathbf{X}_{t}] = \mathbf{B}\mathbf{X}_{t-1}+\mathbf{u}+\mathbf{C}\mathbf{c}_t + \mathbf{w}_t$. The latter are returned by the Kalman filter and smoother. Analogously, the fitted observations at time $t$ are model predictions from the estimated state at time $t$ conditioned on the data: the expected value of the right side of the $\mathbf{y}$ equation without the error term. Like with the states, one can also compute the expected value of the observations at time $t$ conditioned on the data: the expected value of the right side of the $\mathbf{y}$ equation with the error term. The latter would just be equal to the data if there are no missing data, but when there are missing data, this is what is used to estimate their values. The expected value of states and observations are provided via tsSmooth.

observation fitted values

The model predicted $\hat{\mathbf{y}}_t$ is $\mathbf{Z}\mathbf{x}_t^\tau+\mathbf{a} + \mathbf{D}\mathbf{d}_t$, where $\mathbf{x}_t^\tau$ is the expected value of the state at time $t$ conditioned on the data from 1 to $\tau$ ($\tau$ will be $t-1$, $t$ or $T$). Note, if you are using MARSS for estimating the values for missing data, then you want to use tsSmooth() with type="ytT" not fitted(..., type="ytT").

$\mathbf{x}_t^\tau$ is the expected value of the states at time $t$ conditioned on the data from time 1 to $\tau$. If type="ytT", the expected value is conditioned on all the data, i.e. xtT returned by MARSSkf(). If type="ytt1", then the expected value uses only the data up to time $t-1$, i.e. xtt1 returned by MARSSkf(). These are commonly known as the one step ahead predictions for a state-space model. If type="ytt", then the expected value uses the data up to time $t$, i.e. xtt returned by MARSSkf().

If interval="confidence", the standard error and interval is for the predicted $\mathbf{y}$. The standard error is $\mathbf{Z} \mathbf{V}_t^\tau \mathbf{Z}^\top$. If interval="prediction", the standard deviation of new iid $\mathbf{y}$ data sets is returned. The standard deviation of new $\mathbf{y}$ is $\mathbf{Z} \mathbf{V}_t^\tau \mathbf{Z}^\top + \mathbf{R}_t$. $\mathbf{V}_t^\tau$ is conditioned on the data from $t=1$ to $n$. $\tau$ will be either $t$, $t-1$ or $T$ depending on the value of type.

Intervals returned by predict() are not for the data used in the model but rather new data sets. To evaluate the data used to fit the model for residuals analysis or analysis or model inadequacy, you want the model residuals (and residual se's). Use residuals for model residuals and their intervals. The intervals for model residuals are narrower because the predictions for $\mathbf{y}$ were estimated from the model data (so is closer to the data used to estimate the predictions than new independent data will be).

state fitted values

The model predicted $\mathbf{x}_t$ given $\mathbf{x}_{t-1}$ is $\mathbf{B}\mathbf{x}_{t-1}^\tau+\mathbf{u}+\mathbf{C}\mathbf{c}_t$. If you want estimates of the states, rather than the model predictions based on $\mathbf{x}_{t-1}$, go to tsSmooth(). Which function you want depends on your objective and application.

$\mathbf{x}_{t-1}^\tau$ used in the prediction is the expected value of the states at time $t-1$ conditioned on the data from $t=1$ to $t=\tau$. If type="xtT", this is the expected value at time $t-1$ conditioned on all the data, i.e. xtT[,t-1] returned by MARSSkf(). If type="xtt1", it is the expected value conditioned on the data up to time $t-1$, i.e. xtt[,t-1] returned by MARSSkf(). The predicted state values conditioned on data up to $t$ are not provided. This would require the expected value of states at time $t$ conditioned on data up to time $t+1$, and this is not output by the Kalman filter. Only conditioning on data up to $t-1$ and $T$ are output.

If interval="confidence", the standard error of the predicted values (meaning the standard error of the expected value of $\mathbf{X}_t$ given $\mathbf{X}_{t-1}$) is returned. The standard error of the predicted value is $\mathbf{B} \mathbf{V}_{t-1}^\tau\mathbf{B}^\top$. If interval="prediction", the standard deviation of $\mathbf{X}_t$ given $\mathbf{X}_{t-1}$ is output. The latter is $\mathbf{B} \mathbf{V}_{t-1}^\tau \mathbf{B}^\top + \mathbf{Q}$ . $\mathbf{V}_{t-1}^\tau$ is either conditioned on data 1 to $\tau=T$ or $\tau=t-1$ depending on type.

The intervals returned by fitted.marssMLE() are for the state predictions based on the state estimate at $t-1$. These are not typically what one uses or needs–however might be useful for simulation work. If you want confidence intervals for the state estimates, those are provided in tsSmooth. If you want to do residuals analysis (for outliers or model inadequacy), you want the residuals intervals provided in residuals().

Value

If output="data.frame" (the default), a data frame with the following columns is returned. If output="matrix", the columns are returned as matrices in a list. Information computed from the model has a leading "." in the column name.

If interval="none", the following are returned (colnames with . in front are computed values):

.rownames	Names of the data or states.
t	Time step.
y	The data if type is "ytT", "ytt" or "ytt1".
.x	The expected value of $\mathbf{X}_t$ conditioned on all the data if type="xtT" or data up to time $t$ if type="xtt1". From tsSmooth(). This is the expected value of the right-side of the $\mathbf{x}_t$ equation with the errors terms while .fitted is the expected value of the right side without the error term $\mathbf{w}_t$.
.fitted	Predicted values of observations ($\mathbf{y}$) or the states ($\mathbf{x}$). See details.

If interval = "confidence", the following are also returned:

.se	Standard errors of the predictions.
.conf.low	Lower confidence level at alpha = 1-level. The interval is approximated using qnorm(alpha/2)*.se + .fitted
.conf.up	Upper confidence level. The interval is approximated using qnorm(1-alpha/2)*.se + .fitted

The confidence interval is for the predicted value, i.e. $\mathbf{Z}\mathbf{x}_t^\tau+\mathbf{a}$ for $\mathbf{y}$ or $\mathbf{B}\mathbf{x}_{t-1}^\tau+\mathbf{u}$ for $\mathbf{x}$ where $\mathbf{x}_t^\tau$ is the expected value of $\mathbf{X}_t$ conditioned on the data from 1 to $\tau$. ($\tau$ will be $t-1$, $t$ or $T$).

If interval = "prediction", the following are also returned:

.sd	Standard deviation of new $\mathbf{x}_t$ or $\mathbf{y}_t$ iid values.
.lwr	Lower range at alpha = 1-level. The interval is approximated using qnorm(alpha/2)*.sd + .fitted
.upr	Upper range at level. The interval is approximated using qnorm(1-alpha/2)*.sd + .fitted

The prediction interval is for new $\mathbf{x}_t$ or $\mathbf{y}_t$. If you want to evaluate the observed data or the states estimates for outliers then these are not the intervals that you want. For that you need the residuals intervals provided by residuals().

Author(s)

Eli Holmes, NOAA, Seattle, USA.

See Also

MARSSkf(), MARSSresiduals(), residuals(), predict(), tsSmooth()

Examples

dat <- t(harborSeal)
dat <- dat[c(2, 11, 12), ]
fit <- MARSS(dat, model = list(Z = factor(c("WA", "OR", "OR"))))
fitted(fit)

# Example of fitting a stochastic level model to the Nile River flow data
# red line is smooothed level estimate
# grey line is one-step-ahead prediction using xtt1
nile <- as.vector(datasets::Nile)
mod.list <- list(
  Z = matrix(1), A = matrix(0), R = matrix("r"),
  B = matrix(1), U = matrix(0), Q = matrix("q"),
  x0 = matrix("pi")
)
fit <- MARSS(nile, model = mod.list, silent = TRUE)

# Plotting
# There are plot methods for marssMLE objects
library(ggplot2)
autoplot(fit)

# Below shows how to make plots manually but all of these can be made
# with autoplot(fit) or plot(fit)
plot(nile, type = "p", pch = 16, col = "blue")
lines(fitted(fit, type="ytT")$.fitted, col = "red", lwd = 2)
lines(fitted(fit, type="ytt1")$.fitted, col = "grey", lwd = 2)

# Make a plot of the model estimate of y(t), i.e., put a line through the points
# Intervals are for new data not the blue dots 
# (which were used to fit the model so are not new)
library(ggplot2)
d <- fitted(fit, type = "ytT", interval="confidence", level=0.95)
d2 <- fitted(fit, type = "ytT", interval="prediction", level=0.95)
d$.lwr <- d2$.lwr
d$.upr <- d2$.upr
ggplot(data = d) +
  geom_line(aes(t, .fitted), linewidth=1) +
  geom_point(aes(t, y), color = "blue", na.rm=TRUE) +
  geom_ribbon(aes(x = t, ymin = .conf.low, ymax = .conf.up), alpha = 0.3) +
  geom_line(aes(t, .lwr), linetype = 2) +
  geom_line(aes(t, .upr), linetype = 2) +
  facet_grid(~.rownames) +
  xlab("Time Step") + ylab("Count") +
  ggtitle("Blue=data, Black=estimate, grey=CI, dash=prediction interval") +
  geom_text(x=15, y=7, label="The intervals are for \n new data not the blue dots")

---

forecast function for marssMLE objects

Description

MARSS() outputs marssMLE objects. forecast(object), where object is marssMLE object, will return the forecasts of $\mathbf{y}_t$ or $\mathbf{x}_t$ for h steps past the end of the model data. forecast(object) returns a marssPredict object which can be passed to plot.marssPredict for automatic plotting of the forecast. forecast.marssMLE() is used by predict.marssMLE() to generate forecasts.

This is a method for the generic forecast function in the generics package. It is written to mimic the behavior and look of the forecast package.

Usage

## S3 method for class 'marssMLE'
forecast(object, h = 10, level = c(0.8, 0.95), 
     type = c("ytT","xtT", "ytt", "ytt1", "xtt", "xtt1"), 
     newdata = list(y = NULL, c = NULL, d = NULL), 
     interval = c("prediction", "confidence", "none"), 
     fun.kf = c("MARSSkfas", "MARSSkfss"), ...)

Arguments

object	A marssMLE object.
h	Number of steps ahead to forecast. newdata is for the forecast, i.e. for the $h$ time steps after the end of the model data. If there are covariates in the model, $\mathbf{c}_t$ or $\mathbf{d}_t$, then newdata is required. See details.
level	Level for the intervals if interval != "none".
type	The default for observations would be type="ytT" and for the states would be type="xtT", i.e. using all the data. Other possible forecasts are provided for completeness but would in most cases be identical (see details).
newdata	An optional list with matrices for new covariates $\mathbf{c}_t$ or $\mathbf{d}_t$ to use for the forecasts. $\mathbf{c}_t$ or $\mathbf{d}_t$ must be in the original model and have the same matrix rows and columns as used in the MARSS() call but the number of time steps can be different (and should be equal to h).
interval	If interval="confidence", then the standard error and confidence interval of the expected value of $\mathbf{y}_t$ (type="ytT") or $\mathbf{x}_t$ (type="xtT") is returned. interval="prediction" (default) returns prediction intervals which include the uncertainty in the expected value and due to observation error (the $\mathbf{R}$ in the $\mathbf{y}$ equation). Note, in the context of a MARSS model, only confidence intervals are available for the states (the $\mathbf{x}$).
fun.kf	Only if you want to change the default Kalman filter. Can be ignored.
...	Other arguments. Not used.

Details

The type="ytT" forecast for $T+i$ is

$\mathbf{Z}\mathbf{x}_{T+i}^T + \mathbf{a} + \mathbf{D}\mathbf{d}_{T+i}$

where $\mathbf{Z}$, $\mathbf{a}$ and $\mathbf{D}$ are estimated from the data from $t=1$ to $T$. If the model includes $\mathbf{d}_t$ then newdata with d must be passed in. Either confidence or prediction intervals can be shown. Prediction intervals would be the norm for forecasts and show the intervals for new data which based on the conditional variance of $\mathbf{Z}\mathbf{X}_{T+i} + \mathbf{V}_{T+i}$. Confidence intervals would show the variance of the mean of the new data (such as if you ran a simulation multiple times and recorded only the mean observation time series). It is based on the conditional variance of $\mathbf{Z}\mathbf{X}_{T+i}$. The intervals shown are computed with fitted().

The type="xtT" forecast for $T+i$ is

$\mathbf{B}\mathbf{x}_{T+i-1}^T + \mathbf{u} + \mathbf{C}\mathbf{c}_{T+i}$

where $\mathbf{B}$ and $\mathbf{u}$ and $\mathbf{C}$ are estimated from the data from $t=1$ to $T$ (i.e. the estimates in the marssMLE object). If the model includes $\mathbf{c}_t$ then newdata with c must be passed in. The only intervals are confidence intervals which based on the conditional variance of $\mathbf{B}\mathbf{X}_{T+i-1} + \mathbf{W}_{T+i}$. If you pass in data for your forecast time steps, then the forecast will be computed conditioned on the original data plus the data in the forecast period. The intervals shown are computed with the Kalman smoother (or filter if type="xtt" or type="xtt1" specified) via tsSmooth().

If the model has time-varying parameters, the parameter estimates at time $T$ will be used for the whole forecast. If new data c or d are passed in, it must have h time steps.

Note: y in newdata. Data along with covariates can be passed into newdata. In this case, the data in newdata ($T+1$ to $T+h$) are conditioned on for the expected value of $\mathbf{X}_t$ but parameters used are only estimated using the data in the marssMLE object ($t=1$ to $T$). If you include data in newdata, you need to decide how to condition on that new data for the forecast. type="ytT" would mean that the $t=T+i$ forecast is conditioned on all the data, $t=1$ to $T+h$, type="ytt" would mean that the $t=T+i$ forecast is conditioned on the data, $t=1$ to $T+i$, and type="ytt1" would mean that the $t=T+i$ forecast is conditioned on the data, $t=1$ to $T+i-1$. Because MARSS models can be used in all sorts of systems, the $\mathbf{y}$ part of the MARSS model might not be "data" in the traditional sense. In some cases, one of the $\mathbf{y}$ (in a multivariate model) might be a known deterministic process or it might be a simulated future $\mathbf{y}$ that you want to include. In this case the $\mathbf{y}$ rows that are being forecasted are NAs and the $\mathbf{y}$ rows that are known are passed in with newdata.

Value

A list with the following components:

method	The method used for fitting, e.g. "kem".
model	The marssMLE object passed into forecast.marssMLE().
newdata	The newdata list if passed into forecast.marssMLE().
level	The confidence level passed into forecast.marssMLE().
pred	A data frame the forecasts along with the intervals.
type	The type ("ytT" or "xtT") passed into forecast.marssMLE().
t	The time steps used to fit the model (used for plotting).
h	The number of forecast time steps (used for plotting).

Author(s)

Eli Holmes, NOAA, Seattle, USA.

See Also

plot.marssPredict(), predict.marssMLE()

Examples

# More examples are in ?predict.marssMLE

dat <- t(harborSealWA)
dat <- dat[2:4,] #remove the year row
fit <- MARSS(dat, model=list(R="diagonal and equal"))

# 2 steps ahead forecast
fr <- forecast(fit, type="ytT", h=2)
plot(fr)

# forecast and only show last 10 steps of original data
fr <- forecast(fit, h=10)
plot(fr, include=10)

---

Return brief summary information on a MARSS fit

Description

This returns a data frame with brief summary information.

coef.det

The coefficient of determination. This is the squared correlation between the fitted values and the original data points. This is simply a metric for the difference between the data points and the fitted values and should not be used for formal model comparison.

sigma

The sample variance (unbiased) of the data residuals (fitted minus data). This is another simple metric of the difference between the data and fitted values. This is different than the sigma returned by an arima() call for example. That sigma would be akin to sqrt(Q) in the MARSS parameters; 'akin' because MARSS models are multivariate and the sigma returned by arima() is for a univariate model.

df

The number of estimated parameters. Denoted num.params in a marssMLE object.

logLik

The log-likelihood.

AIC

Akaike information criterion.

AICc

Akaike information criterion corrected for small sample size.

AICbb

Non-parametric bootstrap Akaike information criterion if in the marssMLE object.

AICbp

Parametric bootstrap Akaike information criterion if in the marssMLE object.

convergence

0 if converged according to the convergence criteria set. Note the default convergence criteria are high in order to speed up fitting. A number other than 0 means the model did not meet the convergence criteria.

errors

0 if no errors. 1 if some type of error or warning returned.

Usage

## S3 method for class 'marssMLE'
glance(x, ...)

Arguments

x	A marssMLE object
...	Not used.

Examples

dat <- t(harborSeal)
dat <- dat[c(2, 11, 12), ]
fit <- MARSS(dat, model = list(Z = factor(c("WA", "OR", "OR"))))

glance(fit)

---

Harbor Seal Population Count Data (Log counts)

Description

Data sets used in MARSS vignettes in the MARSS-package. These are data sets based on logged count data from Oregon, Washington and California sites where harbor seals were censused while hauled out on land. "harborSealnomiss" is an extrapolated data set where missing values in the original data set have been extrapolated so that the data set can be used to demonstrate fitting population models with different underlying structures.

Usage

data(harborSeal)
data(harborSealWA)

Format

Matrix "harborSeal" contains columns "Years", "StraitJuanDeFuca", "SanJuanIslands", "EasternBays", "PugetSound", "HoodCanal", "CoastalEstuaries", "OlympicPeninsula", "CA.Mainland", "OR.NorthCoast", "CA.ChannelIslands", and "OR.SouthCoast".

Matrix "harborSealnomiss" contains columns "Years", "StraitJuanDeFuca", "SanJuanIslands", "EasternBays", "PugetSound", "HoodCanal", "CoastalEstuaries", "OlympicPeninsula", "OR.NorthCoast", and "OR.SouthCoast". Matrix "harborSealWA" contains columns "Years", "SJF", "SJI", "EBays", "PSnd", and "HC", representing the same five sites as the first five columns of "harborSeal".

Details

Matrix "harborSealWA" contains the original 1978-1999 logged count data for five inland WA sites. Matrix "harborSealnomiss" contains 1975-2003 data for the same sites as well as four coastal sites, where missing values have been replaced with extrapolated values. Matrix "harborSeal" contains the original 1975-2003 LOGGED data (with missing values) for the WA and OR sites as well as a CA Mainland and CA ChannelIslands time series.

Source

Jeffries et al. 2003. Trends and status of harbor seals in Washington State: 1978-1999. Journal of Wildlife Management 67(1):208-219.

Brown, R. F., Wright, B. E., Riemer, S. D. and Laake, J. 2005. Trends in abundance and current status of harbor seals in Oregon: 1977-2003. Marine Mammal Science, 21: 657-670.

Lowry, M. S., Carretta, J. V., and Forney, K. A. 2008. Pacific harbor seal census in California during May-July 2002 and 2004. California Fish and Game 94:180-193.

Hanan, D. A. 1996. Dynamics of Abundance and Distribution for Pacific Harbor Seal, Phoca vitulina richardsi, on the Coast of California. Ph.D. Dissertation, University of California, Los Angeles. 158pp. DFO. 2010. Population Assessment Pacific Harbour Seal (Phoca vitulina richardsi). DFO Can. Sci. Advis. Sec. Sci. Advis. Rep. 2009/011.

Examples

str(harborSealWA)
str(harborSeal)

---

Isle Royale Wolf and Moose Data

Description

Example data set for estimation of species interaction strengths. These are data on the number of wolves and moose on Isle Royal, Michigan. The data are unlogged. The covariate data are the average Jan-Feb, average Apr-May and average July-Sept temperature (Fahrenheit) and precipitation (inches). Also included are 3-year running means of these covariates. The choice of covariates is based on those presented in the Isle Royale 2012 annual report.

Usage

data(isleRoyal)

Format

The data are supplied as a matrix with years in the first column.

Source

Peterson R. O., Thomas N. J., Thurber J. M., Vucetich J. A. and Waite T. A. (1998) Population limitation and the wolves of Isle Royale. In: Biology and Conservation of Wild Canids (eds. D. Macdonald and C. Sillero-Zubiri). Oxford University Press, Oxford, pp. 281-292.

Vucetich, J. A. and R. O. Peterson. (2010) Ecological studies of wolves on Isle Royale. Annual Report 2009-10. School of Forest Resources and Environmental Science, Michigan Technological University, Houghton, Michigan USA 49931-1295

The source for the covariate data is the Western Regional Climate Center (http://www.wrcc.dri.edu) using their data for the NE Minnesota division. The website used was http://www.wrcc.dri.edu/cgi-bin/divplot1_form.pl?2103 and www.wrcc.dri.edu/spi/divplot1map.html.

Examples

str(isleRoyal)

---

Return a diagonal list matrix

Description

Creates a list diagonal matrix where the diagonal can be a combination of numbers and characters. Characters are names of parameters to be estimated.

Usage

ldiag(x, nrow = NA)

Arguments

x	A vector or list of single values
nrow	Rows in square matrix

Details

A diagonal list matrix is returned. The off-diagonals will be 0 and the diagonal will be x . x can be a combination of numbers and characters. If x is numeric, the diagonal will still be list type so that later the diagonal can be replace with characters. See examples.

Value

a square list matrix

Author(s)

Eli Holmes, NOAA, Seattle, USA.

Examples

ldiag(list(0, "b"))
ldiag("a", nrow=3)

# This works
a <- ldiag(1:3)
diag(a) <- "a"
diag(a) <- list("a", 0, 0)

# ldiag() is a convenience function to replace having to 
# write code like this
a <- matrix(list(0), 3, 3)
diag(a) <- list("a", 0, 0)

# diag() alone won't work because it cannot handle mixed number/char lists

# This turns the off-diagonals to character "0" 
a <- diag(1:3)
diag(a) <- "a"

# This would turn our matrix into a list only (not a matrix anymore)
a <- diag(1:3)
diag(a) <- list("a", 0, 0)

# This would return NA on the diagonal
a <- diag("a", 3)

---

Loggerhead Turtle Tracking Data

Description

Data used in MARSS vignettes in the MARSS-package. Tracking data from ARGOS tags on eight individual loggerhead turtles, 1997-2006.

Usage

data(loggerhead)
data(loggerheadNoisy)

Format

Data frames "loggerhead" and "loggerheadNoisy" contain the following columns:

turtle

Turtle name.

day

Day of the month (character).

month

Month number (character).

year

Year (character).

lon

Longitude of observation.

lat

Latitude of observation.

Details

Data frame "loggerhead" contains the original latitude and longitude data. Data frame "loggerheadNoisy" has noise added to the lat and lon data to represent data corrupted by errors.

Source

Gray's Reef National Marine Sanctuary (Georgia) and WhaleNet: http://whale.wheelock.edu/whalenet-stuff/stop_cover_archive.html

Examples

str(loggerhead)
str(loggerheadNoisy)

---

logLik method for MARSS MLE objects

Description

Returns a logLik class object with attributes nobs and df.

Usage

## S3 method for class 'marssMLE'
logLik(object, ...)

Arguments

object	A marssMLE object.
...	Other arguments. Not used.

Value

An object of class logLik.

Author(s)

Eli Holmes, NOAA, Seattle, USA.

See Also

MARSSkf()

Examples

dat <- t(harborSeal)
dat <- dat[c(2, 11, 12), ]
MLEobj <- MARSS(dat, model = list(Z = factor(c("WA", "OR", "OR"))))
logLik(MLEobj)

---

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)

---

Multivariate Dynamic Factor Analysis

Description

The Dynamic Factor Analysis model in MARSS is The argument form="marxss" in a MARSS() function call specifies a MAR-1 model with eXogenous variables model. This is a MARSS(1) model of the form:

$\mathbf{x}_{t} = \mathbf{x}_{t-1} + \mathbf{w}_t, \textrm{ where } \mathbf{W}_t \sim \textrm{MVN}(0,\mathbf{I})$

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

$\mathbf{X}_1 \sim \textrm{MVN}(\mathbf{x0}, 5\mathbf{I})$

Note, by default $\mathbf{x}_1$ is treated as a diffuse prior.

Passing in form="dfa" to MARSS() invokes a helper function to create that model and creates the $\mathbf{Z}$ matrix for the user. $\mathbf{Q}$ is by definition identity, $\mathbf{x}_0$ is zero and $\mathbf{V_0}$ is diagonal with large variance (5). $\mathbf{u}$ is zero, $\mathbf{a}$ is zero, and covariates only enter the $\mathbf{y}$ equation. Because $\mathbf{u}$ and $\mathbf{a}$ are 0, the data should have mean 0 (demeaned) otherwise one is likely to be creating a structurally inadequate model (i.e. the model implies that the data have mean = 0, yet data do not have mean = 0 ).

Arguments

Some arguments are common to all forms: "y" (data), "inits", "control", "method", "form", "fit", "silent", "fun.kf". See MARSS for information on these arguments.

In addition to these, form="dfa" has some special arguments that can be passed in:

• demean Logical. Default is TRUE, which means the data will be demeaned.

• z.score Logical. Default is TRUE, which means the data will be z-scored (demeaned and variance standardized to 1).

• covariates Covariates ($d$) for the $y$ equation. No missing values allowed and must be a matrix with the same number of time steps as the data. An unconstrained $D$ matrix will estimated.

The model argument of the MARSS() call is constrained in terms of what parameters can be changed and how they can be changed. See details below. An additional element, m, can be passed into the model argument that specifies the number of hidden state variables. It is not necessarily for the user to specify Z as the helper function will create a Z appropriate for a DFA model.

Details

The model argument is a list. The following details what list elements can be passed in:

• B "Identity". The standard (and default) DFA model has B="identity". However it can be "identity", "diagonal and equal", "diagonal and unequal" or a time-varying fixed or estimated diagonal matrix.

• U "Zero". Cannot be changed or passed in via model argument.

• Q "Identity". The standard (and default) DFA model has Q="identity". However, it can be "identity", "diagonal and equal", "diagonal and unequal" or a time-varying fixed or estimated diagonal matrix.

• Z Can be passed in as a (list) matrix if the user does not want a default DFA Z matrix. There are many equivalent ways to construct a DFA Z matrix. The default is Zuur et al.'s form (see User Guide).

• A Default="zero". Can be "unequal", "zero" or a matrix.

• R Default="diagonal and equal". Can be set to "identity", "zero", "unconstrained", "diagonal and unequal", "diagonal and equal", "equalvarcov", or a (list) matrix to specify general forms.

• x0 Default="zero". Can be "unconstrained", "unequal", "zero", or a (list) matrix.

• V0 Default=diagonal matrix with 5 on the diagonal. Can be "identity", "zero", or a matrix.

• tinitx Default=0. Can be 0 or 1. Tells MARSS whether x0 is at t=0 or t=1.

• m Default=1. Can be 1 to n (the number of y time-series). Must be integer.

See the User Guide chapter on Dynamic Factor Analysis for examples of of using form="dfa".

Value

A object of class marssMLE. See print() for a discussion of the various output available for marssMLE objects (coefficients, residuals, Kalman filter and smoother output, imputed values for missing data, etc.). See MARSSsimulate() for simulating from marssMLE objects. MARSSboot() for bootstrapping, MARSSaic() for calculation of various AIC related model selection metrics, and MARSSparamCIs() for calculation of confidence intervals and bias.

Usage

MARSS(y, inits = NULL, model = NULL, miss.value = as.numeric(NA), method = "kem", form = "dfa", fit = TRUE, silent = FALSE, control = NULL, fun.kf = "MARSSkfas", demean = TRUE, z.score = TRUE)

Author(s)

Eli Holmes, 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.

See Also

MARSS(), MARSS.marxss()

Examples

## Not run: 
dat <- t(harborSealWA[,-1])
# DFA with 3 states; used BFGS because it fits much faster for this model
fit <- MARSS(dat, model = list(m=3), form="dfa", method="BFGS")

# See the Dynamic Factor Analysis chapter in the User Guide
RShowDoc("UserGuide", package = "MARSS")

## End(Not run)

---

Multivariate AR-1 State-space Model

Description

The form of MARSS models for users is "marxss", the MARSS models with inputs. See MARSS.marxss. In the internal algorithms (e.g. MARSSkem), the "marss" form is used and the $\mathbf{D}\mathbf{d}_t$ are incorporated into the $\mathbf{a}_t$ matrix and $\mathbf{C}\mathbf{c}_t$ are incorporated into the $\mathbf{u}_t$. The $\mathbf{a}$ and $\mathbf{u}$ matrices then become time-varying if the model includes $\mathbf{d}_t$ and $\mathbf{c}_t$.

This is a MARSS(1) model of the marss form:

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

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

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

Note, by default $\mathbf{V0}$ is a matrix of all zeros and thus $\mathbf{x}_1$ or $\mathbf{x}_0$ is treated as an estimated parameter not a diffuse prior. To remove clutter, the rest of the parameters are shown as time-constant (no $t$ subscript) but all parameters can be time-varying.

Note, "marss" is a model form. A model form is defined by a collection of form functions discussed in marssMODEL. These functions are not exported to the user, but are called by MARSS() using the argument form. These internal functions convert the users model list into the vec form of a MARSS model and do extensive error-checking.

Details

See the help page for the MARSS.marxss form for details.

Value

A object of class marssMLE.

Usage

MARSS(y, inits = NULL, model = NULL, miss.value = as.numeric(NA), method = "kem", form = "marxss", fit = TRUE, silent = FALSE, control = NULL, fun.kf = "MARSSkfas", ...)

Author(s)

Eli Holmes, NOAA, Seattle, USA.

See Also

marssMODEL, MARSS.marxss()

Examples

## Not run: 
# See the MARSS man page for examples
?MARSS

# and the Quick Examples chapter in the User Guide
RShowDoc("UserGuide", package = "MARSS")

## End(Not run)

---

Multivariate AR-1 State-space Model with Inputs

Description

The argument form="marxss" in a MARSS() function call specifies a MAR-1 model with eXogenous variables model. This is a MARSS(1) model of the form:

$\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})$

Note, by default $\mathbf{V0}$ is a matrix of all zeros and thus $\mathbf{x}_1$ or $\mathbf{x}_0$ is treated as an estimated parameter not a diffuse prior.

Note, "marxss" is a model form. A model form is defined by a collection of form functions discussed in marssMODEL. These functions are not exported to the user, but are called by MARSS() using the argument form.

Details

The allowed arguments when form="marxss" are 1) the arguments common to all forms: "y" (data), "inits", "control", "method", "form", "fit", "silent", "fun.kf" (see MARSS for information on these arguments) and 2) the argument "model" which is a list describing the MARXSS model (the model list is described below). See the Quick Start Guide guide or the User Guide for examples.

The argument model must be a list. The elements in the list specify the structure for the $\mathbf{B}$, $\mathbf{u}$, $\mathbf{C}$, $\mathbf{c}$, $\mathbf{Q}$, $\mathbf{Z}$, $\mathbf{a}$, $\mathbf{D}$, $\mathbf{d}$, $\mathbf{R}$, $\mathbf{x}_0$, and $\mathbf{V}_0$ in the MARXSS model (above). The list elements can have the following values:

Z

Default="identity". A text string, "identity","unconstrained", "diagonal and unequal", "diagonal and equal", "equalvarcov", or "onestate", or a length n vector of factors specifying which of the m hidden state time series correspond to which of the n observation time series. May be specified as a n x m list matrix for general specification of both fixed and shared elements within the matrix. May also be specified as a numeric n x m matrix to use a custom fixed $\mathbf{Z}$. "onestate" gives a n x 1 matrix of 1s. "identity","unconstrained", "diagonal and unequal", "diagonal and equal", and "equalvarcov" all specify n x n matrices.

B

Default="identity". A text string, "identity", "unconstrained", "diagonal and unequal", "diagonal and equal", "equalvarcov", "zero". Can also be specified as a list matrix for general specification of both fixed and shared elements within the matrix. May also be specified as a numeric m x m matrix to use custom fixed $\mathbf{B}$, but in this case all the eigenvalues of $\mathbf{B}$ must fall in the unit circle.

U, x0

Default="unconstrained". A text string, "unconstrained", "equal", "unequal" or "zero". May be specified as a m x 1 list matrix for general specification of both fixed and shared elements within the matrix. May also be specified as a numeric m x 1 matrix to use a custom fixed $\mathbf{u}$ or $\mathbf{x}_0$. Notice that U is capitalized in the model argument and output lists.

A

Default="scaling". A text string, "scaling","unconstrained", "equal", "unequal" or "zero". May be specified as a n x 1 list matrix for general specification of both fixed and shared elements within the matrix. May also be specified as a numeric n x 1 matrix to use a custom fixed $\mathbf{a}$. Care must be taken when specifying A so that the model is not under-constrained and unsolvable model. The default, "scaling", only applies to $\mathbf{Z}$ matrices that are design matrices (only 1s and 0s and all rows sum to 1). When a column in $\mathbf{Z}$ has multiple 1s, the first row in the $\mathbf{a}$ matrix associated with those $\mathbf{Z}$ rows is 0 and the other associated $\mathbf{a}$ rows have an estimated value. This is used to treat $\mathbf{a}$ as an intercept where one intercept for each $\mathbf{x}$ (hidden state) is fixed at 0 and any other intercepts associated with that $\mathbf{x}$ have an estimated intercept. This ensures a solvable model when $\mathbf{Z}$ is a design matrix. Note in the model argument and output, A is capitalized.

Q

Default="diagonal and unequal". A text string, "identity", "unconstrained", "diagonal and unequal", "diagonal and equal", "equalvarcov", "zero". May be specified as a list matrix for general specification of both fixed and shared elements within the matrix. May also be specified as a numeric g x g matrix to use a custom fixed matrix. Default value of g is m, so $\mathbf{Q}$ is a m x m matrix. g is the number of columns in $\mathbf{G}$ (below).

R

Default="diagonal and equal". A text string, "identity", "unconstrained", "diagonal and unequal", "diagonal and equal", "equalvarcov", "zero". May be specified as a list matrix for general specification of both fixed and shared elements within the matrix. May also be specified as a numeric h x h matrix to use a custom fixed matrix. Default value of h is n, so $\mathbf{R}$ is a n x n matrix. h is the num of columns in $\mathbf{H}$ (below).

V0

Default="zero". A text string, "identity", "unconstrained", "diagonal and unequal", "diagonal and equal", "equalvarcov", "zero". May be specified as a list matrix for general specification of both fixed and shared elements within the matrix. May also be specified as a numeric m x m matrix to use a custom fixed matrix.

D and C

Default="zero". A text string, "identity", "unconstrained", "diagonal and unequal", "diagonal and equal", "equalvarcov", "zero". Can be specified as a list matrix for general specification of both fixed and shared elements within the matrix. May also be specified as a numeric matrix to use custom fixed values. Must have n rows ($\mathbf{D}$) or m rows ($\mathbf{C}$).

d and c

Default="zero". Numeric matrix. No missing values allowed. Must have 1 column or the same number of columns as the data, $\mathbf{y}$. The numbers of rows in $\mathbf{d}$ must be the same as number of columns in $\mathbf{D}$; similarly for $\mathbf{c}$ and $\mathbf{C}$.

G and H

Default="identity". A text string, "identity". Can be specified as a numeric matrix or array for time-varying cases. Must have m rows and g columns ($\mathbf{G}$) or n rows and h columns ($\mathbf{H}$). g is the dim of $\mathbf{Q}$ and h is the dim of $\mathbf{R}$.

tinitx

Default=0. Whether the initial state is specified at t=0 (default) or t=1.

All parameters except $\mathbf{x}_0$ and $\mathbf{V}_0$ may be time-varying. If time-varying, then text shortcuts cannot be used. Enter as an array with the 3rd dimension being time. Time dimension must be 1 or equal to the number of time-steps in the data. See Quick Start guide (RShowDoc("Quick_Start",package="MARSS")) or the User Guide (RShowDoc("UserGuide",package="MARSS")) for examples.Valid model structures for method="BFGS" are the same as for method="kem". See MARSSoptim() for the allowed options for this method.

The default estimation method, method="kem", is the EM algorithm described in the MARSS User Guide. The default settings for the control and inits arguments are set via MARSS:::alldefaults$kem in MARSSsettings.R. The defaults for the model argument are set in MARSS_marxss.R For this method, they are:

• inits = list(B=1, U=0, Q=0.05, Z=1, A=0, R=0.05, x0=-99, V0=0.05, G=0, H=0, L=0, C=0, D=0, c=0, d=0)

• model = list(Z="identity", A="scaling", R="diagonal and equal", B="identity", U="unconstrained", Q="diagonal and unequal", x0="unconstrained", V0="zero", C="zero",D="zero",c=matrix(0,0,1), d=matrix(0,0,1), tinitx=0, diffuse=FALSE)

• control=list(minit=15, maxit=500, abstol=0.001, trace=0, sparse=FALSE, safe=FALSE, allow.degen=TRUE, min.degen.iter=50, degen.lim=1.0e-04, min.iter.conv.test=15, conv.test.deltaT=9, conv.test.slope.tol= 0.5, demean.states=FALSE) You can read about these in MARSS(). If you want to speed up your fits, you can turn off most of the model checking using trace=-1.

• fun.kf = "MARSSkfas"; This sets the Kalman filter function to use. MARSSkfas() is generally more stable as it uses Durban & Koopman's algorithm. But it may dramatically slow down when the data set is large (more than 10 rows of data). Try the classic Kalman filter algorithm to see if it runs faster by setting fun.kf="MARSSkfss". You can read about the two algorithms in MARSSkf.

For method="BFGS", type MARSS:::alldefaults$BFGS to see the defaults.

Value

A object of class marssMLE. See print.marssMLE for a discussion of the various output available for marssMLE objects (coefficients, residuals, Kalman filter and smoother output, imputed values for missing data, etc.). See MARSSsimulate for simulating from marssMLE objects. MARSSboot for bootstrapping, MARSSaic for calculation of various AIC related model selection metrics, and MARSSparamCIs for calculation of confidence intervals and bias. See plot.marssMLE for some default plots of a model fit.

Usage

MARSS(y, inits = NULL, model = NULL, miss.value = as.numeric(NA), method = "kem", form = "marxss", fit = TRUE, silent = FALSE, control = NULL, fun.kf = "MARSSkfas", ...)

Author(s)

Eli Holmes, NOAA, Seattle, USA.

See Also

marssMODEL, MARSS.dfa()

Examples

## Not run: 
#See the MARSS man page for examples
?MARSS

#and the Quick Examples chapter in the User Guide
RShowDoc("UserGuide",package="MARSS")

## End(Not run)

---

Vectorized Multivariate AR-1 State-space Model

Description

The EM algorithm (MARSSkem) in the MARSS package works by converting the more familiar MARSS model in matrix form into the vectorized form which allows general linear constraints (Holmes 2012). The vectorized form is:

$\mathbf{x}(t) = (\mathbf{x}(t-1)^\top \otimes \mathbf{I}_m)(\mathbf{f}_b(t)+\mathbf{D}_b(t)\beta) + (\mathbf{f}_u(t)+\mathbf{D}_u(t)\upsilon) + \mathbf{w}(t), \textrm{ where } \mathbf{W}(t) \sim \textrm{MVN}(0,\textbf{Q}(t))$

$\mathbf{y}(t) = (\mathbf{x}(t)^\top \otimes \mathbf{I}_n)(\mathbf{f}_z(t)+\mathbf{D}_z(t)\zeta) + (\mathbf{f}_a(t)+\mathbf{D}_a(t)\alpha) + \mathbf{v}(t), \textrm{ where } \mathbf{V}(t) \sim \textrm{MVN}(0,\textbf{R}(t))$

$\mathbf{x}(1) \sim \textrm{MVN}(x0, V0) \textrm{ or } \mathbf{x}(0) \sim \textrm{MVN}(x0, V0)$

where $\beta$, $\upsilon$, $\zeta$, and $\alpha$ are column vectors of estimated values, the $\mathbf{f}$ are column vectors of inputs (fixed values), and the $\mathbf{D}$ are perturbation matrices that align the estimated values into the right rows. The $\mathbf{f}$ and $\mathbf{D}$ are potentially time-varying. $\otimes$ means kronecker product and $\mathbf{I}_p$ is a p x p identity matrix.

Normally the user will specify their model in "marxss" form, perhaps with text short-cuts. The "marxss" form is then converted to "marss" form using the conversion function marxss_to_marss(). In "marss" form, the D, d, C, and c information is put in A and U respectively. If there are inputs (d and c), then this will make A and U time-varying. This is unfortunate, because this slows down the EM algorithm considerably due to the unfortunate decision (early on) to store time-varying parameters as 3-dimensional. The functions for the "marss" form (in the file MARSS_marss.R) convert the "marss" form model into vectorized form and prepares the f (fixed) and D (free) matrices that are at the heart of the model specification.

Note, "marss" is a model form. A model form is defined by a collection of form functions discussed in marssMODEL. These functions are not exported to the user, but are called by MARSS() using the argument form. These internal functions convert the users model list into the vectorized form of a MARSS model and do extensive error-checking. "marxss" is also a model form and these models are also stored in vectorized form (See examples below).

Details

See Holmes (2012) for a discussion of MARSS models in vectorized form.

Author(s)

Eli Holmes, NOAA, Seattle, USA.

References

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]

See Also

marssMODEL, MARSS.marss(), MARSS.marxss()

Examples

dat <- t(harborSealWA)
dat <- dat[2:4, ] 
MLEobj <- MARSS(dat)

# free (D) and fixed (f) matrices
names(MLEobj$model$free)
names(MLEobj$model$fixed)
# In marss form, the D, C, d, and c matrices are found in A and U
# If there are inputs, this makes U time-varying
names(MLEobj$marss$free)
names(MLEobj$marss$fixed)

# par is in marss form so does not have values for D, C, d, or c
names(MLEobj$par)
# if you need the par in marxss form, you can use print
tmp <- print(MLEobj, what="par", form="marxss", silent=TRUE)
names(tmp)

---

AIC for MARSS Models

Description

Calculates AIC, AICc, a parametric bootstrap AIC (AICbp) and a non-parametric bootstrap AIC (AICbb). If you simply want the AIC value for a marssMLE object, you can use AIC(fit).

Usage

MARSSaic(MLEobj, output = c("AIC", "AICc"), 
    Options = list(nboot = 1000, return.logL.star = FALSE, 
    silent = FALSE))

Arguments

MLEobj	An object of class marssMLE. This object must have a $par element containing MLE parameter estimates from e.g. MARSSkem().
output	A vector containing one or more of the following: "AIC", "AICc", "AICbp", "AICbb", "AICi", "boot.params". See Details.
Options	A list containing: nboot Number of bootstraps (positive integer) return.logL.star Return the log-likelihoods for each bootstrap? (T/F) silent Suppress printing of the progress bar during AIC bootstraps? (T/F)

Details

When sample size is small, Akaike's Information Criterion (AIC) under-penalizes more complex models. The most commonly used small sample size corrector is AICc, which uses a penalty term of $K n/(n-K-1)$, where $K$ is the number of estimated parameters. However, for time series models, AICc still under-penalizes complex models; this is especially true for MARSS models.

Two small-sample estimators specific for MARSS models have been developed. Cavanaugh and Shumway (1997) developed a variant of bootstrapped AIC using Stoffer and Wall's (1991) bootstrap algorithm ("AICbb"). Holmes and Ward (2010) developed a variant on AICb ("AICbp") using a parametric bootstrap. The parametric bootstrap permits AICb calculation when there are missing values in the data, which Cavanaugh and Shumway's algorithm does not allow. More recently, Bengtsson and Cavanaugh (2006) developed another small-sample AIC estimator, AICi, based on fitting candidate models to multivariate white noise.

When the output argument passed in includes both "AICbp" and "boot.params", the bootstrapped parameters from "AICbp" will be added to MLEobj.

Value

Returns the marssMLE object that was passed in with additional AIC components added on top as specified in the 'output' argument.

Author(s)

Eli Holmes, NOAA, Seattle, USA.

References

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.

Bengtsson, T., and J. E. Cavanaugh. 2006. An improved Akaike information criterion for state-space model selection. Computational Statistics & Data Analysis 50:2635-2654.

Cavanaugh, J. E., and R. H. Shumway. 1997. A bootstrap variant of AIC for state-space model selection. Statistica Sinica 7:473-496.

See Also

MARSSboot()

Examples

dat <- t(harborSealWA)
dat <- dat[2:3, ]
kem <- MARSS(dat, model = list(
  Z = matrix(1, 2, 1),
  R = "diagonal and equal"
))
kemAIC <- MARSSaic(kem, output = c("AIC", "AICc"))

---

Bootstrap MARSS Parameter Estimates

Description

Creates bootstrap parameter estimates and simulated (or bootstrapped) data (if appropriate). This is a base function in the MARSS-package.

Usage

MARSSboot(MLEobj, nboot = 1000, 
    output = "parameters", sim = "parametric", 
    param.gen = "MLE", control = NULL, silent = FALSE)

Arguments

MLEobj	An object of class marssMLE. Must have a $par element containing MLE parameter estimates.
nboot	Number of bootstraps to perform.
output	Output to be returned: "data", "parameters" or "all".
sim	Type of bootstrap: "parametric" or "innovations". See Details.
param.gen	Parameter generation method: "hessian" or "MLE".
control	The options in MLEobj$control are used by default. If supplied here, must contain all of the following: max.iter Maximum number of EM iterations. tol Optional tolerance for log-likelihood change. If log-likelihood decreases less than this amount relative to the previous iteration, the EM algorithm exits. allow.degen Whether to try setting $\mathbf{Q}$ or $\mathbf{R}$ elements to zero if they appear to be going to zero.
silent	Suppresses printing of progress bar.

Details

Approximate confidence intervals (CIs) on the model parameters can be calculated by the observed Fisher Information matrix (the Hessian of the negative log-likelihood function). The Hessian CIs (param.gen="hessian") are based on the asymptotic normality of ML estimates under a large-sample approximation. CIs that are not based on asymptotic theory can be calculated using parametric and non-parametric bootstrapping (param.gen="MLE"). In this case, parameter estimates are generated by the ML estimates from each bootstrapped data set. The MLE method (kem or BFGS) is determined by MLEobj$method.

Stoffer and Wall (1991) present an algorithm for generating CIs via a non-parametric bootstrap for state-space models (sim = "innovations"). The basic idea is that the Kalman filter can be used to generate estimates of the residuals of the model fit. These residuals are then standardized and resampled and used to generate bootstrapped data using the MARSS model and its maximum-likelihood parameter estimates. One of the limitations of the Stoffer and Wall algorithm is that it cannot be used when there are missing data, unless all data at time $t$ are missing. An alternative approach is a parametric bootstrap (sim = "parametric"), in which the ML parameter estimates are used to produce bootstrapped data directly from the state-space model.

Value

A list with the following components:

boot.params	Matrix (number of params x nboot) of parameter estimates from the bootstrap.
boot.data	Array (n x t x nboot) of simulated (or bootstrapped) data (if requested and appropriate).
marss	The marssMODEL object (form="marss") that was passed in via MLEobj$marss.
nboot	Number of bootstraps performed.
output	Type of output returned.
sim	Type of bootstrap.
param.gen	Parameter generation method: "hessian" or "KalmanEM".

Author(s)

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

References

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.

Stoffer, D. S., and K. D. Wall. 1991. Bootstrapping state-space models: Gaussian maximum likelihood estimation and the Kalman filter. Journal of the American Statistical Association 86:1024-1033.

Cavanaugh, J. E., and R. H. Shumway. 1997. A bootstrap variant of AIC for state-space model selection. Statistica Sinica 7:473-496.

See Also

marssMLE, marssMODEL, MARSSaic(), MARSShessian(), MARSSFisherI()

Examples

# nboot is set low in these examples in order to run quickly
# normally nboot would be >1000 at least
dat <- t(kestrel)
dat <- dat[2:3, ]
# maxit set low to speed up the example
kem <- MARSS(dat,
  model = list(U = "equal", Q = diag(.01, 2)),
  control = list(maxit = 50)
)
# bootstrap parameters from a Hessian matrix
hess.list <- MARSSboot(kem, param.gen = "hessian", nboot = 4)

# from resampling the innovations (no missing values allowed)
boot.innov.list <- MARSSboot(kem, output = "all", sim = "innovations", nboot = 4)

# bootstrapped parameter estimates
hess.list$boot.params

---

MARSScv is a wrapper for MARSS that re-fits the model with cross validated data.

Description

MARSScv is a wrapper for MARSS that re-fits the model with cross validated data.

Usage

MARSScv(
  y,
  model = NULL,
  inits = NULL,
  method = "kem",
  form = "marxss",
  silent = FALSE,
  control = NULL,
  fun.kf = c("MARSSkfas", "MARSSkfss"),
  fold_ids = NULL,
  future_cv = FALSE,
  n_future_cv = floor(ncol(y)/3),
  interval = "confidence",
  ...
)

Arguments

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.
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")
inits	A list with the same form as the list output by coef(fit) that specifies initial values for the parameters. See also MARSS.marxss().
method	Estimation method. MARSS provides an EM algorithm (method="kem") (see MARSSkem()) and the BFGS algorithm (method="BFGS") (see MARSSoptim()). Fast TMB fitting provided by the companion package marssTMB.
form	The equation form used in the MARSS() call. The default is "marxss". See MARSS.marxss() or MARSS.dfa()
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.
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.
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).
fold_ids	A n x T matrix of integers, with values assigned by the user to folds. If not included, data are randomly assigned to one of 10 folds
future_cv	Whether or not to use future cross validation (defaults to FALSE), where data up to time T-1 are used to predict data at time T. Data are held out by time slices, and the fold_ids argument is ignored.
n_future_cv	Number of time slices to hold out for future cross validation. Defaults to floor(n_future_cv/3). Predictions are made for the last n_future_cv time steps
interval	uncertainty interval for prediction. Can be one of "confidence" or "prediction", and defaults to "confidence"
...	not used

Value

A list object, containing cv_pred (a matrix of predictions), cv_se (a matrix of SEs), fold_ids (a matrix of fold ids used as data), and df (a dataframe containing the original data, predictions, SEs, and folds)

Examples

dat <- t(harborSealWA)
dat <- dat[2:4, ] # remove the year row
# fit a model with 1 hidden state and 3 observation time series
# cross validation here is random, 10 folds
fit <- MARSScv(dat, model = list(
  Z = matrix(1, 3, 1),
  R = "diagonal and equal"
))

# second, demonstrate passing in pre-specified folds
fold_ids <- matrix(
  sample(1:5, size = nrow(dat) * ncol(dat), replace = TRUE),
  nrow(dat), ncol(dat)
)
fit <- MARSScv(dat, model = list(
  Z = matrix(1, 3, 1),
  R = "diagonal and equal"
), fold_ids = fold_ids)

# third, illustrate future cross validation
fit <- MARSScv(dat, model = list(
  Z = matrix(1, 3, 1),
  R = "diagonal and equal"
), future_cv = TRUE, n_future_cv = 5)

---