MARSSresiduals: MARSS Residuals

In MARSS: Multivariate Autoregressive State-Space Modeling

MARSS Residuals

Description

The normal residuals function is residuals(). MARSSresiduals() returns residuals as a list of matrices while residuals() returns the same information in a data frame. This function calculates the residuals, residuals variance, and standardized residuals for the one-step-ahead (conditioned on data up to t-1), the smoothed (conditioned on all the data), and contemporaneous (conditioned on data up to t) residuals.

Usage

MARSSresiduals(object, ..., type = c("tT", "tt1", "tt"), 
    normalize = FALSE, silent = FALSE, 
    fun.kf = c("MARSSkfas", "MARSSkfss"))

Arguments

object	An object of class marssMLE.
...	Additional arguments to be passed to the residuals functions. For type="tT", Harvey=TRUE can be passed into to use the Harvey et al (1998) algorithm.
type	"tT" for smoothed residuals conditioned on all the data t=1 to T, aka smoothation residuals. "tt1" for one-step-ahead residuals, aka innovations residuals. "tt" for contemporaneous residuals.
normalize	TRUE/FALSE See details.
silent	If TRUE, do not print inversion warnings.
fun.kf	Kalman filter function to use. Can be ignored.

Details

For smoothed residuals, see MARSSresiduals.tT().

For one-step-ahead residuals, see MARSSresiduals.tt1().

For contemporaneous residuals, see MARSSresiduals.tt().

Standardized residuals

Standardized residuals have been adjusted by the variance of the residuals at time t such that the variance of the residuals at time t equals 1. Given the normality assumption, this means that one typically sees +/- 2 confidence interval lines on standardized residuals plots.

std.residuals are Cholesky standardized residuals. These are the residuals multiplied by the inverse of the lower triangle of the Cholesky decomposition of the variance matrix of the residuals:

\hat{\Sigma}_t^{-1/2} \hat{\mathbf{v}}_t.

These residuals are uncorrelated with each other, although they are not necessarily temporally uncorrelated (innovations residuals are temporally uncorrelated).

The interpretation of the Cholesky standardized residuals is not straight-forward when the \mathbf{Q} and \mathbf{R} variance-covariance matrices are non-diagonal. The residuals which were generated by a non-diagonal variance-covariance matrices are transformed into orthogonal residuals in \textrm{MVN}(0,\mathbf{I}) space. For example, if v is 2x2 correlated errors with variance-covariance matrix \mathbf{R}. The transformed residuals (from this function) for the i-th row of v is a combination of the row 1 effect and the row 1 effect plus the row 2 effect. So in this case, row 2 of the transformed residuals would not be regarded as solely the row 2 residual but rather how different row 2 is from row 1, relative to expected. If the errors are highly correlated, then the Cholesky standardized residuals can look rather non-intuitive.

mar.residuals are the marginal standardized residuals. These are the residuals multiplied by the inverse of the diagonal matrix formed from the square-root of the diagonal of the variance matrix of the residuals:

\textrm{dg}(\hat{\Sigma}_t)^{-1/2} \hat{\mathbf{v}}_t,

where dg(A) is the square matrix formed from the diagonal of A, aka diag(diag(A)). These residuals will be correlated if the variance matrix is non-diagonal.

The Block Cholesky standardized residuals are like the Cholesky standardized residuals except that the full variance-covariance matrix is not used, only the variance-covariance matrix for the model or state residuals (respectively) is used for standardization. For the model residuals, the Block Cholesky standardized residuals will be the same as the Cholesky standardized residuals because the upper triangle of the lower triangle of the Cholesky decomposition (which is what we standardize by) is all zero. For type="tt1" and type="tt", the Block Cholesky standardized state residuals will be the same as the Cholesky standardized state residuals because in the former, the model and state residuals are uncorrelated and in the latter, the state residuals do not exist. For type="tT", the model and state residuals are correlated and the Block Cholesky standardized residuals will be different than the Cholesky standardized residuals.

Normalized residuals

If normalize=FALSE, the unconditional variance of \mathbf{V}_t and \mathbf{W}_t are \mathbf{R} and \mathbf{Q} and the model is assumed to be written as

