MAR1 State-Space Model

library(MAR1)
#> Loading required package: leaps
#> Loading required package: bestglm
#> Loading required package: tcltk
library(MARSS)

These examples take 1-2 minutes to run

Construct a MAR model

We will use [run.mar()] arguments to set variables and restrictions. See the [run.mar()] page for information. No restrictions on the restriction matrix (all 0.5).

data(L4.mar)
myvar <- c(0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2) # 8 variates, 3 covariates
myres <- matrix(0.5,
  nrow = length(which(myvar == 1)),
  ncol = length(which(myvar != 0))
)

Fit a MAR1 model first.

run1 <- run.mar(L4.mar, variables = myvar, restrictions = myres, search = "exhaustive")
#>  
#>  
#> searching for best-fit model...
#>        ...BEST-FIT MODEL SELECTED
#>             ( search time:  0 minutes 0.6 seconds )
#>  
#>  
#>  
#> identifying 10 lowest AIC models...
#>        ...TOP MODELS RETAINED 
#>  
#>  
#>  
#> bootstrapping best-fit model...
#>        ...BOOTSTRAPPING COMPLETE
#>  
#>  
#> ════════════════════════════════════════════════════════════════════════════════ 
#>  
#> 

Fit state-space model

This allows use to include observation error. control can be passed in to limit the number of iterations run.

ss.fit <- ss.mar1(L4.mar, run1, control = list(maxit = 50))

Compare to best fit model

ss.fit$B
#>               amphipod chaetognath cladoceran calanoid.lg calanoid.sm
#> amphipod     0.2495807  0.00000000  0.0000000   0.2087382   0.0000000
#> chaetognath  0.0000000  0.55938796  0.0000000   0.0000000  -0.1009282
#> cladoceran   0.0000000  0.00000000  0.4927539   0.0000000   0.0000000
#> calanoid.lg -0.4781403 -0.01429697 -0.1894227   0.6181397   0.0000000
#> calanoid.sm -0.3381246  0.00000000  0.0000000   0.1029780   0.4478054
#> cyclopoid    0.0000000  0.00000000  0.0000000   0.0000000   0.0000000
#> diatom       0.0000000  0.00000000  0.0000000  -0.2615612   0.0000000
#> dino         0.3528929  0.00000000 -0.2273365  -0.1961949   0.4100545
#>              cyclopoid       diatom        dino
#> amphipod    -0.1234572  0.148922126  0.00000000
#> chaetognath  0.0000000  0.000000000  0.00000000
#> cladoceran   0.0000000  0.000000000  0.00000000
#> calanoid.lg -0.1680414 -0.007289215  0.00000000
#> calanoid.sm  0.0000000 -0.017631458 -0.07566676
#> cyclopoid    0.6143973  0.000000000  0.00000000
#> diatom       0.0000000  0.276714637  0.00000000
#> dino        -0.4087143  0.000000000  0.21599674
run1$bestfit$B
#>                amphipod chaetognath cladoceran calanoid.lg calanoid.sm
#> amphipod     0.08151863  0.00000000  0.0000000  0.23997271   0.0000000
#> chaetognath  0.00000000  0.48254589  0.0000000  0.00000000  -0.1039079
#> cladoceran   0.00000000  0.00000000  0.4150929  0.00000000   0.0000000
#> calanoid.lg -0.29967953  0.01557868 -0.1365731  0.42833032   0.0000000
#> calanoid.sm -0.26909543  0.00000000  0.0000000  0.10679836   0.3410686
#> cyclopoid    0.00000000  0.00000000  0.0000000  0.00000000   0.0000000
#> diatom       0.00000000  0.00000000  0.0000000 -0.24443109   0.0000000
#> dino         0.10617918  0.00000000 -0.1988748 -0.03007952   0.2226672
#>              cyclopoid      diatom       dino
#> amphipod    -0.1434129  0.04060264  0.0000000
#> chaetognath  0.0000000  0.00000000  0.0000000
#> cladoceran   0.0000000  0.00000000  0.0000000
#> calanoid.lg -0.1935099  0.01865910  0.0000000
#> calanoid.sm  0.0000000 -0.08207196 -0.0731972
#> cyclopoid    0.4558332  0.00000000  0.0000000
#> diatom       0.0000000  0.18078937  0.0000000
#> dino        -0.2830292  0.00000000  0.1804134

Use a known observation error

R <- diag(0.2, 8)
ss.fit <- ss.mar1(L4.mar, run1, model = list(R = R), control = list(maxit = 50))