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

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?

• Need to understand the difference between the parameters (or function of the parameters) and estimators of them. $\Phi$ and $\text{Var}(\hat\beta)$ are functions of parameters. Commented Dec 31, 2019 at 17:32
• @user158565 I think I totally understand that. According to the SAS documentation I linked, $\Phi$ (the "true" $\Phi$, which is not a random variable) underestimates $\text{Var}(\hat\beta)$. What is wrong with my simulations ? Commented Dec 31, 2019 at 17:38
• I withdraw from it. Commented Dec 31, 2019 at 21:42

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)

• Thanks, interesting answer (+1). But this is not what I meant. The SAS documentation I linked claims: "even $\Phi$ (not $\hat\Phi$) underestimates $\text{Var}(\hat\beta)$ when $\sigma$ is unknown". This is the claim I want to check. Maybe I misunderstand it? Commented Jan 2, 2020 at 8:32
• @StephaneLaurent It is not evident to me how you could validate this statement because if $V$ is unknown you cannot calculate $\Phi$. As far as I can see, in your simulation you calculate $\Phi$ with known $V$ and you compare it with $\mbox{var}(\hat \beta)$. Commented Jan 2, 2020 at 17:37