Estimate of $\text{Var}(\hat\beta)$ in a linear mixed model

Question

Let $Y = X\beta + Zu + \sigma\epsilon$ be a Gaussian linear mixed model.

Let $V = Var(Y)$ be the marginal variance matrix of $Y$. Define the matrix $$ \Phi = {(X'V^{-1}X)}^{-1}. $$ According to this SAS documentation, $\Phi$ underestimates $\text{Var}(\hat\beta)$.

I've checked by simulations, and I don't find this result:

library(lme4)

# simulates a mixed 2-way ANOVA model ####
SimAV2mixed <- function(I, J, Kij, mu=0, alphai, sigmaO=1,
                        sigmaPO=1, sigmaE=1, factor.names=c("Part","Operator"),
                        resp.name="y", keep.intermediate=FALSE){
  Operator <- rep(1:J, each=I)
  Oj <- rep(rnorm(J, 0, sigmaO), each=I)
  Part <- rep(1:I, times=J)
  Pi <- rep(alphai, times=J)
  POij <- rnorm(I*J, 0, sigmaPO)
  simdata0 <- data.frame(Part, Operator, Pi, Oj, POij)
  simdata0$Operator <- factor(simdata0$Operator)
  levels(simdata0$Operator) <- 
    sprintf(paste0("%0", floor(log10(J))+1, "d"), 1:J)
  simdata0$Part <- factor(simdata0$Part)
  levels(simdata0$Part) <- sprintf(paste0("%0", floor(log10(I))+1, "d"), 1:I)
  simdata <- 
    as.data.frame(
      sapply(simdata0, function(v) rep(v, times=Kij), simplify=FALSE))
  Eijk <- rnorm(sum(Kij), 0, sigmaE)
  simdata <- cbind(simdata, Eijk)
  simdata[[resp.name]] <- mu + with(simdata, Oj+Pi+POij+Eijk)
  levels(simdata[,1]) <- paste0("A", levels(simdata[,1]))
  levels(simdata[,2]) <- paste0("B", levels(simdata[,2]))
  names(simdata)[1:2] <- factor.names
  if(!keep.intermediate) simdata <- simdata[,c(factor.names,resp.name)]
  simdata
}

# marginal variance matrix V of y ####
Vfun <- function(vcs, sigma2, reTrms){ # vcs = variance components / sigma2
  G <- Diagonal(x = rep(vcs, times = diff(reTrms$Gp)))
  Zt <- reTrms$Zt
  H <- Diagonal(ncol(Zt)) + t(Zt) %*% G %*% Zt
  sigma2 * H
}

# calculate the Phi matrix ####
Kij <- c(3, 3, 3, 4, 4, 4)
dat <- 
  SimAV2mixed(I=2, J=3, Kij = Kij, mu = 0, alphai = c(3,-3), 
              sigmaO = 3, sigmaPO = 2)
reTrms <- 
  lFormula(y ~ Part + (1|Operator) + (1|Operator:Part), data = dat)$reTrms
V <- Vfun(vcs = c(4,9), sigma2 = 1, reTrms)
Vinv <- chol2inv(chol(V))
X <- model.matrix(~ Part, data = dat)
( Phi <- solve(t(X) %*% Vinv %*% X) )
#           [,1]      [,2]
# [1,]  4.435113 -1.435089
# [2,] -1.435089  2.860920

# simulations ####
set.seed(666)
nsims <- 1000L
betahat <- matrix(NA_real_, nrow = nsims, ncol = 2L)
k <- 0L
while(k < nsims){
  dat <- 
    SimAV2mixed(I=2, J=3, Kij = Kij, mu = 0, alphai = c(3,-3), 
                sigmaO = 3, sigmaPO = 2)
  fit <- 
    suppressMessages(
      lmer(y ~ Part + (1|Operator) + (1|Operator:Part), data = dat, 
           start = setNames(c(2,3), c("Operator:Part", "Operator")), 
           control = lmerControl(optimizer = "Nelder_Mead", 
                                 optCtrl=list(maxfun=100000)), 
           REML = TRUE)
    )
  if(!is.null(fit@optinfo$conv$lme4$code)){ # skip if non-convergence
    print(k)
    next
  }
  k <- k + 1L
  betahat[k,] <- fit@beta
}

cov(betahat)
#           [,1]      [,2]
# [1,]  4.194572 -1.332589
# [2,] -1.332589  3.011642
Phi
#           [,1]      [,2]
# [1,]  4.435113 -1.435089
# [2,] -1.435089  2.860920

As we can see, $\Phi_{1,1}$ is higher than $\text{Var}(\hat\beta)_{1,1}$.

Am I doing something wrong, or am I misunderstanding something?

Answer 1

The point in that the estimator for the marginal covariance matrix of $Y$, i.e., $\hat \Phi = {(X' \hat V^{-1}X)}^{-1}$ will be a biased estimator of $\Phi = {(X'V^{-1}X)}^{-1}$. If I see correctly in your simulation, you compare $\Phi$ with the empirically determined variance of $\hat \beta$, not $\hat \Phi$.

Nonetheless, this bias is not that great. Check also my version below (in the first simulation I use 30 subjects in which case you would expect more bias in the estimation of $\Phi$ than with 200 subjects):

simfun <- function (n = 100, K = 3) {
    # 'n' number of subjects
    # 'K' number of measurements per subject
    t_max <- 15 # maximum follow-up time

    # we constuct a data frame with the design: 
    # everyone has a baseline measurment, 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, 1, 0.9) # fixed effects coefficients
    D11 <- 4 # variance of random intercepts
    D22 <- 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 <- drop(X %*% betas + rowSums(Z * b[DF$id, ]))
    # we simulate normal longitudinal data
    DF$y <- rnorm(n * K, eta_y, sd = 1)
    DF
}

library("lme4")

M <- 1000
betas <- var_betas <- matrix(NA_real_, M, 3)
for (m in seq_len(M)) {
    DF_m <- simfun(n = 30)
    fm_m <- suppressWarnings(lmer(y ~ sex + time + (time | id), data = DF_m))
    betas[m, ] <- fixef(fm_m)
    var_betas[m, ] <- diag(vcov(fm_m))
}

diag(var(betas))
colMeans(var_betas)

#########

M <- 1000
betas <- var_betas <- matrix(NA_real_, M, 3)
for (m in seq_len(M)) {
    DF_m <- simfun(n = 200)
    fm_m <- suppressWarnings(lmer(y ~ sex + time + (time | id), data = DF_m))
    betas[m, ] <- fixef(fm_m)
    var_betas[m, ] <- diag(vcov(fm_m))
}

diag(var(betas))
colMeans(var_betas)