Answer

REML or ML: which should a likelihood ratio test of fixed effects use?

Inspired by a question on Cross Validated · November 5, 2026

multilevel models model selection hypothesis testing

The short answer

Use maximum likelihood (ML). REML scores each model on data with its own fixed effects projected out, so REML fits with different fixed effects are scored on different data and their likelihoods cannot be compared. Refit the models with ML for the test, then refit the chosen model with REML to report it. REML comparisons are fine only when the fixed effects are identical.

What REML does differently

Maximum likelihood (ML) estimates the fixed effects and the variance components together from the raw outcome. Its variance estimates are biased downward, for the same reason that dividing a sum of squares by n instead of n - 1 underestimates a variance: the fitted fixed effects use up some of the information.

Restricted (or residual) maximum likelihood, REML, fixes this. It first removes everything the fixed-effect design can explain and then maximizes the likelihood of what is left, a set of error contrasts that do not depend on the fixed effects at all. That is why REML is the default in nlme, lme4 and most other mixed-model software: its variance components are less biased, and the standard errors built on them are more trustworthy, especially with few groups.

The price is that the restricted likelihood belongs to a particular fixed-effect design. Add or remove a predictor and the error contrasts change, so the model is now being scored on a different transformed data set.

Why REML likelihoods of different fixed-effect models cannot be compared

A likelihood ratio test needs both models to assign probabilities to the same data. With REML and different fixed effects, that is not true, so twice the difference in log-likelihoods has no chi-square distribution and its p-value means nothing. It is not even consistently too big or too small. The REML log-likelihood contains a term for the determinant of the fixed-effect design (weighted by the fitted covariance), so it changes if you merely change the units of a predictor, even though the fitted model is identical.

Software differs in how much it protects you. In R, anova() for nlme fits warns that REML comparisons are not meaningful when the fixed effects differ, and anova() for lme4 fits refits the models by ML before testing. Generic functions such as lmtest::lrtest() simply take whatever logLik() returns, so they report a statistic and a p-value without any warning. The same caution applies to AIC and BIC: do not rank REML fits that differ in their fixed effects by information criteria either.

See it in R and Python

The simulation has 30 people measured 6 times each, with a random intercept per person. The question is whether a second predictor z belongs in the model alongside x. Both models are fitted with ML and with REML using the nlme package, which ships with R. Then x is rescaled to thousands of years, which changes nothing about the fit, and the tests are repeated. The Python tab writes out the ML and REML log-likelihoods of the random-intercept model with NumPy and SciPy, using the same REML form as nlme.

library(nlme)   # ships with R
set.seed(677636)
J <- 30; n <- 6                          # 30 people, 6 visits each
id <- factor(rep(1:J, each = n))
x <- rnorm(J * n, 50, 10)                # e.g. age in years
z <- rnorm(J * n)                        # a second predictor (true effect 1.5)
y <- 100 + 0.5 * x + 1.5 * z + rnorm(J, 0, 4)[id] + rnorm(J * n, 0, 5)
d <- data.frame(y, x, z, id)

lrt <- function(small, big) {            # likelihood ratio test, 1 df
  LR <- 2 * (as.numeric(logLik(big)) - as.numeric(logLik(small)))
  round(c(LR = LR, p = pchisq(LR, df = 1, lower.tail = FALSE)), 4)
}

# Same two fixed-effect structures, fitted by ML and by REML
ml0   <- lme(y ~ x,     random = ~ 1 | id, data = d, method = "ML")
ml1   <- lme(y ~ x + z, random = ~ 1 | id, data = d, method = "ML")
reml0 <- lme(y ~ x,     random = ~ 1 | id, data = d, method = "REML")
reml1 <- lme(y ~ x + z, random = ~ 1 | id, data = d, method = "REML")
lrt(ml0, ml1)      # valid test of z
lrt(reml0, reml1)  # not a valid test

