CS1.3 — REML vs ML Estimation

1 Before you start: concepts this chapter builds on

This chapter explains how the mixed model parameters (\(\boldsymbol{\beta}\), \(\sigma_b^2\), \(\sigma^2\)) are estimated using likelihood. Two background ideas are essential: what a likelihood is, and why dividing by \(n\) gives biased variance estimates. Expand each section below if the concept feels unfamiliar — otherwise skip ahead.

Note1. Likelihood and log-likelihood

What is a likelihood?

The likelihood \(L(\theta)\) is the probability of observing your data, viewed as a function of the parameters \(\theta\) rather than of the data. For a single normally distributed observation \(y\):

\[L(\mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\!\left(-\frac{(y - \mu)^2}{2\sigma^2}\right)\]

Maximum likelihood estimation (MLE) finds the parameter values that make the observed data most probable — the \(\theta\) that maximises \(L(\theta)\).

Why log? For \(n\) independent observations, \(L(\theta) = \prod_{i=1}^n f(y_i;\, \theta)\). Products of small probabilities cause numerical underflow and are hard to differentiate. The log converts the product to a sum without changing the location of the maximum:

\[\ell(\theta) = \log L(\theta) = \sum_{i=1}^n \log f(y_i;\, \theta)\]

The peak is exactly \(\bar{y}\) — this is why the MLE for \(\mu\) in a normal model is the sample mean.

Note2. Why dividing by n underestimates variance — and what this has to do with REML

The standard sample variance divides by \(n - 1\):

\[s^2 = \frac{\sum(y_i - \bar{y})^2}{n - 1}\]

But the MLE for \(\sigma^2\) divides by \(n\):

\[\hat{\sigma}^2_{\text{ML}} = \frac{\sum(y_i - \bar{y})^2}{n}\]

The ML version is biased downward: it plugs in the estimated \(\bar{y}\) without accounting for the 1 degree of freedom spent estimating it. Dividing by \(n - 1\) corrects for this.

The connection to REML: in a mixed model there are \(p\) fixed-effect parameters (not just 1 mean) consuming degrees of freedom. REML applies the same correction — it integrates out \(\boldsymbol{\beta}\) before estimating the variance components \(\sigma_b^2\) and \(\sigma^2\), recovering the lost degrees of freedom. ML does not, so it underestimates the variance components, especially when the number of clusters \(m\) is small.

Note3. The random intercept model — recap from CS1.2

The model fitted in this chapter:

\[Y_{ij} = \beta_0 + \beta_1 \cdot \text{treat}_{ij} + \beta_2 \cdot \text{dbp1}_{ij} + b_i + \varepsilon_{ij}\]

where \(b_i \sim \mathcal{N}(0,\, \sigma_b^2)\) and \(\varepsilon_{ij} \sim \mathcal{N}(0,\, \sigma^2)\).

Symbol	What it is	Estimated how?
\(\beta_0, \beta_1, \beta_2\)	Fixed effects — same for every centre	Directly by ML or REML
\(\sigma_b^2\)	Between-centre variance	From the likelihood — this is what ML underestimates
\(\sigma^2\)	Within-centre (residual) variance	From the likelihood
\(b_i\)	Random intercept for centre \(i\)	Predicted afterward (BLUP), not estimated

What this chapter is about: ML and REML produce the same \(\hat{\boldsymbol{\beta}}\) but differ in how they estimate \(\sigma_b^2\) and \(\sigma^2\). REML is more accurate.

---

2 Concept: How Are Mixed Model Parameters Estimated?

The random intercept model has three sets of unknowns that need to be estimated from data:

• \(\boldsymbol{\beta}\) — the fixed-effect coefficients (treatment effect, baseline association, …)
• \(\sigma_b^2\) — the between-cluster variance
• \(\sigma^2\) — the within-cluster (residual) variance

\(\sigma_b^2\) and \(\sigma^2\) together are called the variance components. Both ML and REML find the values that make the observed data most plausible — they differ only in how they handle the fixed effects while doing so.

The fixed effects β (treatment effect, baseline association, ...)

The between-cluster variance σ²_b

The individual random intercepts b_i for each cluster

The residual variance σ²

The b_i are random effects — they are not estimated as parameters but predicted afterward from the fitted model (as BLUPs). The three unknowns estimated by the likelihood are β, σ²_b, and σ². ML and REML both estimate these three, but differ in how they handle β while estimating the variance components.

2.1 ML: Maximum Likelihood

ML estimates \(\boldsymbol{\beta}\) and the variance components simultaneously by maximising the likelihood of the data. The problem: while fitting \(\boldsymbol{\beta}\), ML “uses up” degrees of freedom without accounting for this. The result is that ML systematically underestimates variance components — the same bias as dividing by \(n\) instead of \(n-1\) when computing a sample variance (see prerequisite section 2 above).

The bias decreases as the number of clusters \(m\) grows, but is substantial when \(m\) is small (< 30).

The number of observations per cluster

The degrees of freedom used up when estimating the fixed effects β

The correlation between residuals and random effects

The shape of the likelihood function near its peak

Estimating the p fixed effects β "uses up" p degrees of freedom, just like estimating the mean uses 1 degree of freedom in a simple variance calculation. ML treats β̂ as if it were the true β — it does not deduct these degrees of freedom from the variance estimation. With few clusters (small m) and many fixed effects (large p), the downward bias is most severe.

2.2 REML: Restricted Maximum Likelihood

REML fixes the underestimation by removing \(\boldsymbol{\beta}\) from the picture before estimating the variance components. Concretely: it maximises a version of the likelihood in which the fixed effects have been projected out, so the \(p\) degrees of freedom spent estimating \(\boldsymbol{\beta}\) are correctly accounted for. This yields unbiased variance component estimates.

The trade-off: because REML has projected out \(\boldsymbol{\beta}\), two REML likelihoods are only comparable if the fixed effects are identical. You cannot use a REML-based likelihood ratio test to compare models with different fixed effects.

REML becomes numerically unstable and prone to convergence failures whenever the fixed-effect structure differs between the compared models

The fixed-effect estimates β̂ are systematically biased under REML when the number of predictors in the model exceeds a certain threshold

Models with different fixed effects have projected out different subspaces — REML likelihoods from such models describe different restricted quantities

REML requires both models to have the exact same number of estimated fixed-effect parameters to produce comparable likelihood values

Think of it this way: REML answers the question "given that β is integrated out, what variance components fit best?" If two models have different fixed effects, they integrate out different things — the resulting likelihoods are not on the same scale. ML keeps β in the likelihood explicitly, so two ML likelihoods always describe the same quantity and can be compared via LRT.

ImportantWhen to use REML vs ML

Purpose	Use
Estimating variance components (\(\sigma_b^2\), \(\sigma^2\)) for reporting	REML (default in lme4)
Comparing models with different fixed effects via LRT	ML (REML=FALSE) — REML likelihoods are not comparable when fixed effects differ
Comparing models with same fixed effects, different random effects	REML is approximately correct
Final model after fixed-effects selection	Refit with REML

Standard workflow: (1) select fixed effects using ML-based LRT, (2) refit the chosen model with REML for final variance component estimates.

Compare the two REML fits using anova()

Refit both models with REML=FALSE (ML) and then compare with anova()

Compare AIC from the default REML fits — AIC is invariant to the estimation method

Fit both with REML and subtract log-likelihoods manually

REML likelihoods are not comparable when fixed effects differ: different fixed effects mean different quantities were integrated out of the likelihood. Always use ML (REML=FALSE) when comparing models with different fixed effect structures. Note: lme4 will warn you and auto-refit with ML if you try anova() on REML fits with different fixed effects.

2.3 AIC and BIC: alternatives to the LRT

When comparing models, two information criteria appear in anova() output and summarise fit vs complexity in a single number:

\[\text{AIC} = -2\ell + 2k \qquad \text{BIC} = -2\ell + k\log n\]

where \(k\) = number of estimated parameters, \(n\) = total observations. Lower is better for both.

• AIC (Akaike Information Criterion): penalises each extra parameter by 2
• BIC (Bayesian Information Criterion): penalises each extra parameter by \(\log n\) — more strongly than AIC when \(n > 8\), so it favours simpler models in large samples

Use cases: - Cross-check alongside the LRT: if AIC and the LRT agree, more confidence in the result - Comparing non-nested models (different predictor sets, not one a special case of the other): the LRT is not valid, but AIC/BIC differences still give guidance

WarningAIC/BIC require consistent estimation method

Compare AIC or BIC only across models fitted with the same method (both ML, or both REML). An AIC from a REML fit and an AIC from an ML fit are not on the same scale and cannot be compared.

2.4 Boundary testing — testing whether a variance component is zero

A common question: “Is the random intercept variance significantly different from zero?” This is tested with a likelihood ratio test comparing a model with \(\sigma_b^2 > 0\) against one with \(\sigma_b^2 = 0\) (i.e., ordinary linear regression).

The catch: under \(H_0\), the parameter \(\sigma_b^2 = 0\) sits on the boundary of the parameter space (variances cannot be negative). The standard chi-squared distribution for the LRT statistic is not correct here — it gives a p-value that is too conservative.

Practical fix: halve the p-value from the LRT. This is an approximation, but standard practice:

anova(fit_lm, fit_lmer)   # LRT statistic
# Halve the p-value for the boundary correction

Note

The boundary correction applies when testing \(H_0: \sigma_b^2 = 0\). For testing fixed effects (treatment, covariates), the standard chi-squared approximation is fine.

p = 0.08 — no correction needed for variance tests

p = 0.16 — double it, because the LRT is anti-conservative at boundaries

p = 0.04 — halve it, because the standard χ² approximation is too conservative here

Cannot determine without knowing the true σ²_b

When the null hypothesis places a parameter exactly on the boundary of its valid range (variances cannot be negative), the standard chi-squared p-value is approximately twice what it should be — it is too conservative. Halving gives p = 0.04, which would be significant at α = 0.05. This correction applies only to boundary tests of variance components (H₀: σ² = 0), not to tests of fixed effects.

2.5 Visualising the bias of ML

# Simulate many datasets, compare ML vs REML estimates of sigma_b.
# 20 groups × 4 obs per group: enough for stable sigma_b estimation
# while still showing ML's downward bias (which decreases as m grows).
set.seed(1)
n_sim <- 500
ml_est <- reml_est <- numeric(n_sim)
for(i in 1:n_sim){
  d <- data.frame(
    g = factor(rep(1:20, each = 4)),
    y = rep(rnorm(20, 0, 3), each = 4) + rnorm(80, 0, 2)
  )
  fit_ml   <- lmer(y ~ 1 + (1|g), data=d, REML=FALSE)
  fit_reml <- lmer(y ~ 1 + (1|g), data=d, REML=TRUE)
  ml_est[i]   <- as.data.frame(VarCorr(fit_ml))$sdcor[1]
  reml_est[i] <- as.data.frame(VarCorr(fit_reml))$sdcor[1]
}
cat("True sigma_b = 3\n")

True sigma_b = 3

cat("ML mean estimate:   ", round(mean(ml_est), 3), "\n")

ML mean estimate:    2.878 

cat("REML mean estimate: ", round(mean(reml_est), 3), "\n")

REML mean estimate:  2.962 

TipWhat to notice

ML systematically underestimates \(\sigma_b\) — its mean is below the true value of 3. REML is closer. The gap between them shrinks as the number of groups \(m\) increases: with large \(m\), the degrees-of-freedom loss from estimating \(\boldsymbol{\beta}\) becomes negligible relative to the total information. With \(m = 20\) groups the bias is moderate; try changing rep(1:20, each=4) to rep(1:5, each=4) (and rnorm(80,...) to rnorm(20,...)) to see it become severe.

ML and REML estimates would converge — smaller m makes estimation easier

The bias would disappear — fewer groups means fewer parameters to estimate

The gap between ML and REML would widen — ML's bias is worst when m is small

REML would also become biased — it needs at least 20 groups to correct properly

ML's underestimation of variance components is most severe when the number of groups m is small, because the fixed effects consume a large fraction of the available information. With m = 5 groups, you would also expect to see frequent singular fits (σ̂_b = 0), since 5 groups barely provide enough information to separate between-group from within-group variance.

---

3 Worked Example: Fixed Effects Model Comparison

We compare two fixed effects structures. For this, we must use ML.

# Fit two models with different fixed effects using ML
fit1_ml <- lmer(dbp ~ treat + dbp1 + (1|center), data=dat, REML=FALSE)
fit2_ml <- lmer(dbp ~ dbp1 + (1|center),          data=dat, REML=FALSE)

# Likelihood ratio test
anova(fit2_ml, fit1_ml)

Data: dat
Models:
fit2_ml: dbp ~ dbp1 + (1 | center)
fit1_ml: dbp ~ treat + dbp1 + (1 | center)
        npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)    
