Skip to main content
Notice removed Draw attention by Shafee
Bounty Ended with statmerkur's answer chosen by Shafee
edited tags
Link
Notice added Draw attention by Shafee
Bounty Started worth 50 reputation by Shafee
Source Link

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

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.