Standard Error of Variance Component from the output of GLMMadaptive::mixed_model

Question

I am using the {GLMMadaptive} package to fit a mixed effect random slope model. And I need to extract the standard error of variance components from the output of GLMMadaptive::mixed_model(). And from the package documentation, it seems that I can use the vcov() method to extract the variance of the random components. But I am confused by the returned values.

Consider the following reprex from the package article Methods for MixMod Objects

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

Now the estimated variance-covariance matrix of the maximum likelihood estimates of random effects is returned using the vcov() method,

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

Then that package article says about the returned value of vcov,

The elements of this covariance matrix that correspond to the elements of the covariance matrix of the random effects (i.e., the elements D_xx) are on the log-Cholesky scale.

I am really confused by the above line. What does it mean by The elements of this covariance matrix are on the log-Cholesky scale.?

Please Note that my end goal is to get the standard errors of the estimated random effects that is, $Se(D_{11})$, $Se(D_{12})$, $Se(D_{22})$. So do I need to apply any transformation on the resulting matrix to get these? If so, how?

---

Note that: I am aware of this Q/A thread from stackoverflow that contains very useful discussion about this, but that uses {lme4} package. But my issue is specifically about the {GLMMadaptive} package.

Answer 1

You could use the delta method to calculate standard errors (of the estimators of the entries of the random effects' covariance matrix) on the variance-covariance scale.

Note, however, that there are better alternatives to quantify uncertainty (e.g., in terms of differently constructed confidence intervals) of such estimators.

In the code below, the function D_chol_to_D (corresponding to $g$ in the linked StatLect article) is the transformation from the log-Cholesky parameterization of the entries of the random effects' covariance matrix to the original variances and covariances. The function numDeriv::jacobian computes a numerical approximation to the value of the Jacobian matrix of $g$ (denoted by $J_g$ in the StatLect article) at the MLE D_chol_entries of the log-Cholesky parameterized random effects' covariance matrix.

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