fit2_ml    4 465.08 475.34 -228.54    457.08                         
fit1_ml    5 419.64 432.46 -204.82    409.64 47.441  1  5.668e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TipReading the anova() output

Each row is one model. The key columns:

Column	Meaning
Df	Number of estimated parameters in this model
AIC	Lower = better fit after penalising complexity (penalty = 2 per parameter)
BIC	Like AIC but stronger penalty (\(\log n\) per parameter) — more conservative in large samples
logLik	Log-likelihood at the MLE — higher = better raw fit
Chisq	LRT statistic = \(2(\ell_{\text{larger}} - \ell_{\text{smaller}})\)
Chi Df	Degrees of freedom = difference in number of parameters
Pr(>Chisq)	p-value from \(\chi^2_{df}\) — small means the extra parameters significantly improve fit

A small p-value supports the more complex model. For AIC: differences < 2 are negligible; > 10 is strong evidence for the model with lower AIC.

         Df    AIC    BIC   logLik   Chisq Chi Df Pr(>Chisq)
fit2_ml   5  498.3  514.2  -244.1
fit1_ml   6  490.8  508.7  -239.4   9.50      1   0.00206 **

fit2_ml is better — it has lower df and simpler structure

Adding treat significantly improves fit: χ²(1) = 9.50, p = 0.002; prefer fit1_ml

