MARSS Smoothed Residuals
MARSSresiduals_tT.RdCalculates the standardized (or auxiliary) smoothed residuals sensu Harvey, Koopman and Penzer (1998). The expected values and variance for missing (or left-out) data are also returned (Holmes 2014). Not exported. Access this function with MARSSresiduals(object, type="tT"). At time \(t\) (in the returned matrices), the model residuals are for time \(t\), while the state residuals are for the transition from \(t\) to \(t+1\) following the convention in Harvey, Koopman and Penzer (1998).
Usage
MARSSresiduals.tT(object, Harvey = FALSE, normalize = FALSE,
silent = FALSE, fun.kf = c("MARSSkfas", "MARSSkfss"))Arguments
- object
An object of class
marssMLE.- Harvey
TRUE/FALSE. Use the Harvey et al. (1998) algorithm or use the Holmes (2014) algorithm. The values are the same except for missing values.
- normalize
TRUE/FALSE See details.
- silent
If TRUE, don't print inversion warnings.
- fun.kf
Kalman filter function to use. Can be ignored.
Value
A list with the following components
- model.residuals
The the observed smoothed model residuals: data minus the model predictions conditioned on all observed data. This is different than the Kalman filter innovations which use on the data up to time \(t-1\) for the predictions. See details.
- state.residuals
The smoothed state residuals \(\mathbf{x}_{t+1}^T - \mathbf{Z} \mathbf{x}_{t}^T - \mathbf{u}\). The last time step will be NA because the last step would be for T to T+1 (past the end of the data).
- residuals
The residuals conditioned on the observed data. Returned as a (n+m) x T matrix with
model.residualsin rows 1 to n andstate.residualsin rows n+1 to n+m. NAs will appear in rows 1 to n in the places where data are missing.- var.residuals
The joint variance of the model and state residuals conditioned on observed data. Returned as a (n+m) x (n+m) x T matrix. For Harvey=FALSE, this is Holmes (2014) equation 57. For Harvey=TRUE, this is the residual variance in eqn. 24, page 113, in Harvey et al. (1998). They are identical except for missing values, for those Harvey=TRUE returns 0s. For the state residual variance, the last time step will be all NA because the last step would be for T to T+1 (past the end of the data).
- 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. The model standardized residuals associated with the missing data are replaced with NA.- 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. The model marginal residuals associated with the missing data are replaced with NA.- 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 residuals (values in
model.residuals). For unobserved data, this will be 0 if \(\mathbf{R}\) is diagonal but non-zero if \(\mathbf{R}\) is non-diagonal. See details.- 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 details.
- msg
Any warning messages. This will be printed unless Object$control$trace = -1 (suppress all error messages).
Details
This function returns the raw, the Cholesky standardized and the marginal standardized smoothed model and state residuals. 'smoothed' means conditioned on all the observed data and a set of parameters. These are the residuals presented in Harvey, Koopman and Penzer (1998) pages 112-113, with the addition of the values for unobserved data (Holmes 2014). If Harvey=TRUE, the function uses the algorithm on page 112 of Harvey, Koopman and Penzer (1998) to compute the conditional residuals and variance of the residuals. If Harvey=FALSE, the function uses the equations in the technical report (Holmes 2014). Unlike the innovations residuals, the smoothed residuals are autocorrelated (section 4.1 in Harvey and Koopman 1992) and thus an ACF test on these residuals would not reveal model inadequacy.
The residuals matrix has a value for each time step. The residuals in column \(t\) rows 1 to n are the model residuals associated with the data at time \(t\). The residuals in rows n+1 to n+m are the state residuals associated with the transition from \(\mathbf{x}_{t}\) to \(\mathbf{x}_{t+1}\), not the transition from \(\mathbf{x}_{t-1}\) to \(\mathbf{x}_{t}\). Because \(\mathbf{x}_{t+1}\) does not exist at time \(T\), the state residuals and associated variances at time \(T\) are NA.
Below the conditional residuals and their variance are discussed. The random variables are capitalized and the realizations from the random variables are lower case. The random variables are \(\mathbf{X}\), \(\mathbf{Y}\), \(\mathbf{V}\) and \(\mathbf{W}\). There are two types of \(\mathbf{Y}\). The observed \(\mathbf{Y}\) that are used to estimate the states \(\mathbf{x}\). These are termed \(\mathbf{Y}^{(1)}\). The unobserved \(\mathbf{Y}\) are termed \(\mathbf{Y}^{(2)}\). These are not used to estimate the states \(\mathbf{x}\) and we may or may not know the values of \(\mathbf{y}^{(2)}\). Typically we treat \(\mathbf{y}^{(2)}\) as unknown but it may be known but we did not include it in our model fitting. Note that the model parameters \(\Theta\) are treated as fixed or known. The 'fitting' does not involve estimating \(\Theta\); it involves estimating \(\mathbf{x}\). All MARSS parameters can be time varying but the \(t\) subscripts are left off parameters to reduce clutter.
Model residuals
\(\mathbf{v}_{t}\) is the difference between the data and the predicted data at time \(t\) given \(\mathbf{x}_{t}\):
$$ \mathbf{v}_{t} = \mathbf{y}_{t} - \mathbf{Z} \mathbf{x}_{t} - \mathbf{a} - \mathbf{D}\mathbf{d}_t$$
\(\mathbf{x}_{t}\) is unknown (hidden) and our data are one realization of \(\mathbf{y}_{t}\). The observed model residuals \(\hat{\mathbf{v}}_{t}\) are the difference between the observed data and the predicted data at time \(t\) using the fitted model. MARSSresiduals.tT fits the model using all the data, thus
$$ \hat{\mathbf{v}}_{t} = \mathbf{y}_{t} - \mathbf{Z}\mathbf{x}_{t}^T - \mathbf{a} - \mathbf{D}\mathbf{d}_t$$
where \(\mathbf{x}_{t}^T\) is the expected value of \(\mathbf{X}_{t}\) conditioned on the data from 1 to \(T\) (all the data), i.e. the Kalman smoother estimate of the states at time \(t\). \(\mathbf{y}_{t}\) are your data and missing values will appear as NA in the observed model residuals. These are returned as model.residuals and rows 1 to \(n\) of residuals.
res1 and res2 in the code below will be the same.
dat = t(harborSeal)[2:3,]
fit = MARSS(dat)
Z = coef(fit, type="matrix")$Z
A = coef(fit, type="matrix")$A
res1 = dat - Z %*% fit$states - A %*% matrix(1,1,ncol(dat))
res2 = MARSSresiduals(fit, type="tT")$model.residualsState residuals
\(\mathbf{w}_{t+1}\) are the difference between the state at time \(t+1\) and the expected value of the state at time \(t+1\) given the state at time \(t\):
$$ \mathbf{w}_{t+1} = \mathbf{x}_{t+1} - \mathbf{B} \mathbf{x}_{t} - \mathbf{u} - \mathbf{C}\mathbf{c}_{t+1}$$
The estimated state residuals \(\hat{\mathbf{w}}_{t+1}\) are the difference between estimate of \(\mathbf{x}_{t+1}\) minus the estimate using \(\mathbf{x}_{t}\).
$$ \hat{\mathbf{w}}_{t+1} = \mathbf{x}_{t+1}^T - \mathbf{B}\mathbf{x}_{t}^T - \mathbf{u} - \mathbf{C}\mathbf{c}_{t+1}$$
where \(\mathbf{x}_{t+1}^T\) is the Kalman smoother estimate of the states at time \(t+1\) and \(\mathbf{x}_{t}^T\) is the Kalman smoother estimate of the states at time \(t\).
The estimated state residuals \(\mathbf{w}_{t+1}\) are returned in state.residuals and rows \(n+1\) to \(n+m\) of residuals. state.residuals[,t] is \(\mathbf{w}_{t+1}\) (notice time subscript difference). There are no NAs in the estimated state residuals as an estimate of the state exists whether or not there are associated data.
res1 and res2 in the code below will be the same.
dat <- t(harborSeal)[2:3,]
TT <- ncol(dat)
fit <- MARSS(dat)
B <- coef(fit, type="matrix")$B
U <- coef(fit, type="matrix")$U
statestp1 <- MARSSkf(fit)$xtT[,2:TT]
statest <- MARSSkf(fit)$xtT[,1:(TT-1)]
res1 <- statestp1 - B %*% statest - U %*% matrix(1,1,TT-1)
res2 <- MARSSresiduals(fit, type="tT")$state.residuals[,1:(TT-1)]Note that the state residual at the last time step (not shown) will be NA because it is the residual associated with \(\mathbf{x}_T\) to \(\mathbf{x}_{T+1}\) and \(T+1\) is beyond the data. Similarly, the variance matrix at the last time step will have NAs for the same reason.
Variance of the residuals
In a state-space model, \(\mathbf{X}\) and \(\mathbf{Y}\) are stochastic, and the model and state residuals are random variables \(\hat{\mathbf{V}}_{t}\) and \(\hat{\mathbf{W}}_{t+1}\). To evaluate the residuals we observed (with \(\mathbf{y}^{(1)}\)), we use the joint distribution of \(\hat{\mathbf{V}}_{t}, \hat{\mathbf{W}}_{t+1}\) across all the different possible data sets that our MARSS equations with parameters \(\Theta\) might generate. Denote the matrix of \(\hat{\mathbf{V}}_{t}, \hat{\mathbf{W}}_{t+1}\), as \(\widehat{\mathcal{E}}_{t}\). That distribution has an expected value (mean) and variance:
$$ \textrm{E}[\widehat{\mathcal{E}}_{t}] = 0; \textrm{var}[\widehat{\mathcal{E}}_{t}] = \hat{\Sigma}_{t} $$
Our observed residuals (returned in residuals) are one sample from this distribution.
To standardize the observed residuals, we will use \( \hat{\Sigma}_{t} \). \( \hat{\Sigma}_{t} \) is returned in var.residuals. Rows/columns 1 to \(n\) are the conditional variances of the model residuals and rows/columns \(n+1\) to \(n+m\) are the conditional variances of the state residuals. The off-diagonal blocks are the covariances between the two types of residuals.
Standardized residuals
MARSSresiduals will return the Cholesky standardized residuals sensu Harvey et al. (1998) in std.residuals for outlier and shock detection. These are the model and state residuals multiplied by the inverse of the lower triangle of the Cholesky decomposition of var.residuals (note chol() in R returns the upper triangle thus a transpose is needed). The standardized model residuals are set to NA when there are missing data. The standardized state residuals however always exist since the expected value of the states exist without data. The calculation of the standardized residuals for both the observations and states requires the full residuals variance matrix. Since the state residuals variance is NA at the last time step, the standardized residual in the last time step will be all NA (for both model and state residuals).
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 R. The transformed residuals (from this function) for the i-th row of \(\mathbf{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 transformed residuals can look rather non-intuitive.
The marginal standardized residuals are returned in mar.residuals. These are the model and state residuals multiplied by the inverse of the diagonal matrix formed by the square root of the diagonal of var.residuals. These residuals will be correlated (across the residuals at time \(t\)) but are easier to interpret when \(\mathbf{Q}\) and \(\mathbf{R}\) are 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 the state residuals, the Block Cholesky standardization will be different because Block Cholesky standardization treats the model and state residuals as independent (which they are not in the smoothations case).
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).
MARSSresiduals.tT returns the residuals defined as in the first equations. To get the residuals defined as Harvey et al. (1998) define them (second equations), then use normalize=TRUE. In that case the unconditional variance of residuals will be \(\mathbf{I}\) instead of \(\mathbf{Q}\) and \(\mathbf{R}\).
Missing or left-out data
\( \textrm{E}[\widehat{\mathcal{E}}_{t}] \) and \( \textrm{var}[\widehat{\mathcal{E}}_{t}] \) are for the distribution across all possible \(\mathbf{X}\) and \(\mathbf{Y}\). We can also compute the expected value and variance conditioned on a specific value of \(\mathbf{Y}\), the one we observed \(\mathbf{y}^{(1)}\) (Holmes 2014). If there are no missing values, this is not very interesting as \(\textrm{E}[\hat{\mathbf{V}}_{t}|\mathbf{y}^{(1)}]=\hat{\mathbf{v}}_{t}\) and \(\textrm{var}[\hat{\mathbf{V}}_{t}|\mathbf{y}^{(1)}] = 0\). If we have data that are missing because we left them out, however, \(\textrm{E}[\hat{\mathbf{V}}_{t}|\mathbf{y}^{(1)}]\) and \(\textrm{var}[\hat{\mathbf{V}}_{t}|\mathbf{y}^{(1)}]\) are the values we need to evaluate whether the left-out data are unusual relative to what you expect given the data you did collect.
E.obs.residuals is the conditional expected value \(\textrm{E}[\hat{\mathbf{V}}|\mathbf{y}^{(1)}]\) (notice small \(\mathbf{y}\)). It is
$$\textrm{E}[\mathbf{Y}_{t}|\mathbf{y}^{(1)}] - \mathbf{Z}\mathbf{x}_t^T - \mathbf{a} $$
It is similar to \(\hat{\mathbf{v}}_{t}\). The difference is the \(\mathbf{y}\) term. \(\textrm{E}[\mathbf{Y}^{(1)}_{t}|\mathbf{y}^{(1)}] \) is \(\mathbf{y}^{(1)}_{t}\) for the non-missing values. For the missing values, the value depends on \(\mathbf{R}\). If \(\mathbf{R}\) is diagonal, \(\textrm{E}[\mathbf{Y}^{(2)}_{t}|\mathbf{y}^{(1)}] \) is \(\mathbf{Z}\mathbf{x}_t^T + \mathbf{a}\) and the expected residual value is 0. If \(\mathbf{R}\) is non-diagonal however, it will be non-zero.
var.obs.residuals is the conditional variance \(\textrm{var}[\hat{\mathbf{V}}|\mathbf{y}^{(1)}]\) (eqn 24 in Holmes (2014)). For the non-missing values, this variance is 0 since \(\hat{\mathbf{V}}|\mathbf{y}^{(1)}\) is a fixed value. For the missing values, \(\hat{\mathbf{V}}|\mathbf{y}^{(1)}\) is not fixed because \(\mathbf{Y}^{(2)}\) is a random variable. For these values, the variance of \(\hat{\mathbf{V}}|\mathbf{y}^{(1)}\) is determined by the variance of \(\mathbf{Y}^{(2)}\) conditioned on \(\mathbf{Y}^{(1)}=\mathbf{y}^{(1)}\). This variance matrix is returned in var.obs.residuals. The variance of \(\hat{\mathbf{W}}|\mathbf{y}^{(1)}\) is 0 and thus is not included.
The variance \(\textrm{var}[\hat{\mathbf{V}}_{t}|\mathbf{Y}^{(1)}] \) (uppercase \( \mathbf{Y} \)) returned in the 1 to \(n\) rows/columns of var.residuals may also be of interest depending on what you are investigating with regards to missing values. For example, it may be of interest in a simulation study or cases where you have multiple replicated \(\mathbf{Y}\) data sets. var.residuals would allow you to determine if the left-out residuals are unusual with regards to what you would expect for left-out data in that location of the \(\mathbf{Y}\) matrix but not specifically relative to the data you did collect. If \(\mathbf{R}\) is non-diagonal and the \(\mathbf{y}^{(1)}\) and \(\mathbf{y}^{(2)}\) are highly correlated, the variance of \(\textrm{var}[\hat{\mathbf{V}}_{t}|\mathbf{Y}^{(1)}] \) and variance of \(\textrm{var}[\hat{\mathbf{V}}_{t}|\mathbf{y}^{(1)}] \) for the left-out data would be quite different. In the latter, the variance is low because \(\mathbf{y}^{(1)} \) has strong information about \(\mathbf{y}^{(2)} \). In the former, we integrate over \(\mathbf{Y}^{(1)} \) and the variance could be high (depending on the parameters).
Note, if Harvey=TRUE then the rows and columns of var.residuals corresponding to missing values will be NA. This is because the Harvey et al. algorithm does not compute the residual variance for missing values.
References
Harvey, A., S. J. Koopman, and J. Penzer. 1998. Messy time series: a unified approach. Advances in Econometrics 13: 103-144 (see page 112-113). Equation 21 is the Kalman eqns. Eqn 23 and 24 is the backward recursion to compute the smoothations. This function uses the MARSSkf output for eqn 21 and then implements the backwards recursion in equation 23 and equation 24. Pages 120-134 discuss the use of standardized residuals for outlier and structural break detection.
de Jong, P. and J. Penzer. 1998. Diagnosing shocks in time series. Journal of the American Statistical Association 93: 796-806. This one shows the same equations; see eqn 6. This paper mentions the scaling based on the inverse of the sqrt (Cholesky decomposition) of the variance-covariance matrix for the residuals (model and state together). This is in the right column, half-way down on page 800.
Koopman, S. J., N. Shephard, and J. A. Doornik. 1999. Statistical algorithms for models in state space using SsfPack 2.2. Econometrics Journal 2: 113-166. (see pages 147-148).
Harvey, A. and S. J. Koopman. 1992. Diagnostic checking of unobserved-components time series models. Journal of Business & Economic Statistics 4: 377-389.
Holmes, E. E. 2014. Computation of standardized residuals for (MARSS) models. Technical Report. arXiv:1411.0045.
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 residuals
state.resids1 <- MARSSresiduals(fit, type="tT")$state.residuals
#this is the same as hatx_t-(hatx_{t-1}+u)
states <- fit$states
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 the state residuals to a variance of 1
Q <- coef(fit,type="matrix")$Q
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
#Cholesky standardized (by joint variance) model & state residuals
MARSSresiduals(fit, type="tT")$std.residuals
#> [,1] [,2] [,3] [,4] [,5]
#> CoastalEstuaries -0.10919851 -0.69616161 -0.21283712 0.4543333 NA
#> OR.NorthCoast NA NA -0.07076508 NA NA
#> X.CoastalEstuaries 0.08910975 0.65993119 1.31238202 1.2287727 0.8394832
#> X.OR.NorthCoast 0.02307109 0.02307109 0.11400130 0.1341658 0.1341658
#> [,6] [,7] [,8] [,9] [,10]
#> CoastalEstuaries -1.3656819 1.1475600 1.53880072 -0.3874468 -1.32282146
#> OR.NorthCoast NA NA -1.04427738 0.7703940 -0.01732641
#> X.CoastalEstuaries 1.7524119 1.8939213 -0.02706276 -0.8374093 0.36257668
#> X.OR.NorthCoast 0.1341658 0.1341658 0.56014137 0.9033785 0.34052208
#> [,11] [,12] [,13] [,14] [,15]
#> CoastalEstuaries 1.3763064 -1.7637571 1.7612261 0.07191435 -0.7363352
#> OR.NorthCoast 0.3274781 -1.7406921 0.7755500 0.02395583 0.9645447
#> X.CoastalEstuaries 0.3250536 0.5973221 0.5965871 -0.69579724 -0.1463872
#> X.OR.NorthCoast 0.1296784 1.0854444 1.7385635 1.19758302 0.5287523
#> [,16] [,17] [,18] [,19] [,20]
#> CoastalEstuaries NA 1.1655754 0.01865192 1.0861153 -1.6220364
#> OR.NorthCoast 0.2346844 -1.5211471 1.63783211 0.5638654 0.1711639
#> X.CoastalEstuaries 0.2397652 -0.3929017 -1.33417849 -2.0924123 -1.7157320
#> X.OR.NorthCoast -0.2826585 0.5877872 0.50877259 -0.9812779 -1.4814897
#> [,21] [,22] [,23] [,24] [,25]
#> CoastalEstuaries -0.3916652 0.1715189 0.4373641 NA -0.8768965
#> OR.NorthCoast -0.9572601 0.2281407 NA 1.8844534 -0.8653923
#> X.CoastalEstuaries -0.3025196 -0.1537960 -0.7602172 -0.8768965 0.0000000
#> X.OR.NorthCoast -0.9453488 -0.5602184 -0.6002422 -1.7441471 -2.4955040
#> [,26] [,27] [,28] [,29] [,30]
#> CoastalEstuaries NA NA NA NA NA
#> OR.NorthCoast -1.057835 -1.3091723 0.3772807 0.7476442 NA
#> X.CoastalEstuaries 0.000000 0.0000000 0.0000000 0.0000000 NA
#> X.OR.NorthCoast -1.148606 0.4686348 1.1198848 0.0000000 NA
# Returns residuals in a data frame in long form
residuals(fit, type="tT")
#> type .rownames name t value .fitted .resids .sigma
#> 1 ytT CoastalEstuaries model 1 7.434848 7.443643 -8.794738e-03 0.08053900
#> 2 ytT CoastalEstuaries model 2 7.462789 7.516263 -5.347372e-02 0.07681222
#> 3 ytT CoastalEstuaries model 3 7.641084 7.657249 -1.616452e-02 0.07594786
#> 4 ytT CoastalEstuaries model 4 7.851661 7.818901 3.276033e-02 0.07210638
#> 5 ytT CoastalEstuaries model 5 NA 7.938669 NA NA
#> 6 ytT CoastalEstuaries model 6 7.959975 8.058437 -9.846296e-02 0.07209801
#> 7 ytT CoastalEstuaries model 7 8.391176 8.304090 8.708663e-02 0.07588852
#> 8 ytT CoastalEstuaries model 8 8.555837 8.438403 1.174343e-01 0.07631544
#> 9 ytT CoastalEstuaries model 9 8.392990 8.422577 -2.958727e-02 0.07636473
#> 10 ytT CoastalEstuaries model 10 8.343554 8.444578 -1.010244e-01 0.07637036
#> 11 ytT CoastalEstuaries model 11 8.700847 8.595738 1.051090e-01 0.07637036
#> 12 ytT CoastalEstuaries model 12 8.477828 8.612517 -1.346888e-01 0.07636473
#> 13 ytT CoastalEstuaries model 13 8.935904 8.801495 1.344087e-01 0.07631544
#> 14 ytT CoastalEstuaries model 14 8.824089 8.818632 5.457470e-03 0.07588848
#> 15 ytT CoastalEstuaries model 15 8.775704 8.828792 -5.308803e-02 0.07209765
#> 16 ytT CoastalEstuaries model 16 NA 8.906824 NA NA
#> 17 ytT CoastalEstuaries model 17 9.068892 8.984857 8.403523e-02 0.07209763
#> 18 ytT CoastalEstuaries model 18 8.956866 8.955451 1.415464e-03 0.07588838
#> 19 ytT CoastalEstuaries model 19 9.007122 8.924236 8.288641e-02 0.07631456
#> 20 ytT CoastalEstuaries model 20 8.663196 8.787051 -1.238542e-01 0.07635722
#> 21 ytT CoastalEstuaries model 21 8.778326 8.808212 -2.988622e-02 0.07630554
#> 22 ytT CoastalEstuaries model 22 8.880586 8.867583 1.300266e-02 0.07580893
#> 23 ytT CoastalEstuaries model 23 8.941545 8.910330 3.121569e-02 0.07137231
#> 24 ytT CoastalEstuaries model 24 NA 8.913168 NA NA
#> 25 ytT CoastalEstuaries model 25 8.870242 8.916006 -4.576419e-02 0.05218881
#> 26 ytT CoastalEstuaries model 26 NA 8.977353 NA NA
#> 27 ytT CoastalEstuaries model 27 NA 9.038700 NA NA
#> 28 ytT CoastalEstuaries model 28 NA 9.100047 NA NA
#> 29 ytT CoastalEstuaries model 29 NA 9.161394 NA NA
#> 30 ytT CoastalEstuaries model 30 NA 9.222741 NA NA
#> 31 ytT OR.NorthCoast model 1 NA 6.322973 NA NA
#> 32 ytT OR.NorthCoast model 2 NA 6.375274 NA NA
#> 33 ytT OR.NorthCoast model 3 6.423247 6.427575 -4.327605e-03 0.06115453
#> 34 ytT OR.NorthCoast model 4 NA 6.484465 NA NA
#> 35 ytT OR.NorthCoast model 5 NA 6.541356 NA NA
#> 36 ytT OR.NorthCoast model 6 NA 6.598247 NA NA
#> 37 ytT OR.NorthCoast model 7 NA 6.655137 NA NA
#> 38 ytT OR.NorthCoast model 8 6.638568 6.712028 -7.346024e-02 0.07034552
#> 39 ytT OR.NorthCoast model 9 6.906755 6.846831 5.992331e-02 0.07778267
#> 40 ytT OR.NorthCoast model 10 6.916715 6.918080 -1.364562e-03 0.07875619
#> 41 ytT OR.NorthCoast model 11 7.016610 6.990775 2.583472e-02 0.07888991
#> 42 ytT OR.NorthCoast model 12 6.898715 7.036070 -1.373552e-01 0.07890839
#> 43 ytT OR.NorthCoast model 13 7.288244 7.227045 6.119938e-02 0.07891095
#> 44 ytT OR.NorthCoast model 14 7.355002 7.353112 1.890386e-03 0.07891130
#> 45 ytT OR.NorthCoast model 15 7.553287 7.477173 7.611352e-02 0.07891135
#> 46 ytT OR.NorthCoast model 16 7.539027 7.520508 1.851926e-02 0.07891133
#> 47 ytT OR.NorthCoast model 17 7.424165 7.544201 -1.200355e-01 0.07891116
#> 48 ytT OR.NorthCoast model 18 7.824446 7.695205 1.292412e-01 0.07890993
#> 49 ytT OR.NorthCoast model 19 7.753624 7.709134 4.448958e-02 0.07890106
#> 50 ytT OR.NorthCoast model 20 7.689371 7.675877 1.349403e-02 0.07883686
#> 51 ytT OR.NorthCoast model 21 7.553287 7.628308 -7.502167e-02 0.07837125
#> 52 ytT OR.NorthCoast model 22 7.677400 7.660308 1.709214e-02 0.07491929
#> 53 ytT OR.NorthCoast model 23 NA 7.674180 NA NA
#> 54 ytT OR.NorthCoast model 24 7.829233 7.688052 1.411804e-01 0.07491851
#> 55 ytT OR.NorthCoast model 25 7.484369 7.552186 -6.781780e-02 0.07836654
#> 56 ytT OR.NorthCoast model 26 7.404888 7.488249 -8.336153e-02 0.07880388
#> 57 ytT OR.NorthCoast model 27 7.409742 7.512726 -1.029840e-01 0.07866342
#> 58 ytT OR.NorthCoast model 28 7.675546 7.646429 2.911729e-02 0.07717672
#> 59 ytT OR.NorthCoast model 29 7.798113 7.749249 4.886326e-02 0.06535630
#> 60 ytT OR.NorthCoast model 30 NA 7.800245 NA NA
#> 61 xtT X.CoastalEstuaries state 1 7.443643 7.443613 1.127322e-02 0.08664366
#> 62 xtT X.CoastalEstuaries state 2 7.516263 7.504990 7.963885e-02 0.08501149
#> 63 xtT X.CoastalEstuaries state 3 7.657249 7.577610 1.003050e-01 0.08338868
#> 64 xtT X.CoastalEstuaries state 4 7.818901 7.718596 5.842127e-02 0.06959195
#> 65 xtT X.CoastalEstuaries state 5 7.938669 7.880248 5.842127e-02 0.06959195
#> 66 xtT X.CoastalEstuaries state 6 8.058437 8.000016 1.843052e-01 0.08336761
#> 67 xtT X.CoastalEstuaries state 7 8.304090 8.119785 7.296579e-02 0.08481890
#> 68 xtT X.CoastalEstuaries state 8 8.438403 8.365437 -7.717274e-02 0.08498541
#> 69 xtT X.CoastalEstuaries state 9 8.422577 8.499750 -3.934571e-02 0.08500464
#> 70 xtT X.CoastalEstuaries state 10 8.444578 8.483924 8.981293e-02 0.08500662
#> 71 xtT X.CoastalEstuaries state 11 8.595738 8.505925 -4.456791e-02 0.08500464
#> 72 xtT X.CoastalEstuaries state 12 8.612517 8.657085 1.276304e-01 0.08498540
#> 73 xtT X.CoastalEstuaries state 13 8.801495 8.673864 -4.420982e-02 0.08481888
#> 74 xtT X.CoastalEstuaries state 14 8.818632 8.862842 -5.118714e-02 0.08336747
#> 75 xtT X.CoastalEstuaries state 15 8.828792 8.879979 1.668539e-02 0.06959052
#> 76 xtT X.CoastalEstuaries state 16 8.906824 8.890139 1.668539e-02 0.06959052
#> 77 xtT X.CoastalEstuaries state 17 8.984857 8.968171 -9.075282e-02 0.08336743
#> 78 xtT X.CoastalEstuaries state 18 8.955451 9.046204 -9.256248e-02 0.08481854
#> 79 xtT X.CoastalEstuaries state 19 8.924236 9.016798 -1.985319e-01 0.08498247
#> 80 xtT X.CoastalEstuaries state 20 8.787051 8.985583 -4.018557e-02 0.08497931
#> 81 xtT X.CoastalEstuaries state 21 8.808212 8.848398 -1.976336e-03 0.08478766
#> 82 xtT X.CoastalEstuaries state 22 8.867583 8.869559 -1.860011e-02 0.08309542
#> 83 xtT X.CoastalEstuaries state 23 8.910330 8.928930 -5.850906e-02 0.06672288
#> 84 xtT X.CoastalEstuaries state 24 8.913168 8.971677 -5.850906e-02 0.06672288
#> 85 xtT X.CoastalEstuaries state 25 8.916006 8.974515 1.179612e-16 0.00000000
#> 86 xtT X.CoastalEstuaries state 26 8.977353 8.977353 1.179612e-16 0.00000000
#> 87 xtT X.CoastalEstuaries state 27 9.038700 9.038700 1.179612e-16 0.00000000
#> 88 xtT X.CoastalEstuaries state 28 9.100047 9.100047 1.179612e-16 0.00000000
#> 89 xtT X.CoastalEstuaries state 29 9.161394 9.161394 1.179612e-16 0.00000000
#> 90 xtT X.CoastalEstuaries state 30 9.222741 9.222741 NA NA
#> 91 xtT X.OR.NorthCoast state 1 6.322973 6.321668 1.305044e-03 0.05656620
#> 92 xtT X.OR.NorthCoast state 2 6.375274 6.373969 1.305044e-03 0.05656620
#> 93 xtT X.OR.NorthCoast state 3 6.427575 6.426270 5.894950e-03 0.04393781
#> 94 xtT X.OR.NorthCoast state 4 6.484465 6.478570 5.894950e-03 0.04393781
#> 95 xtT X.OR.NorthCoast state 5 6.541356 6.535461 5.894950e-03 0.04393781
#> 96 xtT X.OR.NorthCoast state 6 6.598247 6.592352 5.894950e-03 0.04393781
#> 97 xtT X.OR.NorthCoast state 7 6.655137 6.649242 5.894950e-03 0.04393781
#> 98 xtT X.OR.NorthCoast state 8 6.712028 6.706133 8.380770e-02 0.07121834
#> 99 xtT X.OR.NorthCoast state 9 6.846831 6.763024 2.025237e-02 0.07420662
#> 100 xtT X.OR.NorthCoast state 10 6.918080 6.897827 2.169964e-02 0.07461057
#> 101 xtT X.OR.NorthCoast state 11 6.990775 6.969075 -5.700953e-03 0.07466628
#> 102 xtT X.OR.NorthCoast state 12 7.036070 7.041771 1.399795e-01 0.07467398
#> 103 xtT X.OR.NorthCoast state 13 7.227045 7.087065 7.507078e-02 0.07467505
#> 104 xtT X.OR.NorthCoast state 14 7.353112 7.278041 7.306581e-02 0.07467520
#> 105 xtT X.OR.NorthCoast state 15 7.477173 7.404107 -7.661036e-03 0.07467521
#> 106 xtT X.OR.NorthCoast state 16 7.520508 7.528169 -2.730277e-02 0.07467514
#> 107 xtT X.OR.NorthCoast state 17 7.544201 7.571504 1.000082e-01 0.07467463
#> 108 xtT X.OR.NorthCoast state 18 7.695205 7.595197 -3.706648e-02 0.07467093
#> 109 xtT X.OR.NorthCoast state 19 7.709134 7.746200 -8.425263e-02 0.07464418
#> 110 xtT X.OR.NorthCoast state 20 7.675877 7.760130 -9.856455e-02 0.07445051
#> 111 xtT X.OR.NorthCoast state 21 7.628308 7.726873 -1.899573e-02 0.07303521
#> 112 xtT X.OR.NorthCoast state 22 7.660308 7.679304 -3.712384e-02 0.06184810
#> 113 xtT X.OR.NorthCoast state 23 7.674180 7.711304 -3.712384e-02 0.06184810
#> 114 xtT X.OR.NorthCoast state 24 7.688052 7.725176 -1.868614e-01 0.07303346
#> 115 xtT X.OR.NorthCoast state 25 7.552186 7.739048 -1.149331e-01 0.07443696
#> 116 xtT X.OR.NorthCoast state 26 7.488249 7.603182 -2.651891e-02 0.07454540
#> 117 xtT X.OR.NorthCoast state 27 7.512726 7.539245 8.270704e-02 0.07395312
#> 118 xtT X.OR.NorthCoast state 28 7.646429 7.563722 5.182492e-02 0.06931762
#> 119 xtT X.OR.NorthCoast state 29 7.749249 7.697424 1.110223e-16 0.00000000
#> 120 xtT X.OR.NorthCoast state 30 7.800245 7.800245 NA NA
#> .std.resids
#> 1 -0.10919851
#> 2 -0.69616161
#> 3 -0.21283712
#> 4 0.45433331
#> 5 NA
#> 6 -1.36568194
#> 7 1.14756000
#> 8 1.53880072
#> 9 -0.38744679
#> 10 -1.32282146
#> 11 1.37630639
#> 12 -1.76375707
#> 13 1.76122614
#> 14 0.07191435
#> 15 -0.73633516
#> 16 NA
#> 17 1.16557539
#> 18 0.01865192
#> 19 1.08611528
#> 20 -1.62203640
#> 21 -0.39166517
#> 22 0.17151891
#> 23 0.43736411
#> 24 NA
#> 25 -0.87689650
#> 26 NA
#> 27 NA
#> 28 NA
#> 29 NA
#> 30 NA
#> 31 NA
#> 32 NA
#> 33 -0.07076508
#> 34 NA
#> 35 NA
#> 36 NA
#> 37 NA
#> 38 -1.04427738
#> 39 0.77039405
#> 40 -0.01732641
#> 41 0.32747812
#> 42 -1.74069213
#> 43 0.77554996
#> 44 0.02395583
#> 45 0.96454465
#> 46 0.23468445
#> 47 -1.52114708
#> 48 1.63783211
#> 49 0.56386545
#> 50 0.17116394
#> 51 -0.95726008
#> 52 0.22814069
#> 53 NA
#> 54 1.88445336
#> 55 -0.86539229
#> 56 -1.05783527
#> 57 -1.30917232
#> 58 0.37728072
#> 59 0.74764422
#> 60 NA
#> 61 0.08910975
#> 62 0.65993119
#> 63 1.31238202
#> 64 1.22877268
#> 65 0.83948322
#> 66 1.75241192
#> 67 1.89392130
#> 68 -0.02706276
#> 69 -0.83740925
#> 70 0.36257668
#> 71 0.32505359
#> 72 0.59732212
#> 73 0.59658708
#> 74 -0.69579724
#> 75 -0.14638718
#> 76 0.23976524
#> 77 -0.39290173
#> 78 -1.33417849
#> 79 -2.09241230
#> 80 -1.71573201
#> 81 -0.30251961
#> 82 -0.15379600
#> 83 -0.76021717
#> 84 -0.87689650
#> 85 0.00000000
#> 86 0.00000000
#> 87 0.00000000
#> 88 0.00000000
#> 89 0.00000000
#> 90 NA
#> 91 0.02307109
#> 92 0.02307109
#> 93 0.11400130
#> 94 0.13416576
#> 95 0.13416576
#> 96 0.13416576
#> 97 0.13416576
#> 98 0.56014137
#> 99 0.90337853
#> 100 0.34052208
#> 101 0.12967839
#> 102 1.08544439
#> 103 1.73856348
#> 104 1.19758302
#> 105 0.52875235
#> 106 -0.28265847
#> 107 0.58778717
#> 108 0.50877259
#> 109 -0.98127794
#> 110 -1.48148973
#> 111 -0.94534880
#> 112 -0.56021838
#> 113 -0.60024219
#> 114 -1.74414711
#> 115 -2.49550401
#> 116 -1.14860580
#> 117 0.46863478
#> 118 1.11988476
#> 119 0.00000000
#> 120 NA