# Linear mixed model with heteroscedastic residuals
mod <- galamm(
formula = y ~ x + (1 | id),
weights = ~ (1 | item),
data = hsced
)
# Extract covariance matrix for fixed regression coefficients
vcov(mod, parm = "beta")
# and then for weights, which gives us the variance.
vcov(mod, parm = "weights")
Run the code above in your browser using DataLab