# Same models with x measured in thousands of years: identical fits
ml0k   <- lme(y ~ I(x / 1000),     random = ~ 1 | id, data = d, method = "ML")
reml0k <- lme(y ~ I(x / 1000),     random = ~ 1 | id, data = d, method = "REML")
round(c(ML = logLik(ml0k) - logLik(ml0), REML = logLik(reml0k) - logLik(reml0)), 4)
lrt(ml0k, ml1)     # ML test unchanged
lrt(reml0k, reml1) # REML "test" changes with the units of x

# After choosing the model, report REML variance components
sds <- function(m) as.numeric(VarCorr(m)[, "StdDev"])
round(rbind(ML = sds(ml1), REML = sds(reml1)), 3)  # person SD, residual SD
round(summary(reml1)$tTable["z", ], 4)

Python

import numpy as np
from scipy import optimize, stats

rng = np.random.default_rng(677636)
J, n = 30, 6                                  # 30 people, 6 visits each
idx = np.repeat(np.arange(J), n)
x = rng.normal(50, 10, J * n)                 # e.g. age in years
z = rng.normal(size=J * n)                    # a second predictor (true effect 1.5)
y = 100 + 0.5 * x + 1.5 * z + rng.normal(0, 4, J)[idx] + rng.normal(0, 5, J * n)

def loglik(theta, X, reml):                   # random-intercept model, beta profiled out
    s2u, s2e = np.exp(theta)
    c = s2u / (s2e + n * s2u)                 # V^-1 v = (v - c * person sum of v) / s2e
    Vinv = lambda v: (v - c * v.reshape((J, n) + v.shape[1:]).sum(axis=1)[idx]) / s2e
    XtViX = X.T @ Vinv(X)
    beta = np.linalg.solve(XtViX, X.T @ Vinv(y))
    r = y - X @ beta
    N, p = X.shape
    logdetV = J * ((n - 1) * np.log(s2e) + np.log(s2e + n * s2u))
    ll = -0.5 * (logdetV + r @ Vinv(r))
    if reml:                                  # same REML form as nlme
        return ll - 0.5 * (np.linalg.slogdet(XtViX)[1] + (N - p) * np.log(2 * np.pi)), beta, XtViX
    return ll - 0.5 * N * np.log(2 * np.pi), beta, XtViX

def fit(X, reml):
    res = optimize.minimize(lambda t: -loglik(t, X, reml)[0], np.log([10.0, 10.0]),
                            method="Nelder-Mead", options={"xatol": 1e-8, "fatol": 1e-10, "maxiter": 5000})
    ll, beta, XtViX = loglik(res.x, X, reml)
    return {"ll": ll, "beta": beta, "cov": np.linalg.inv(XtViX), "sd": np.sqrt(np.exp(res.x))}

def lrt(small, big):                          # likelihood ratio test, 1 df
    LR = 2 * (big["ll"] - small["ll"])
    return round(LR, 4), round(stats.chi2.sf(max(LR, 0), 1), 4)

one = np.ones(J * n)
X0, X1, X0k = np.column_stack([one, x]), np.column_stack([one, x, z]), np.column_stack([one, x / 1000])
ml0, ml1, reml0, reml1 = fit(X0, False), fit(X1, False), fit(X0, True), fit(X1, True)
print("ML LRT:", lrt(ml0, ml1))               # valid test of z
print("REML LRT:", lrt(reml0, reml1))         # not a valid test

# Same models with x measured in thousands of years: identical fits
ml0k, reml0k = fit(X0k, False), fit(X0k, True)
print("logLik change, ML and REML:", round(ml0k["ll"] - ml0["ll"], 4) + 0.0, round(reml0k["ll"] - reml0["ll"], 4))
print("ML LRT, rescaled:", lrt(ml0k, ml1))    # unchanged
print("REML LRT, rescaled:", lrt(reml0k, reml1))  # changes with the units of x

