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
<- list("z11", 0, 0, "z21", "z22", 0, "z31", "z32", "z33",
Z_vals "z41", "z42", "z43", "z51", "z52", "z53")
<- matrix(Z_vals, nrow = N_ts, ncol = 3, byrow = TRUE)
ZZ 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
<- "zero"
aa ## 'DD' and 'd' are for covariates
<- "zero" # matrix(0,mm,1)
DD <- "zero" # matrix(0,1,wk_last)
dd ## 'RR' is var-cov matrix for obs errors
<- "diagonal and unequal" RR
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
<- 3
mm ## 'BB' is identity: 1's along the diagonal & 0's elsewhere
<- "identity" # diag(mm)
BB ## 'uu' is a column vector of 0's
<- "zero" # matrix(0, mm, 1)
uu ## 'CC' and 'cc' are for covariates
<- "zero" # matrix(0, mm, 1)
CC <- "zero" # matrix(0, 1, wk_last)
cc ## 'QQ' is identity
<- "identity" # diag(mm) QQ
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:
- A list of specifications for the model’s vectors and matrices;
- A list of any initial values –
MARSS
will pick its own otherwise; - A list of control parameters for the
MARSS()
function.
## list with specifications for model vectors/matrices
<- list(Z = ZZ, A = aa, D = DD, d = dd, R = RR, B = BB,
mod_list U = uu, C = CC, c = cc, Q = QQ)
## list with model inits
<- list(x0 = matrix(rep(0, mm), mm, 1))
init_list ## list with model control parameters
<- list(maxit = 3000, allow.degen = TRUE) con_list
Now we can fit the model.
## fit MARSS
<- MARSS(y = dat, model = mod_list, inits = init_list,
dfa_1 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.