\mathbf{y}_t = \mathbf{Z} \mathbf{x}_t + \mathbf{a} + \mathbf{v}_t

\mathbf{x}_t = \mathbf{B} \mathbf{x}_{t-1} + \mathbf{u} + \mathbf{w}_t

If normalize=TRUE, the model is assumed to be written as

\mathbf{y}_t = \mathbf{Z} \mathbf{x}_t + \mathbf{a} + \mathbf{H}\mathbf{v}_t

\mathbf{x}_t = \mathbf{B} \mathbf{x}_{t-1} + \mathbf{u} + \mathbf{G}\mathbf{w}_t

with the variance of \mathbf{V}_t and \mathbf{W}_t equal to \mathbf{I} (identity).

Missing or left-out data

See the discussion of residuals for missing and left-out data in MARSSresiduals.tT().

Value

A list of the following components

model.residuals	The model residuals (data minus model predicted values) as a n x T matrix.
state.residuals	The state residuals. This is the state residual for the transition from t=t to t+1 thus the last time step will be NA (since T+1 is past the data). State residuals do not exist for the type="tt" case (since this would required the expected value of \mathbf{X}_t conditioned on data to t+1).
residuals	The residuals as a (n+m) x T matrix with model.residuals on top and state.residuals below.
var.residuals	The variance of the model residuals and state residuals as a (n+m) x (n+m) x T matrix with the model residuals variance in rows/columns 1 to n and state residuals variances in rows/columns n+1 to n+m. The last time step will be all NA since the state residual is for t=t to t+1.
std.residuals	The Cholesky standardized residuals as a (n+m) x T matrix. This is residuals multiplied by the inverse of the lower triangle of the Cholesky decomposition of var.residuals.
mar.residuals	The marginal standardized residuals as a (n+m) x T matrix. This is residuals multiplied by the inverse of the diagonal matrix formed by the square-root of the diagonal of var.residuals.
bchol.residuals	The Block Cholesky standardized residuals as a (n+m) x T matrix. This is model.residuals multiplied by the inverse of the lower triangle of the Cholesky decomposition of var.residuals[1:n,1:n,] and state.residuals multiplied by the inverse of the lower triangle of the Cholesky decomposition of var.residuals[(n+1):(n+m),(n+1):(n+m),].
E.obs.residuals	The expected value of the model residuals conditioned on the observed data. Returned as a n x T matrix. For observed data, this will be the observed model residuals. For unobserved data, this will be 0 if \mathbf{R} is diagonal but non-zero if \mathbf{R} is non-diagonal. See MARSSresiduals.tT().
var.obs.residuals	The variance of the model residuals conditioned on the observed data. Returned as a n x n x T matrix. For observed data, this will be 0. See MARSSresiduals.tT().
msg	Any warning messages. This will be printed unless Object$control$trace = -1 (suppress all error messages).

Author(s)

Eli Holmes, NOAA, Seattle, USA.

References

Holmes, E. E. 2014. Computation of standardized residuals for (MARSS) models. Technical Report. arXiv:1411.0045.

See also the discussion and references in MARSSresiduals.tT(), MARSSresiduals.tt1() and MARSSresiduals.tt().

See Also

residuals.marssMLE(), MARSSresiduals.tT(), MARSSresiduals.tt1(), plot.marssMLE()

Examples

  dat <- t(harborSeal)
  dat <- dat[c(2,11),]
  fit <- MARSS(dat)
  
  #state smoothed residuals
  state.resids1 <- MARSSresiduals(fit, type="tT")$state.residuals
  #this is the same as
  states <- fit$states
  Q <- coef(fit, type="matrix")$Q
  state.resids2 <- states[,2:30]-states[,1:29]-matrix(coef(fit,type="matrix")$U,2,29)
  #compare the two
  cbind(t(state.resids1[,-30]), t(state.resids2))

  #normalize to variance of 1
  state.resids1 <- MARSSresiduals(fit, type="tT", normalize=TRUE)$state.residuals
  state.resids2 <- (solve(t(chol(Q))) %*% state.resids2)
  cbind(t(state.resids1[,-30]), t(state.resids2))

  #one-step-ahead standardized residuals
  MARSSresiduals(fit, type="tt1")$std.residuals

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