# After choosing the model, report REML variance components
print("ML SDs:", np.round(ml1["sd"], 3), "REML SDs:", np.round(reml1["sd"], 3))
b, se = reml1["beta"][2], np.sqrt(reml1["cov"][2, 2])
print("z: b, SE, t, p:", round(b, 4), round(se, 4), round(b / se, 4),
      round(2 * stats.t.sf(abs(b / se), J * n - J - 2), 4))

R output from one seeded run (same numbers on a second run):
     LR       p 
23.5849  0.0000 
     LR       p 
23.1958  0.0000 
    ML   REML 
0.0000 6.9078 
     LR       p 
23.5849  0.0000 
    LR      p 
9.3803 0.0022 
      [,1]  [,2]
ML   3.747 4.697
REML 3.828 4.728
    Value Std.Error        DF   t-value   p-value 
   1.8667    0.3750  148.0000    4.9779    0.0000

The figures below come from one seeded run of the R code. Python uses a different random number generator, so its figures differ (its estimate for z happens to be smaller), but every pattern described here is the same, including the exact shift of 6.9078.

The ML test is the valid one. Comparing the ML fits gives a statistic of 23.58 on 1 degree of freedom, p < .001, so z improves the model.
The REML number looks plausible, which is the danger. The REML difference is 23.20, close to the ML value here, so nothing seems wrong. But it is not a test statistic.
Changing units moves the REML likelihood. Expressing x in thousands of years leaves the ML log-likelihood unchanged but raises the REML log-likelihood of the smaller model by 6.9078, exactly log(1000). The ML test still gives 23.58; the REML "test" drops to 9.38 purely because of the units. In the Python run the same shift makes it negative.
REML still wins for the final estimates. The ML fit puts the person SD at 3.747 and the residual SD at 4.697; REML gives slightly larger values, 3.828 and 4.728, because it accounts for the three estimated fixed effects.

The workflow to use

To compare fixed effects, refit with ML. In nlme, update(model, method = "ML"); in lme4, refit = TRUE (the default in anova()) or REML = FALSE. Keep the data, random effects and correlation structure the same, and make sure both models use exactly the same rows, since missing values in an added predictor can silently drop observations.
Run the likelihood ratio test or compare AIC on those ML fits. The chi-square approximation can be optimistic with few groups; a Wald F test with Kenward-Roger or Satterthwaite degrees of freedom on the REML fit is a reasonable alternative for testing a single term.
Refit the chosen model with REML and report its fixed effects, standard errors, confidence intervals and variance components. That is standard practice: the comparison and the final estimates answer different questions.
To compare random effects or correlation structures, REML is fine as long as the fixed effects are identical. When the test is whether a variance is zero, the usual chi-square p-value is too large because zero is at the edge of the allowed range; with one variance tested, halving it is a common correction.

If the fixed-effect comparison is part of a larger search over many candidate models, remember that the final p-values no longer have their nominal meaning; see why stepwise regression is a bad idea. For help choosing an analysis in the first place, try our test chooser, and the DASS blog has more on regression and model building.

How to report a likelihood ratio test in APA style (7th edition)

Say in the method section that models were compared with ML and that the final model was estimated with REML, then give the chi-square statistic with its degrees of freedom. Using the example above:

"Nested models were compared with likelihood ratio tests on maximum likelihood fits, and the final model was re-estimated by REML with the nlme package in R. Adding z improved the fit of the random-intercept model, χ²(1) = 23.58, p < .001. In the final model, z predicted the outcome, b = 1.87, SE = 0.375, t(148) = 4.98, p < .001; the between-person SD was 3.83 and the residual SD was 4.73."

When several models are compared, a table with each model's log-likelihood, number of parameters, AIC and the χ² test against the previous model is easier to read than a list of statistics in the text.

Related tools and guides

Working with your own data?

General answers only go so far. Send us your situation and we'll reply by email within two business days.

Ask your question

Written by AskStats with AI assistance. This is general information, not advice for your specific data or study. When the results matter, check your approach with a qualified statistician.