The result is inconclusive — AIC difference of 7.5 is too small to matter

We cannot conclude anything — REML fits are required for a valid LRT

The LRT statistic is 9.50 on 1 df (the 6 − 5 = 1 extra parameter in fit1_ml). p = 0.002 is well below 0.05 → adding treat significantly improves the model. AIC also decreases by 7.5 (498.3 − 490.8), further supporting fit1_ml. This comparison is valid because both models were fitted with ML (REML=FALSE). After selecting fit1_ml, the next step is to refit with REML for final estimates.

# After choosing the fixed effects structure, refit with REML for final estimates
fit_final <- lmer(dbp ~ treat + dbp1 + (1|center), data=dat, REML=TRUE)
summary(fit_final)

Linear mixed model fit by REML ['lmerMod']
Formula: dbp ~ treat + dbp1 + (1 | center)
   Data: dat

REML criterion at convergence: 414

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.0144 -0.5577 -0.0356  0.5597  3.4241 

Random effects:
 Groups   Name        Variance Std.Dev.
 center   (Intercept) 6.435    2.537   
 Residual             3.304    1.818   
Number of obs: 96, groups:  center, 8

Fixed effects:
               Estimate Std. Error t value
(Intercept)    82.88300    2.19974  37.678
treatTreatment -2.91250    0.37162  -7.837
dbp1            0.38181    0.02218  17.216

