GLMMadaptive::mixed_model的输出中的方差成分标准误差

Question

我正在使用{GLMMadaptive}软件包来拟合混合效应随机斜率模型。而且我需要从GLMMadaptive::mixed_model()的输出中提取方差分量的标准误差。根据软件包文档, 似乎可以使用vcov()方法来提取随机分量的方差。但是我对返回的值感到困惑。

考虑以下来自软件包文章MixMod对象的方法的示例。

library(GLMMadaptive)

set.seed(1234)
n <- 100 # number of subjects
K <- 8 # number of measurements per subject
t_max <- 15 # maximum follow-up time

# we construct a data frame with the design: 
# everyone has a baseline measurement, and then measurements at random follow-up times
DF <- data.frame(id = rep(seq_len(n), each = K),
                 time = c(replicate(n, c(0, sort(runif(K - 1, 0, t_max))))),
                 sex = rep(gl(2, n/2, labels = c("male", "female")), each = K))

# design matrices for the fixed and random effects
X <- model.matrix(~ sex * time, data = DF)
Z <- model.matrix(~ time, data = DF)

betas <- c(-2.13, -0.25, 0.24, -0.05) # fixed effects coefficients
D11 <- 0.48 # variance of random intercepts
D22 <- 0.1 # variance of random slopes

# we simulate random effects
b <- cbind(rnorm(n, sd = sqrt(D11)), rnorm(n, sd = sqrt(D22)))
# linear predictor
eta_y <- as.vector(X %*% betas + rowSums(Z * b[DF$id, ]))
# we simulate binary longitudinal data
DF$y <- rbinom(n * K, 1, plogis(eta_y))

fm <- mixed_model(fixed = y ~ sex * time, random = ~ time | id, data = DF,
                  family = binomial())

现在使用vcov()方法返回随机效应的最大似然估计的方差-协方差矩阵。

vcov(fm, parm = "var-cov")

#>             D_11        D_12       D_22
#> D_11  0.42942062 -0.09963969  0.5065884
#> D_12 -0.09963969  0.03701847 -0.2117451
#> D_22  0.50658839 -0.21174511  1.3651870

那个软件包文章中提到了vcov的返回值，

与随机效应协方差矩阵对应的矩阵元素（即D_xx元素）以对数Cholesky尺度表示。

我对上述句子感到非常困惑。"The elements of this covariance matrix are on the log-Cholesky scale."是什么意思？

请注意，我的最终目标是获得估计的随机效应的标准误差，即SE(D_11)、SE(D_12)、SE(D_22)。那么，我需要对结果矩阵进行任何转换才能获得这些值吗？如果需要转换，应该如何进行？

---

请注意：我知道这个问答主题包含了关于这个问题非常有用的讨论，但是那个使用了{lme4}软件包。而我的问题特别涉及到{GLMMadaptive}软件包。

- shafee

---

能不能解释一下这些标准误差在这里有什么用呢？方差的不确定性分布并不对称。使用confint(fm, "var-cov")不是更有用吗？ - undefined

@Roland，是的，我确实知道这一点来自这里。但实际上我确实需要标准误差（即使它们不是很可靠的不确定性指标）:) - undefined

Answer 1

下面的答案是 我的回答 到 Cross Validated 上的转贴问题 的副本。

您可以使用 Delta 方法 来计算方差协方差尺度上的标准误差（随机效应协方差矩阵估计值的估计量）。

然而，请注意，有更好的选择来量化这些估计值的不确定性（例如，以不同构建的置信区间为单位）。

在下面的代码中，函数 D_chol_to_D（对应于链接的 StatLect 文章中的 g）是将随机效应协方差矩阵的对数Cholesky参数化转换为原始方差和协方差的函数。函数 numDeriv::jacobian 计算 g 的 Jacobian 矩阵的数值近似值（在 StatLect 文章中用 J_g 表示），该值基于对数Cholesky参数化随机效应协方差矩阵的 MLE D_chol_entries。

library("numDeriv")

# estimated covariance matrix of random effects
D <- fm$D

# transform from covariance matrix to entries of cholesky factor with 
# log-transformed main diagonal
D_chol_entries <- GLMMadaptive:::chol_transf(D)

D_chol_to_D <- function(x) {
  
  # transform from entries of cholesky factor with log-transformed main diagonal
  # to covariance matrix
  D <- GLMMadaptive:::chol_transf(x)
  
  D[upper.tri(D, diag = TRUE)]
}

J <- jacobian(D_chol_to_D, D_chol_entries)

# estimated covariance matrix of D_chol_entries
V_chol <- vcov(fm, parm = "var-cov")

# estimated covariance matrix of entries of D
V <- J %*% V_chol %*% t(J)

se <- sqrt(diag(V))

cat("--- D_11 ---\nEstimate:", D[1, 1], "\nStd. Error:", se[1], fill = TRUE)
# --- D_11 ---
# Estimate: 0.4639974 
# Std. Error: 0.5992368
cat("--- D_22 ---\nEstimate:", D[2, 2], "\nStd. Error:", se[3], fill = TRUE)
# --- D_22 ---
# Estimate: 0.06304059 
# Std. Error: 0.02843676

- statmerkur