7.8 Set up the hypotheses as different models

Only the \(\mathbf{Z}\) matrices change for our model. We will set up a base model list used for all models.

mod.list <- list(
  B = "identity",
  U = "unequal",
  Q = "equalvarcov",
  Z = "placeholder",
  A = "scaling",
  R = "diagonal and equal",
  x0 = "unequal",
  tinitx = 0
)

Then we set up the \(\mathbf{Z}\) matrices using the factor short-cut.

Z.models <- list(
  H1 = factor(c("pnw", "pnw", rep("ps", 4), "ca", "ca", "pnw", "pnw", "ps")),
  H2 = factor(c(rep("coast", 2), rep("ps", 4), rep("coast", 4), "ps")),
  H3 = factor(c(rep("N", 6), "S", "S", "N", "S", "N")),
  H4 = factor(c("nc", "nc", "is", "is", "ps", "ps", "sc", "sc", "nc", "sc", "is")),
  H5 = factor(rep("pan", 11)),
  H6 = factor(1:11) # site
)
names(Z.models) <-
  c("stock", "coast+PS", "N+S", "NC+strait+PS+SC", "panmictic", "site")

7.8.1 Fit the models

We loop through the models, fit and store the results:

out.tab <- NULL
fits <- list()
for (i in 1:length(Z.models)) {
    mod.list$Z <- Z.models[[i]]
    fit <- MARSS::MARSS(sealData, model = mod.list, silent = TRUE, 
        control = list(maxit = 1000))
    out <- data.frame(H = names(Z.models)[i], logLik = fit$logLik, 
        AICc = fit$AICc, num.param = fit$num.params, m = length(unique(Z.models[[i]])), 
        num.iter = fit$numIter, converged = !fit$convergence)
    out.tab <- rbind(out.tab, out)
    fits <- c(fits, list(fit))
}

We will use AICc and AIC weights to summarize the data support for the different hypotheses. First we will sort the fits based on AICc:

min.AICc <- order(out.tab$AICc)
out.tab.1 <- out.tab[min.AICc, ]

Next we add the \(\Delta\)AICc values by subtracting the lowest AICc:

out.tab.1 <- cbind(out.tab.1, delta.AICc = out.tab.1$AICc - out.tab.1$AICc[1])

Relative likelihood is defined as \(\,\text{exp}(-\Delta \mathrm{AICc}/2)\).

out.tab.1 <- cbind(out.tab.1, rel.like = exp(-1 * out.tab.1$delta.AICc/2))

The AIC weight for a model is its relative likelihood divided by the sum of all the relative likelihoods.

out.tab.1 <- cbind(out.tab.1, AIC.weight = out.tab.1$rel.like/sum(out.tab.1$rel.like))

Let’s look at the model weights (out.tab.1):

               H delta.AICc AIC.weight converged
 NC+strait+PS+SC       0.00      0.979      TRUE
            site       7.65      0.021      TRUE
             N+S      36.97      0.000      TRUE
           stock      47.02      0.000      TRUE
        coast+PS      48.78      0.000      TRUE
       panmictic      71.67      0.000      TRUE