ML and REML Comparisons for Random-Effect Structures

Why the likelihood calculation matters

A richer random-effects structure can describe more sources of variation, but its likelihood value is comparable with that of a simpler structure only when both values refer to the same data and the same probability model conventions. The distinction between maximum likelihood (ML) and restricted maximum likelihood (REML) is especially important: ML is based on the response likelihood, whereas exact REML in a Gaussian linear mixed model is based on likelihood information in contrasts that eliminate the fixed effects.

This tutorial concerns comparability, not the reference distribution for a likelihood-ratio test. A positive difference in optimized likelihoods can be calculated only after checking that the fits are comparable; its inferential calibration is a separate question. In particular, variance-component tests may involve parameter-space boundaries, as in the earlier material on variance components and boundary tests.

ML for random-effects structures

Write a GLMM in terms of responses \(Y_i\), fixed-effects design matrix \(X\), fixed-effect vector \(\beta\), and random effects \(b\). Conditional on \(b\), the responses follow the specified exponential-family response model, with linear predictor \(X\beta+Zb\) (including any stated offset). A usual random-effects assumption is \(b\sim N(0,G_\psi)\), independent of the response sampling given \(b\). The marginal likelihood integrates out the random effects:

ML estimates \(\beta\) and the covariance parameters \(\psi\) by maximizing this marginal log likelihood. A likelihood value based only on the conditional response density evaluated at an estimated or predicted \(b\) is not this marginal likelihood. The integration may be evaluated numerically or approximated, for example by a Laplace approximation; the method used must be reported and applied consistently to both fits.

For a nested comparison, let model 0 be the smaller random-effects structure and model 1 the larger. If the fixed-effects specification and all other likelihood conventions are the same, the parameter space of model 0 is contained in that of model 1. Therefore the optimized ML log likelihoods satisfy:

This inequality follows because maximizing over a parameter space that includes every value allowed by the smaller model cannot produce a lower maximum. A materially negative difference is a warning to check convergence, model nesting, data alignment, and whether the two reported likelihoods use the same convention. Tiny negative differences can arise from numerical optimization error.

What REML compares

Exact REML is naturally defined for the Gaussian linear mixed model \(y\sim N(X\beta,V_\theta)\), where \(V_\theta\) is the marginal covariance matrix determined by the residual and random-effects parameters \(\theta\). Let \(X\) have full column rank \(p\). Choose an \(n\times(n-p)\) matrix \(K\) whose columns are orthonormal and span the orthogonal complement of the columns of \(X\). Thus \(K'K=I\) and \(K'X=0\). The contrast vector \(K'y\) has a normal distribution with mean zero and covariance \(K'V_\theta K\). Its log likelihood is the exact restricted log likelihood:

The contrast likelihood eliminates the unknown fixed effects because \(K'X\beta=0\). It depends on the fixed-effects column space, since that space determines which directions are removed, rather than on a particular scaling of its basis columns. Consequently, exact REML likelihoods for random-effects structures are comparable when the data and the fixed-effects column space are the same. Changing the fixed-effects column space changes the contrasts and hence changes the likelihood being evaluated.

For a full-column-rank \(X\), the same contrast likelihood can be written using \(P_\theta=V_\theta^{-1}-V_\theta^{-1}X(X'V_\theta^{-1}X)^{-1}X'V_\theta^{-1}\):

The \(-\log|X'X|\) term inside the brackets is important. If the alternative determinant formula is written without it, that expression is not invariant to rescaling the columns of \(X\). The inverse and determinants involving \(X\) also require full column rank. If a design is rank deficient, use a full-rank basis for its column space or work directly with orthonormal contrasts. When the same design matrix parameterization is used in every fit, the \(\log|X'X|\) term is constant across those fits and cancels in their difference; it must not be dropped when claiming the formula is invariant to a change of parameterization.

A practical consequence is that REML likelihoods are not the appropriate basis for comparing different fixed-effects column spaces. Such comparisons alter the contrast construction. For Gaussian mixed models, use ML when fixed-effects specifications differ; for random-effects comparisons with the same fixed-effects space, REML can be used. The introductory material on ML and REML comparisons for fixed effects covers the corresponding fixed-effects comparison issue.

GLMMs, diagnostics, and alternatives

For a non-Gaussian GLMM, the exact Gaussian contrast construction above generally does not provide a REML likelihood: linear contrasts of a non-Gaussian response do not generally have the required normal likelihood. Some software offers a restricted or adjusted likelihood approximation for GLMMs, but the construction can depend on the approximation and implementation. A reported GLMM “REML” value should therefore not automatically be treated as the exact Gaussian restricted likelihood or compared across packages and methods. For GLMM random-effects comparisons, marginal ML is the usual likelihood basis; if an approximate restricted criterion is used, its definition and comparability conditions must be established for that specific method.

Before comparing fits, check the following:

• Same response data: confirm identical observations, response coding, missing-data handling, offsets, and prior or sampling weights.
• Same model conventions: use the same response distribution, link, integration method, approximation settings, and likelihood normalization.
• Correct nesting: verify that the smaller random-effects model is obtained by constraints within the larger model, not by changing the fixed-effects part or data.
• Successful fitting: inspect convergence messages, gradients, optimizer sensitivity, and whether estimated covariance matrices are singular or near singular.
• Numerical accuracy: for approximate marginal likelihoods, check whether quadrature or integration settings materially change the reported values.

AIC or another information criterion can be an alternative for comparing non-nested candidate structures when their likelihoods refer to the same response data and likelihood convention. It does not repair incompatible likelihoods, and a REML-based criterion should not be used to compare different fixed-effects spaces. Predictive validation answers a different question—how well models predict new data—and may be preferable when prediction is the goal. Neither an information criterion nor predictive performance replaces checking model assumptions and the scientific meaning of the random-effects structure.

The most important limitation is that likelihood comparability is not the same as inferential validity. A valid difference in ML or REML values does not alone determine the test’s null distribution, establish that the selected structure is scientifically adequate, or remove sensitivity to sparse data and weakly identified covariance parameters. Those issues require separate assessment.