Correlation of Fixed Effects:
            (Intr) trtTrt
treatTrtmnt -0.035       
dbp1        -0.905 -0.055

The two-step procedure — (1) select fixed effects with ML, (2) estimate with REML — is standard practice.

Keep the ML fit — the fixed effect estimates are the same under ML and REML

Refit with REML to obtain unbiased variance component estimates for reporting

Refit with REML only if the ICC exceeds 0.1

Use BIC to choose between the ML and REML fits

The standard two-step procedure: (1) select fixed effects using ML-based LRT, (2) refit the chosen model with REML for final estimation. The REML variance components are less biased and should be what you report. Fixed effect estimates are nearly identical between ML and REML in large samples.

---

4 Guided Practice

R is loading in the background — the exercises below will be interactive once the status bar disappears. You can read the concept sections above immediately.

TipHints

• lmer(y ~ x + group + (1|site), data=d2, REML=FALSE)
• anova(fit_red_ml, fit_full_ml) — always put the smaller model first
• After selecting: lmer(y ~ x + group + (1|site), data=d2, REML=TRUE)

---

5 Retrieval Challenge

NoteModel answers

Q1: No — REML likelihoods are not comparable when fixed effects differ, because the restricted likelihood depends on the column space of \(\mathbf{X}\). Use ML for fixed-effects comparison.

Q2: The colleague is right. lme4 automatically refits both models with ML when you call anova() on REML fits, and prints a warning. Read the output carefully — it tells you what it did.

---

6 Why-Questions

CautionQ: “REML is unbiased for variance components — so why not always use REML?”

REML integrates out the fixed effects from the likelihood. This means the REML likelihood is not a function of \(\boldsymbol{\beta}\) — you can’t use it to compare models where \(\boldsymbol{\beta}\) differs (different fixed effects mean different things were integrated out). ML keeps \(\boldsymbol{\beta}\) explicit, so nested fixed-effects models are genuinely comparable via LRT.

Analogy: if you use \(s^2 = \frac{\sum(x_i - \bar{x})^2}{n-1}\) (unbiased), you can’t directly compare it across datasets with different means without accounting for what was subtracted. ML is the “raw” version; REML applies a correction that depends on \(\mathbf{X}\).

CautionQ: “The LRT for fixed effects in a mixed model — is the chi-squared approximation exact?”

No. The LRT statistic \(2(\ell_{full} - \ell_{reduced})\) is asymptotically \(\chi^2_{df}\) under \(H_0\), but:

1. The approximation requires large \(n\) and large \(m\) (number of groups). With few groups, it can be anti-conservative.
2. For testing a single fixed effect, the \(F\)-test with Kenward-Roger or Satterthwaite denominator df (from lmerTest) is generally preferred in small samples.

Use library(lmerTest) and refit — it replaces \(t\)-values with \(t\)-tests with approximate df.

---

7 Spaced Review

NoteReview answer

Complete pooling: fast, simple — but biased SEs when ICC > 0. No pooling: correct inference for these specific groups — but wastes df, no generalisation. Partial pooling (mixed): efficient, generalisable, shrinks small-group estimates — preferred when groups are a sample from a population.

NoteFurther reading

Brown, H. & Prescott, R. (2015). Applied Mixed Models in Medicine (3rd ed.). Wiley. Companion datasets and code: wiley.com/go/brown/applied_mixed. — Sections 5.3–5.4.