4
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – user158565 Dec 31 '19 at 17:32
  • $\begingroup$ @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 ? $\endgroup$ – Stéphane Laurent Dec 31 '19 at 17:38
  • $\begingroup$ I withdraw from it. $\endgroup$ – user158565 Dec 31 '19 at 21:42
2
$\begingroup$

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)
$\endgroup$
  • 1
    $\begingroup$ 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? $\endgroup$ – Stéphane Laurent Jan 2 at 8:32
  • $\begingroup$ @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)$. $\endgroup$ – Dimitris Rizopoulos Jan 2 at 17:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.