10.6 Fitting DFA models with the MARSS package

The MARSS package is designed to work with the fully specified matrix form of the multivariate state-space model we wrote out in Sec 3. Thus, we will need to create a model list with forms for each of the vectors and matrices. Note that even though some of the model elements are scalars and vectors, we will need to specify everything as a matrix (or array for time series of matrices).

Notice that the code below uses some of the MARSS shortcuts for specifying forms of vectors and matrices. We will also use the matrix(list(),nrow,ncol) trick we learned previously.

10.6.1 The observation model

Here we will fit the DFA model above where we have R N_ts observed time series and we want 3 hidden states. Now we need to set up the observation model for MARSS. Here are the vectors and matrices for our first model where each nutrient follows its own process. Recall that we will need to set the elements in the upper R corner of \(\mathbf{Z}\) to 0. We will assume that the observation errors have different variances and they are independent of one another.

## 'ZZ' is loadings matrix
Z_vals <- list("z11", 0, 0, "z21", "z22", 0, "z31", "z32", "z33", 
    "z41", "z42", "z43", "z51", "z52", "z53")
ZZ <- matrix(Z_vals, nrow = N_ts, ncol = 3, byrow = TRUE)
ZZ
     [,1]  [,2]  [,3] 
[1,] "z11" 0     0    
[2,] "z21" "z22" 0    
[3,] "z31" "z32" "z33"
[4,] "z41" "z42" "z43"
[5,] "z51" "z52" "z53"
## 'aa' is the offset/scaling
aa <- "zero"
## 'DD' and 'd' are for covariates
DD <- "zero"  # matrix(0,mm,1)
dd <- "zero"  # matrix(0,1,wk_last)
## 'RR' is var-cov matrix for obs errors
RR <- "diagonal and unequal"

10.6.2 The process model

We need to specify the explicit form for all of the vectors and matrices in the full form of the MARSS model we defined in Sec 3.1. Note that we do not have to specify anything for the states \((\mathbf{x})\) – those are elements that MARSS will identify and estimate itself based on our definitions of the other vectors and matrices.

## number of processes
mm <- 3
## 'BB' is identity: 1's along the diagonal & 0's elsewhere
BB <- "identity"  # diag(mm)
## 'uu' is a column vector of 0's
uu <- "zero"  # matrix(0, mm, 1)
## 'CC' and 'cc' are for covariates
CC <- "zero"  # matrix(0, mm, 1)
cc <- "zero"  # matrix(0, 1, wk_last)
## 'QQ' is identity
QQ <- "identity"  # diag(mm)

10.6.3 Fit the model in MARSS

Now it’s time to fit our first DFA model To do so, we need to create three lists that we will need to pass to the MARSS() function:

  1. A list of specifications for the model’s vectors and matrices;
  2. A list of any initial values – MARSS will pick its own otherwise;
  3. A list of control parameters for the MARSS() function.
## list with specifications for model vectors/matrices
mod_list <- list(Z = ZZ, A = aa, D = DD, d = dd, R = RR, B = BB, 
    U = uu, C = CC, c = cc, Q = QQ)
## list with model inits
init_list <- list(x0 = matrix(rep(0, mm), mm, 1))
## list with model control parameters
con_list <- list(maxit = 3000, allow.degen = TRUE)

Now we can fit the model.

## fit MARSS
dfa_1 <- MARSS(y = dat, model = mod_list, inits = init_list, 
    control = con_list)
Success! abstol and log-log tests passed at 246 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 246 iterations. 
Log-likelihood: -692.9795 
AIC: 1425.959   AICc: 1427.42   
 
                            Estimate
Z.z11                         0.2738
Z.z21                         0.4487
Z.z31                         0.3170
Z.z41                         0.4107
Z.z51                         0.2553
Z.z22                         0.3608
Z.z32                        -0.3690
Z.z42                        -0.0990
Z.z52                        -0.3793
Z.z33                         0.0185
Z.z43                        -0.1404
Z.z53                         0.1317
R.(Cryptomonas,Cryptomonas)   0.1638
R.(Diatoms,Diatoms)           0.2913
R.(Greens,Greens)             0.8621
R.(Unicells,Unicells)         0.3080
R.(Other.algae,Other.algae)   0.5000
x0.X1                         0.2218
x0.X2                         1.8155
x0.X3                        -4.8097
Initial states (x0) defined at t=0

Standard errors have not been calculated. 
Use MARSSparamCIs to compute CIs and bias estimates.