8.3 Observation-error only model
We can estimate the effect of the covariates using a process-error only model, an observation-error only model, or a model with both types of error. An observation-error only model is a multivariate regression, and we will start here so you see the relationship of MARSS model to more familiar linear regression models.
In a standard multivariate linear regression, we only have an observation model with independent errors (the state process does not appear in the model): \[\begin{equation} \mathbf{y}_t = \mathbf{a} + \mathbf{D}\mathbf{d}_t + \mathbf{v}_t, \text{ where } \mathbf{v}_t \sim \text{MVN}(0,\mathbf{R}) \tag{8.2} \end{equation}\] The elements in \(\mathbf{a}\) are the intercepts and those in \(\mathbf{D}\) are the slopes (effects). We have dropped the \(t\) subscript on \(\mathbf{a}\) and \(\mathbf{D}\) because these will be modeled as time-constant. Writing this out for the two plankton and the two covariates we get: \[\begin{equation} \begin{split} \begin{bmatrix} y_{g} \\ y_{bg} \end{bmatrix}_t &= \begin{bmatrix} a_1 \\ a_2 \end{bmatrix} + \begin{bmatrix} \beta_{\mathrm{g,temp}}&\beta_{\mathrm{g,tp}} \\ \beta_{\mathrm{bg,temp}}&\beta_{\mathrm{bg,tp}} \end{bmatrix} \begin{bmatrix} \mathrm{temp} \\ \mathrm{tp} \end{bmatrix}_{t} + \begin{bmatrix} v_{1} \\ v_{2} \end{bmatrix}_t \end{split} \tag{8.3} \end{equation}\]
Let’s fit this model with MARSS. The \(\mathbf{x}\) part of the model is irrelevant so we want to fix the parameters in that part of the model. We won’t set \(\mathbf{B}=0\) or \(\mathbf{Z}=0\) since that might cause numerical issues for the Kalman filter. Instead we fix them as identity matrices and fix \(\mathbf{x}_0=0\) so that \(\mathbf{x}_t=0\) for all \(t\).
<- U <- x0 <- "zero"
Q <- Z <- "identity"
B <- covariates
d <- "zero"
A <- "unconstrained"
D <- dat # to show relationship between dat & the equation
y <- list(B = B, U = U, Q = Q, Z = Z, A = A, D = D,
model.list d = d, x0 = x0)
<- MARSS(y, model = model.list) kem
Success! algorithm run for 15 iterations. abstol and log-log tests passed.
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
Algorithm ran 15 (=minit) iterations and convergence was reached.
Log-likelihood: -276.4287
AIC: 562.8573 AICc: 563.1351
Estimate
R.diag 0.706
D.(Greens,Temp) 0.367
D.(Bluegreens,Temp) 0.392
D.(Greens,TP) 0.058
D.(Bluegreens,TP) 0.535
Initial states (x0) defined at t=0
Standard errors have not been calculated.
Use MARSSparamCIs to compute CIs and bias estimates.
We set A="zero"
because the data and covariates have been demeaned. Of course, one can do multiple regression in R using, say, lm()
, and that would be much, much faster. The EM algorithm is over-kill here, but it is shown so that you see how a standard multivariate linear regression model is written as a MARSS model in matrix form.