MARSS model with multiple species and multiple sites

Question

I've been trying to figure out how to properly format a MARSS model (using the MARSS package in R) so that I can model multiple species at multiple sites simultaneously, if possible. In my full study design, I have yearly samples of 9 species from 9 different sites across 30 years. I've been able to get the models to run for one site with multiple species and for multiple sites and only one species (see reprex below, though I'm not entirely sure those are coded correctly), but haven't been able to figure out how to put it all together or if that's even something that can be done. I could also just run separate models for each site if that's an appropriate approach. Ideally I'd like to get separate B-matrices for each of my sites so that I can compare community stability metrics between sites. If anyone has any thoughts (or if I'm going about this wrong) I'd really appreciate it!

library(MARSS)

dat <- data.frame(site = rep(1:2,each = 20),          # two independent sites
                  year = rep(c(1:20),2),              # sample years (1-20)
                  spp1 = sample(50,40, replace = T),  # simulated count data species 1
                  spp2 = sample(20,40, replace = T),  # simulated count data species 2
                  spp3 = sample(30,40, replace = T))  # simulated count data species 3

# Put count data into wide format (rows = species, columns = years)
site1 <- t(dat[1:20,3:5])
site2 <- t(dat[21:40,3:5])

rownames(site1) <- paste0("site1_", rownames(site1))
rownames(site2) <- paste0("site2_", rownames(site2))

y <- rbind(site1,site2)

# RUN THE MODEL FOR ONE SPECIES AND BOTH SITES
y1 <- rbind(y[1,], y[4,])

# Z-matrix formatting
z.model1 <- matrix(c(1,0,0,1),ncol = 2)

# A-matrix formatting
a.model1 <- matrix(list(0),2,1)
a.model1[2,1] <- c("b")

# U-matrix formatting
u.model1 <- matrix(list(0),2,1)
u.model1[1:2,1] <- c("a","b")

# Put together all parts of model
model.1 <- list(Z = z.model1,                
                A = a.model1,                
                Q = "diagonal and unequal", 
                R = "diagonal and equal",                     
                U = u.model1,                
                tinitx = 1) 

# Run the model
m1.out <- MARSS(y1, model = model.1, method = "BFGS")


# NOW RUN THE MODEL WITH ALL THREE SPECIES BUT ONLY ONE SITE
y2 <- y[1:3,]

# Z-matrix formatting
z.model2 <- matrix(0,3,3)
diag(z.model2) <- 1

# U-matrix formatting
u.model2 <- matrix(list(0),3,1)
u.model2[1:3,1] <- c("a","b","c")

# B-matrix formatting
b.model2 <- matrix(c("b11","b12","b13",
                     "b21","b22","b23",
                     "b31","b32","b33"),
                   ncol = 3)

# Put together all parts of model
model.2 <- list(Z = z.model2,                
                B = b.model2,
                Q = "diagonal and unequal",  
                R = "diagonal and equal",                     
                U = u.model2,                
                tinitx = 1) 

# Run the model
m2.out <- MARSS(y2, model = model.2, method = "BFGS")

Answer 1

I may have just figured out my question (see code below). I'll leave this question up in case anyone knows if I did this incorrectly.

# Using the same data as in the original question
y <- rbind(site1,site2)
y

# Create needed matrices (doing this manually for clarity)
# B-matrix formatting
b.model <- matrix(list(0),6,6)

b.model[1:3,1:3] <- c("A.b11","A.b12","A.b13",   # Each site has a 3x3 "matrix" embedded within
                      "A.b21","A.b22","A.b23",   # the overall 6x6 matrix. Not sure this is an
                      "A.b31","A.b32","A.b33")   # appropriate way of doing this, the plan would  
                                                 # be to isolate out the individual B-matrices      
b.model[4:6,4:6] <- c("B.b11","B.b12","B.b13",   # for each site from the "overall B-matrix"
                      "B.b21","B.b22","B.b23",
                      "B.b31","B.b32","B.b33")

# Z-matrix formatting
z.model <- matrix(0,6,6)
diag(z.model) <- 1        # same as "identity"
z.model

# Q-matrix formatting
q.model <- matrix(list(0),6,6)
diag(q.model) <- c("Qsp1A","Qsp2A","Qsp3A",
                   "Qsp1B","Qsp2B","Qsp3B")   # same as "diagonal & unequal"
q.model

# R-matrix formatting
r.model <- matrix(0,6,6)   # same as "zero" (for now considering detection as being perfect)

# U-matrix formatting
u.model <- matrix(c("Usp1A","Usp2A","Usp3A",   # same as "unconstrained"
                    "Usp1B","Usp2B","Usp3B"),
                  6,1) 

# Put together all parts of model
model.0 <- list(Z = z.model,                
                B = b.model,
                Q = q.model,  
                R = r.model,                     
                U = u.model,                
                tinitx = 0) 

# Run the model
m0.out <- MARSS(y, model = model.0, method = "BFGS")

m0.B <- coef(m0.out, type = "matrix")$B[1:6, 1:6]
m0.B

# Separating out individual site's B-matrices
# Again, not sure this is appropriate/ok to do...
m0.B.siteA <- m0.B[1:3,1:3]
m0.B.siteB <- m0.B[4:6,4:6]

max(abs(eigen(m0.B.siteA)$values))
max(abs(eigen(m0.B.siteB)$values))