Why is REML default if it inflates MSE?

Question

Within the mixed effects model world, REML has become the method of choice in order to correct for the downward bias in variance components. For years, I accepted this rationale without thinking about the potential effects this bias-correction might have on the estimator's variability.

Recently, however, I bumped into the article "On the inefficiency of the restricted maximum likelihood" (Longford, 2015), in which it is shown that "this unbiasedness is accompanied in some balanced designs by an inflation of the mean squared error". The author therefore aims to "encourage reevaluation of the uncritical preference for REML".

This makes me wonder; Why is the REML method of choice when it in fact increases the MSE, which is supposed to reflect the general quality of the estimator (i.e., low bias and low variability). Unbiasedness is a useful property, but in combination with large variability our estimate could still be located very far from the population value. From an applied researcher's perspective, wouldn't we prefer a slightly biased estimator, with lower variability?

I'd be happy with any intuitive explanations, as well as more detailed references to articles or books. Any reading material I've come across usually highlights the pros (i.e. how REML reduces bias), but the cons (e.g. MSE/variability) are hardly ever mentioned.

Thank you in advance!

Answer 1

I'll have to be honest that I have not previously heard that REML results in larger mean squared error "in some balanced designs". I don#t have access to the article, but I thought I would run some simulations to investigate.

Here we simulate data for a mixed model with a balanced design: 20 groups A, 5 observations per group, with 2 fixed effects, X and Z. We do this 100 times:

dt <- expand.grid(A = 1:20, reps = 1:5)
betas <- c(10, 2, 4)

n_sim <- 100
vec_mse_reml <- numeric(n_sim)
vec_mse_ml <- numeric(n_sim)

for (i in 1:n_sim) {
  set.seed(i)
  dt$X <- rnorm(nrow(dt))
  dt$Z <- rnorm(nrow(dt))
  X <- model.matrix(~ X + Z, dt)      # Model matrix for fixed effects

  dt$Y <- 1    # just used to obtain Z in the next 3 lines
  myFormula <- "Y ~ X + Z + (1 | A)"
  foo <- lFormula(eval(myFormula), dt)
  Z <- t(as.matrix(foo$reTrms$Zt))  # model matrix for random effects

  b <- rnorm(20, 0, 2)      # random intercepts

  # now simulate Y using the general mixed model equation :
  dt$Y <- X %*% betas + Z %*% b + rnorm(nrow(dt), 0, 5)

  m0_reml <- lmer(myFormula, dt, REML = TRUE)
  m0_ml <- lmer(myFormula, dt, REML = FALSE)

  vec_mse_reml[i] <- mean(dt$Y - predict(m0_reml, dt)^2)
  vec_mse_ml[i] <- mean(dt$Y - predict(m0_ml, dt)^2)
}

And now we look at the overall mean and standard deviations:

> mean(vec_mse_reml); sd(vec_mse_reml)
[1] -110.8568
[1] 14.99271
> mean(vec_mse_ml);sd(vec_mse_ml)
[1] -110.7117
[1] 14.97813

So, with these simulations there is a tiny difference in both the mse and the sdard deviation of it. I would contend that in this case the difference is so small to be irrelevant. It would be interesting to know what models Longford was referring to.