In nwfsc-timeseries/MARSS: Multivariate Autoregressive State-Space Modeling

knitr::opts_chunk$set(echo = FALSE)

This is a vignette to run the plotting functions and check that outputs look ok.

Load MARSS and check that I am checking the right version.

library(MARSS)
library(ggplot2)
packageVersion("MARSS")
plottypes <- eval(formals(MARSS:::autoplot.marssMLE)$plot.type)
plottypes <- plottypes[!(plottypes %in% c("all", "residuals"))]

marssMLE objects

# Little harder model 

dat <- t(harborSealWA)
dat <- dat[2:4, ] # remove the year row
Qvals <- list("unconstrained", "diagonal and equal", "equalvarcov", "zero")[[1]]
Bvals <- c("identity", "diagonal and unequal")[1]
Rvals <- list("unconstrained", "diagonal and equal", "equalvarcov", "zero")[[2]]

for (Q in Qvals) {
  for (B in Bvals) {
    for (R in Rvals) {
      if (B == "diagonal and unequal" && (is.list(Q) || is.list(R))) next
      mod <- list(Q = Q, Z = "identity", R = R, B = B, U = "zero", x0 = dat[, 1, drop = FALSE] * 1.1)
      if (B != "identity" && R != "zero") mod$tinitx <- 1
      if (Q == "zero" && R == "zero") next
      if (Q == "zero" && B != "identity") next
      kemfit1 <- MARSS(dat, model = mod, fun.kf = "MARSSkfss", silent = TRUE)
      kemfit2 <- MARSS(dat, model = mod, fun.kf = "MARSSkfas", silent = TRUE)
      for(plt in plottypes){
      cat(paste0("\n", plt, " Q:", Q, " R:", R, " B:", B))
      conf.int <- TRUE
      if(plt %in% c("ytt", "ytt1")) conf.int <- FALSE
      a1 <- autoplot(kemfit1, plot.type=plt, silent=TRUE, conf.int = conf.int)[[1]]
      a1 <- autoplot(kemfit2, plot.type=plt, silent=TRUE, conf.int = conf.int)[[1]]
      }
    }
  }
}

Easy model

set.seed(123)
test.num <- "easy1"
dat <- cumsum(rnorm(20))
mod.list <- list(tinitx = 1, U = "zero", R = "unconstrained", Q = "unconstrained", B = "unconstrained", x0 = matrix(dat[1] * 1.0001))
kemfit1 <- MARSS(dat, model = mod.list, silent = TRUE)
kemfit2 <- MARSS(dat, model = list(tinitx = 1, U = "zero", R = "unconstrained", Q = "unconstrained", B = "unconstrained", x0 = matrix(dat[1] * 1.0001)), fun.kf = "MARSSkfss", silent = TRUE)

for(plt in plottypes){
  conf.int <- TRUE
  if(plt %in% c("ytt", "ytt1")) conf.int <- FALSE
    a1 <- autoplot(kemfit1, plot.type=plt, silent=TRUE, conf.int = conf.int)[[1]]
    a1 <- autoplot(kemfit2, plot.type=plt, silent=TRUE, conf.int = conf.int)[[1]]
}

test B unconstrained

Q <- "diagonal and unequal"
B <- "unconstrained"
R <- "diagonal and unequal"
mod <- list(Q = Q, Z = "identity", R = R, B = B, U = "zero", x0 = "unequal", tinitx = 1)
kemfit1 <- MARSS(dat, model = mod, fun.kf = "MARSSkfss", silent = TRUE)
kemfit2 <- MARSS(dat, model = mod, fun.kf = "MARSSkfas", silent = TRUE)

for(plt in plottypes){
cat(paste0("\n", plt, " Q:", Q, " R:", R, " B:", B))
if (plt=="residuals") cat("MARSSkfss")
plot(kemfit1, plot.type=plt, silent=TRUE)
if (plt=="residuals") cat("MARSSkfas")
plot(kemfit2, plot.type=plt, silent=TRUE)
}

test Q with some zeros

Q <- ldiag(list("q1", 0, "q2"))
B <- "identity"
R <- "diagonal and equal"
mod <- list(Q = Q, Z = "identity", R = R, B = B, U = "zero", x0 = "unequal", tinitx = 1)
kemfit1 <- MARSS(dat, model = mod, fun.kf = "MARSSkfss", silent = TRUE)
kemfit2 <- MARSS(dat, model = mod, fun.kf = "MARSSkfas", silent = TRUE)

for(plt in plottypes){
cat(paste0("\n", plt, " Q:", Q, " R:", R, " B:", B))
if (plt=="residuals") cat("MARSSkfss")
plot(kemfit1, plot.type=plt, silent=TRUE)
if (plt=="residuals") cat("MARSSkfas")
plot(kemfit2, plot.type=plt, silent=TRUE)
}

test R with some zeros

R <- ldiag(list("r1", 0, "r2"))
B <- "identity"
Q <- "diagonal and equal"
mod <- list(Q = Q, Z = "identity", R = R, B = B, U = "zero", x0 = "unequal", tinitx = 1)
kemfit1 <- MARSS(dat, model = mod, fun.kf = "MARSSkfss", silent = TRUE)
kemfit2 <- MARSS(dat, model = mod, fun.kf = "MARSSkfas", silent = TRUE)

for(plt in plottypes){
cat(paste0("\n", plt, " Q:", Q, " R:", R, " B:", B))
if (plt=="residuals") cat("MARSSkfss")
plot(kemfit1, plot.type=plt, silent=TRUE)
if (plt=="residuals") cat("MARSSkfas")
plot(kemfit2, plot.type=plt, silent=TRUE)
}