Cluster-robust SEs in State Space Model

Hi all,

I am using the State Space Model with multiple subjects and would like to compute cluster-robust standard errors to account for the multi-level structure. I am aware of the function imxRobustSE, but it says it's "for those people who know what they're doing" and I don't think that applies to me. Moreover, the function "does not work with multigroup models in which the groups themselves contain subgroups" and "does not correctly handle multilevel data". Am I correct in assuming that this applies to multi-subject SSM and I thus cannot compute cluster-robust SEs?

Below the code of the model I am using:

# A matrix (= phi, regression effects)
amat <- mxMatrix(type = "Full", nrow = 2, ncol = 2,
                 free = TRUE,
                 values = c(.5, .5,
                            .5, .5),
                 name = "A",
                 labels = c("phi11", "phi12",
                            "phi21", "phi22"),
                 byrow = TRUE)

# B matrix (= covariates on latent variables, in this case 0s)
bmat <- mxMatrix("Zero", nrow = 2, ncol = 1, 
                 free = FALSE,
                 name = "B")

# C matrix (= lambda, factor loadings)
clab <- matrix(c(paste0("l", 1:8, "_1"), paste0("l", 1:8, "_2")), ncol = 2)
cdim <- list(paste0("v", 1:8), c("F1", "F2"))

cmat <- mxMatrix("Full", nrow = 8, ncol = 2,
                 free = c(FALSE, FALSE,
                          TRUE, FALSE,
                          TRUE, FALSE,
                          TRUE, FALSE,
                          FALSE, FALSE,
                          FALSE, TRUE,
                          FALSE, TRUE,
                          FALSE, TRUE),
                 values = c(1, 0,
                            .5, 0,
                            .5, 0,
                            .5, 0,
                            0, 1,
                            0, .5,
                            0, .5,
                            0, .5),
                 name = "C",
                 dimnames = cdim,
                 labels = clab,
                 byrow = TRUE)

# D matrix (covariates on observed variables, in this case 0s)
dmat <- mxMatrix("Zero", nrow = 8, ncol = 1,
                 name = "D")

# Q matirx (= zeta, innovation variances)
qmat <- mxMatrix("Symm", 2, 2,
                 free = TRUE,
                 values = c(1.5,
                            .5, 1.5),
                 name = "Q",
                 labels = c("zeta1", 
                            "zeta12", "zeta2"),
                 byrow = TRUE)

# R matrix (= theta, residual variances)
rlab <- paste0("resid", 1:8)
rmat <- mxMatrix("Diag", nrow = 8, ncol = 8,
                 free = TRUE,
                 .2, name = "R",
                 labels = rlab)

# x0 and P0
xmat <- mxMatrix("Zero", nrow = 2, ncol = 1, name = "x0")
pmat <- mxMatrix("Full", 2, 2, TRUE,
                 values = c(1, 1,
                            1, 1),
                 name = "P0")

umat <- mxMatrix("Zero", nrow = 1, ncol = 1, name = "u")

fullmodel <- mxModel(model = "LVAR",
                     amat, bmat, cmat, dmat, qmat, rmat, xmat,
                     pmat, umat,
                     mxData(observed = data[data$id ==1,], type = "raw"),
                     mxExpectationStateSpace("A", "B", "C",
                                             "D", "Q", "R", "x0", "P0", "u"),
                     mxFitFunctionML())

indivmodels = list()
modNames <- paste0("indiv", unique(data$id))
for(i in unique(data$id)){
  data_i <- data[data$id == i, 1:8]
  bmat_i <- bmat
  # bmat_i$labels <- paste0(bmat$labels, "_", i)
  indivmodels[[i]] <- mxModel(model = modNames[i],
                              amat, bmat_i, cmat, dmat, qmat, rmat, xmat,
                              pmat, umat,
                              mxData(observed = data[data$id ==1,], type = "raw"),
                              mxExpectationStateSpace("A", "B", "C",
                                                      "D", "Q", "R", "x0", "P0", "u"),
                              mxFitFunctionML(),
                              mxData(data_i, type = "raw"))
}
multiSubjModel <- mxModel(name = "multiMod", indivmodels, mxFitFunctionMultigroup(modNames))
multiSubjRun <- mxRun(multiSubjModel)
summary(multiSubjRun)