Блог

Я пытаюсь добавить 95% доверительные интервалы к анализу DFA MARSS. Мой код

 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 amp; 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()
  

Я ожидаю, что значения загрузки необработанного фактора будут относительно близки к середине 95% ДИ, но результирующие столбики ошибок не совпадают с точками.

Вопрос 1: Корректны ли вращение varimax и матричная математика для определения загрузок коэффициентов DFA в MARSS?
Вопрос 2: Есть ли более простой (и правильный!!) как добавить доверительные интервалы к факторным загрузкам?

Вы хотите повернуть верхнюю и нижнюю Z-матрицы. К сожалению, ваш вопрос привел к обнаружению ошибки в coef() функции, которая затрудняет их получение. Но этот код является взломом, позволяющим обойти это. Он использует внутреннюю функцию, которую coef () использует для получения матриц параметров.

 # 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))