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.
MARSSresiduals(object, ..., type=c("tT", "tt1", "tt"),
normalize=FALSE, silent=FALSE,
fun.kf=c("MARSSkfas", "MARSSkfss"))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.
"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.
TRUE/FALSE See details.
If TRUE, do not print inversion warnings.
Kalman filter function to use. Can be ignored.
A list of the following components
The model residuals (data minus model predicted values) as a n x T matrix.
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 \(X_t\) conditioned on data to \(t+1\)).
The residuals as a (n+m) x T matrix with model.residuals on top and state.residuals below.
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\).
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.
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.
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),].
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().
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().
Any warning messages. This will be printed unless Object$control$trace = -1 (suppress all error messages).
For smoothed residuals, see MARSSresiduals.tT().
For one-step-ahead residuals, see MARSSresiduals.tt1().
For contemporaneous residuals, see MARSSresiduals.tt().
Standardized residuals
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.
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 MVN(0,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
$$\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().
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.
residuals.marssMLE, MARSSresiduals.tT, MARSSresiduals.tt1, plot.marssMLE
# NOT RUN {
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
# }
Run the code above in your browser using DataLab