MARSS Residuals
MARSSresiduals.RdThe 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.
Arguments
- object
An object of class
marssMLE.- ...
Additional arguments to be passed to the residuals functions. For type="tT",
Harvey=TRUEcan 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.
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.residualson top andstate.residualsbelow.- 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
residualsmultiplied by the inverse of the lower triangle of the Cholesky decomposition ofvar.residuals.- mar.residuals
The marginal standardized residuals as a (n+m) x T matrix. This is
residualsmultiplied by the inverse of the diagonal matrix formed by the square-root of the diagonal ofvar.residuals.- bchol.residuals
The Block Cholesky standardized residuals as a (n+m) x T matrix. This is
model.residualsmultiplied by the inverse of the lower triangle of the Cholesky decomposition ofvar.residuals[1:n,1:n,]andstate.residualsmultiplied by the inverse of the lower triangle of the Cholesky decomposition ofvar.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).
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().
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().
Examples
dat <- t(harborSeal)
dat <- dat[c(2,11),]
fit <- MARSS(dat)
#> Success! abstol and log-log tests passed at 26 iterations.
#> Alert: conv.test.slope.tol is 0.5.
#> Test with smaller values (<0.1) to ensure convergence.
#>
#> MARSS fit is
#> Estimation method: kem
#> Convergence test: conv.test.slope.tol = 0.5, abstol = 0.001
#> Estimation converged in 26 iterations.
#> Log-likelihood: 11.74016
#> AIC: -9.480311 AICc: -6.3692
#>
#> Estimate
#> R.diag 0.0115
#> U.X.CoastalEstuaries 0.0613
#> U.X.OR.NorthCoast 0.0510
#> Q.(X.CoastalEstuaries,X.CoastalEstuaries) 0.0147
#> Q.(X.OR.NorthCoast,X.OR.NorthCoast) 0.0122
#> x0.X.CoastalEstuaries 7.3823
#> x0.X.OR.NorthCoast 6.2707
#> Initial states (x0) defined at t=0
#>
#> Standard errors have not been calculated.
#> Use MARSSparamCIs to compute CIs and bias estimates.
#>
#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))
#> X.CoastalEstuaries X.OR.NorthCoast X.CoastalEstuaries X.OR.NorthCoast
#> [1,] 1.127322e-02 1.305044e-03 1.127322e-02 1.305044e-03
#> [2,] 7.963885e-02 1.305044e-03 7.963885e-02 1.305044e-03
#> [3,] 1.003050e-01 5.894950e-03 1.003050e-01 5.894950e-03
#> [4,] 5.842127e-02 5.894950e-03 5.842127e-02 5.894950e-03
#> [5,] 5.842127e-02 5.894950e-03 5.842127e-02 5.894950e-03
#> [6,] 1.843052e-01 5.894950e-03 1.843052e-01 5.894950e-03
#> [7,] 7.296579e-02 5.894950e-03 7.296579e-02 5.894950e-03
#> [8,] -7.717274e-02 8.380770e-02 -7.717274e-02 8.380770e-02
#> [9,] -3.934571e-02 2.025237e-02 -3.934571e-02 2.025237e-02
#> [10,] 8.981293e-02 2.169964e-02 8.981293e-02 2.169964e-02
#> [11,] -4.456791e-02 -5.700953e-03 -4.456791e-02 -5.700953e-03
#> [12,] 1.276304e-01 1.399795e-01 1.276304e-01 1.399795e-01
#> [13,] -4.420982e-02 7.507078e-02 -4.420982e-02 7.507078e-02
#> [14,] -5.118714e-02 7.306581e-02 -5.118714e-02 7.306581e-02
#> [15,] 1.668539e-02 -7.661036e-03 1.668539e-02 -7.661036e-03
#> [16,] 1.668539e-02 -2.730277e-02 1.668539e-02 -2.730277e-02
#> [17,] -9.075282e-02 1.000082e-01 -9.075282e-02 1.000082e-01
#> [18,] -9.256248e-02 -3.706648e-02 -9.256248e-02 -3.706648e-02
#> [19,] -1.985319e-01 -8.425263e-02 -1.985319e-01 -8.425263e-02
#> [20,] -4.018557e-02 -9.856455e-02 -4.018557e-02 -9.856455e-02
#> [21,] -1.976336e-03 -1.899573e-02 -1.976336e-03 -1.899573e-02
#> [22,] -1.860011e-02 -3.712384e-02 -1.860011e-02 -3.712384e-02
#> [23,] -5.850906e-02 -3.712384e-02 -5.850906e-02 -3.712384e-02
#> [24,] -5.850906e-02 -1.868614e-01 -5.850906e-02 -1.868614e-01
#> [25,] 1.179612e-16 -1.149331e-01 1.179612e-16 -1.149331e-01
#> [26,] 1.179612e-16 -2.651891e-02 1.179612e-16 -2.651891e-02
#> [27,] 1.179612e-16 8.270704e-02 1.179612e-16 8.270704e-02
#> [28,] 1.179612e-16 5.182492e-02 1.179612e-16 5.182492e-02
#> [29,] 1.179612e-16 1.110223e-16 1.179612e-16 1.110223e-16
#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))
#> X.CoastalEstuaries X.OR.NorthCoast X.CoastalEstuaries X.OR.NorthCoast
#> [1,] 9.303347e-02 1.182461e-02 9.303347e-02 1.182461e-02
#> [2,] 6.572282e-01 1.182461e-02 6.572282e-01 1.182461e-02
#> [3,] 8.277781e-01 5.341238e-02 8.277781e-01 5.341238e-02
#> [4,] 4.821278e-01 5.341238e-02 4.821278e-01 5.341238e-02
#> [5,] 4.821278e-01 5.341238e-02 4.821278e-01 5.341238e-02
#> [6,] 1.520998e+00 5.341238e-02 1.520998e+00 5.341238e-02
#> [7,] 6.021581e-01 5.341238e-02 6.021581e-01 5.341238e-02
#> [8,] -6.368764e-01 7.593565e-01 -6.368764e-01 7.593565e-01
#> [9,] -3.247047e-01 1.835007e-01 -3.247047e-01 1.835007e-01
#> [10,] 7.411909e-01 1.966139e-01 7.411909e-01 1.966139e-01
#> [11,] -3.678015e-01 -5.165463e-02 -3.678015e-01 -5.165463e-02
#> [12,] 1.053284e+00 1.268313e+00 1.053284e+00 1.268313e+00
#> [13,] -3.648463e-01 6.801938e-01 -3.648463e-01 6.801938e-01
#> [14,] -4.224274e-01 6.620275e-01 -4.224274e-01 6.620275e-01
#> [15,] 1.376980e-01 -6.941435e-02 1.376980e-01 -6.941435e-02
#> [16,] 1.376980e-01 -2.473822e-01 1.376980e-01 -2.473822e-01
#> [17,] -7.489474e-01 9.061444e-01 -7.489474e-01 9.061444e-01
#> [18,] -7.638818e-01 -3.358483e-01 -7.638818e-01 -3.358483e-01
#> [19,] -1.638406e+00 -7.633879e-01 -1.638406e+00 -7.633879e-01
#> [20,] -3.316358e-01 -8.930639e-01 -3.316358e-01 -8.930639e-01
#> [21,] -1.630993e-02 -1.721146e-01 -1.630993e-02 -1.721146e-01
#> [22,] -1.534994e-01 -3.363680e-01 -1.534994e-01 -3.363680e-01
#> [23,] -4.828523e-01 -3.363680e-01 -4.828523e-01 -3.363680e-01
#> [24,] -4.828523e-01 -1.693095e+00 -4.828523e-01 -1.693095e+00
#> [25,] 9.734875e-16 -1.041374e+00 9.734875e-16 -1.041374e+00
#> [26,] 9.734875e-16 -2.402799e-01 9.734875e-16 -2.402799e-01
#> [27,] 9.734875e-16 7.493838e-01 9.734875e-16 7.493838e-01
#> [28,] 9.734875e-16 4.695701e-01 9.734875e-16 4.695701e-01
#> [29,] 9.734875e-16 1.005940e-15 9.734875e-16 1.005940e-15
#one-step-ahead standardized residuals
MARSSresiduals(fit, type="tt1")$std.residuals
#> [,1] [,2] [,3] [,4] [,5]
#> CoastalEstuaries -0.05418676 -0.206286276 0.566365934 1.005768 NA
#> OR.NorthCoast NA NA -0.001882139 NA NA
#> X.CoastalEstuaries -0.20628628 0.566365934 1.005767860 0.000000 0.2204673
#> X.OR.NorthCoast 0.00000000 -0.001882139 0.000000000 0.000000 0.0000000
#> [,6] [,7] [,8] [,9] [,10]
#> CoastalEstuaries 0.2204673 2.0405379 1.2425972 -0.7982010 -0.8746660
#> OR.NorthCoast NA NA -0.1395806 1.1554718 0.1779561
#> X.CoastalEstuaries 2.0405379 1.2425972 -0.7982010 -0.8746660 1.3133883
#> X.OR.NorthCoast 0.0000000 -0.1395806 1.1554718 0.1779561 0.3438440
#> [,11] [,12] [,13] [,14] [,15]
#> CoastalEstuaries 1.3133883 -1.1010233 1.7849809 -0.3351368 -0.71142906
#> OR.NorthCoast 0.3438440 -0.8332411 1.6166181 0.6909545 1.09517141
#> X.CoastalEstuaries -1.1010233 1.7849809 -0.3351368 -0.7114291 0.00000000
#> X.OR.NorthCoast -0.8332411 1.6166181 0.6909545 1.0951714 0.03597429
#> [,16] [,17] [,18] [,19] [,20]
#> CoastalEstuaries NA 0.5726369 -0.76775282 -0.31659312 -2.3132481
#> OR.NorthCoast 0.03597429 -0.9304752 1.64162505 -0.08266883 -0.6865966
#> X.CoastalEstuaries 0.57263695 -0.7677528 -0.31659312 -2.31324808 -0.4943444
#> X.OR.NorthCoast -0.93047525 1.6416251 -0.08266883 -0.68659660 -1.3199941
#> [,21] [,22] [,23] [,24] [,25]
#> CoastalEstuaries -0.49434441 0.05447392 0.01643296 NA -0.8768965
#> OR.NorthCoast -1.31999410 -0.07484734 NA 0.2166126 -2.1429904
#> X.CoastalEstuaries 0.05447392 0.01643296 0.00000000 -0.8768965 0.0000000
#> X.OR.NorthCoast -0.07484734 0.00000000 0.21661261 -2.1429904 -1.5204468
#> [,26] [,27] [,28] [,29] [,30]
#> CoastalEstuaries NA NA NA NA NA
#> OR.NorthCoast -1.5204468 -0.8263151 0.9151561 0.7476442 NA
#> X.CoastalEstuaries 0.0000000 0.0000000 0.0000000 0.0000000 NA
#> X.OR.NorthCoast -0.8263151 0.9151561 0.7476442 0.0000000 NA