rstate-space

How to add confidence intervals to MARSS package DFA factor loadings

I'm trying to add 95% confidence intervals to MARSS DFA analysis. My code

library(MARSS) ## version 3.10.12
library(dplyr)
library(tidyr)
library(ggplot2)

data(lakeWAplankton, package = "MARSS")

all_dat <- lakeWAplanktonTrans
yr_frst <- 1980
yr_last <- 1989
plank_dat <- all_dat[all_dat[, "Year"] >= yr_frst & all_dat[, "Year"] <= yr_last, ]

phytoplankton <- c("Cryptomonas", "Diatoms", "Greens", "Unicells", 
                   "Other.algae")
## get only the phytoplankton
dat_1980 <- plank_dat[, phytoplankton]

## transpose data so time goes across columns
dat <- t(dat_1980)

mod_list = list(m = 3, R = "diagonal and unequal")
con_list <- list(maxit = 3000, allow.degen = TRUE)
dfa_temp <- MARSS(dat, model = mod_list, form = "dfa", z.score = FALSE, 
                  control = con_list)

dfa_temp <- MARSSparamCIs(dfa_temp)


up = coef(dfa_temp, what = "par.upCI")$Z
lo = coef(dfa_temp, what = "par.lowCI")$Z
raw = coef(dfa_temp, what = "par")$Z
ci_list <- list(up = up, lo = lo, raw = raw)

f = dfa_temp$marss$fixed$Z
d = dfa_temp$marss$free$Z

dims = attr(dfa_temp$marss, "model.dims")$Z
par_names <- rownames(dfa_temp$call$data)

delem = d
attr(delem, "dim") = attr(delem, "dim")[1:2]

felem = f
attr(felem, "dim") = attr(felem, "dim")[1:2]

par_mat <- lapply(1:3, function(x) matrix(felem + delem %*% ci_list[[x]],
                                          dims[1], dims[2]))
names(par_mat) <- names(ci_list)

h_inv = varimax(par_mat$raw)$rotmat
h_inv_lo = varimax(par_mat$lo)$rotmat
h_inv_up = varimax(par_mat$up)$rotmat
z_rot = data.frame(var = "raw", par_mat$raw %*% h_inv)
z_rot_lo = data.frame(var = "lo", par_mat$lo %*% h_inv_lo)
z_rot_up = data.frame(var = "up", par_mat$up %*% h_inv_up)

loads <- rbind(z_rot, z_rot_lo, z_rot_up)
colnames(loads)[2:4] <- c("trend_1", "trend_2", "trend_3")
loads$plank <- rep(row.names(dat), 3)

loads_long <- pivot_longer(loads, cols = c("trend_1", "trend_2", "trend_3"),
                                  names_to = "trend", values_to = "val")

loads_dat <- pivot_wider(loads_long, names_from = var, values_from = val)

ggplot(data = loads_dat, 
       aes(x = plank, 
           ymin = lo,
           ymax = up,
           y = raw)) +
  geom_hline(yintercept = 0, color = "black") +
  geom_pointrange() +
  facet_wrap(~trend) +
  coord_flip() +
  theme_bw()

DFA Factor loadings plot

I expect the raw factor loading values to be relatively close to the middle of the 95% CI, but the resulting error bars do not align with the points.

Question 1: Is the varimax rotation and matrix math correct to come up with MARSS DFA factor loadings?
Question 2: Is there an easier (and correct!!) way to add confidence intervals to factor loadings?

Solution

  • You want to rotate the upper and lower Z matrices. Unfortunately, your question lead to a discovery of a bug in the coef() function that makes it hard to get those. But this code is a hack to get around that. It uses the internal function that coef() uses to get the parameter matrices.

    # add the par.lowCI and par.upCI to the fit
dfa_temp <- MARSSparamCIs(dfa_temp)

# get the rotation matrix for the Z
Z <- coef(dfa_temp, type = "matrix")$Z
H.inv <- varimax(Z)$rotmat

# Get the upZ, lowZ
# normally this would work coef(dfa_temp, type="matrix", what="par.lowCI")
# but there is a bug and it is not respecting type="matrix"
# I use the internal function parmat() which is using the $par element
#  just need to replace $par with the upper and lower pars
tmp <- dfa_temp; tmp$par <- tmp$par.upCI
Z.up <- MARSS:::parmat(tmp)$Z
tmp <- dfa_temp; tmp$par <- tmp$par.lowCI
Z.low <- MARSS:::parmat(tmp)$Z

Z.rot <- Z %*% H.inv
Z.rot.up <- Z.up %*% H.inv
Z.rot.low <- Z.low %*% H.inv

df <- data.frame(Z=as.vector(Z.rot), Zup=as.vector(Z.rot.up), Zlow=as.vector(Z.rot.low))