Compare DFA models with spatial structure

[slarge]

I want to identify common trends in environmental conditions in estuaries across the Northeast US. Each estuary has about 4 stations where conditions are measured monthly. These stations are very close to each other and I want to test different hypotheses of spatial structure (similar to ATSA chapter 7.7. I think that I have set up the Z matrices in a way that could test this, but I am not very confident! I tried to create an example from the Lake WA data. From the code (below), the "two distinct sites" hypothesis seems to be supported. Is there a better way to nest subpopulations within a DFA framework? Thanks in advance for any feedback.

` library(MARSS) library(dplyr) library(tidyr)

load the data (there are 3 datasets contained here)

data(lakeWAplankton, package = "MARSS")

we want lakeWAplanktonTrans, which has been transformed

so the 0s are replaced with NAs and the data z-scored

all_dat <- lakeWAplanktonTrans %>% as_tibble() %>% filter(Year >= 1980, Year <= 1989) %>% mutate(date = as.Date(paste0(Year, "-", Month, "-1"))) %>% select(date, Cryptomonas, Diatoms, Greens, Unicells, Other.algae) %>% pivot_longer(cols = c("Cryptomonas", "Diatoms", "Greens", "Unicells", "Other.algae"))

Create a new variable "site" with some additional error

site_dat <- bind_rows(all_dat %>% mutate(site = "A"), all_dat %>% group_by(name) %>% rowwise() %>% mutate(site = "B", value = value + rnorm(1, mean = 0, sd = 1)))

site_wide <- site_dat %>% mutate(name_site = paste0(name, "_", site)) %>% select(-name, -site) %>% pivot_wider(names_from = date, values_from = value)

site_names <- site_wide[[1]] site_wide <- as.matrix(site_wide[, -1]) row.names(site_wide) <- site_names N_ts <- dim(site_wide)[1]

H1 two distinct sites

Z_vals <- list("z11", 0, 0, "z11", "z12", 0, "z11", "z12", "z13", "z11", "z12", "z13", "z11", "z12", "z13", "z21", "z22", "z23", "z21", "z22", "z23", "z21", "z22", "z22", "z21", "z22", "z23", "z21", "z22", "z23") Z1 <- matrix(Z_vals, nrow = N_ts, ncol = 3, byrow = TRUE)

H2 one site

Z_vals <- list("z011", 0, 0, "z021", "z022", 0, "z031", "z032", "z033", "z041", "z042", "z043", "z051", "z052", "z053", "z061", "z062", "z063", "z071", "z072", "z073", "z081", "z082", "z083", "z091", "z092", "z093", "z101", "z102", "z103") Z2 <- matrix(Z_vals, nrow = N_ts, ncol = 3, byrow = TRUE)

aa <- "zero" ## 'aa' is the offset/scaling

DD <- "zero" # matrix(0,mm,1) ## 'DD' and 'd' are for covariates dd <- "zero" # matrix(0,1,wk_last) ## 'RR' is var-cov matrix for obs errors RR <- "diagonal and unequal" mm <- 3 ## number of processes BB <- "identity" # diag(mm) ## 'BB' is identity: 1's along the diagonal & 0's elsewhere uu <- "zero" # matrix(0, mm, 1) ## 'uu' is a column vector of 0's CC <- "zero" # matrix(0, mm, 1) ## 'CC' and 'cc' are for covariates cc <- "zero" # matrix(0, 1, wk_last) QQ <- "identity" # diag(mm) ## 'QQ' is identity

list with specifications for model vectors/matrices

mod_list <- list(Z = Z1, A = aa, D = DD, d = dd, R = RR, B = BB, U = uu, C = CC, c = cc, Q = QQ) init_list <- list(x0 = matrix(rep(0, mm), mm, 1)) con_list <- list(maxit = 3000, allow.degen = TRUE)

fit MARSS H1

dfa_1 <- MARSS(y = site_wide, model = mod_list, inits = init_list, form = "dfa", control = con_list)

fit MARSS H2

mod_list$Z <- Z2 dfa_2 <- MARSS(y = site_wide, model = mod_list, inits = init_list, control = con_list)

dfa_1$AICc ## 3308.385 dfa_2$AICc ## 3320.462`