Comparing REML and ML in GWAS Null Models

Why the Choice of Likelihood Matters

A mixed-model GWAS null model estimates variance components before testing each genetic marker. Those estimates determine the covariance structure used to account for relatedness and other modeled sources of variation. Maximum likelihood (ML) and restricted maximum likelihood (REML) use the same Gaussian model but treat the fixed effects differently when estimating that covariance structure. As a result, the estimated variance components can differ, particularly when the number of fixed effects is large relative to the sample size.

The distinction is not that one method accounts for fixed effects and the other ignores them. Both estimate fixed effects. The distinction is how their likelihoods use the information remaining after those effects are fitted. That distinction matters in a GWAS because the null model contains nuisance fixed effects—commonly an intercept and prespecified covariates—but excludes the marker currently being tested.

Model and Assumptions

The model assumes that the mean is correctly represented by \(X\beta\), that the covariance has the stated form, and that the Gaussian distribution is an adequate working model. The variance components are nonnegative, and \(V\) must be positive definite for the usual likelihood calculations. In practice, the observations, phenotype, covariates, and relatedness matrix must also be aligned consistently. A mismatch in sample order or an inappropriate covariance matrix cannot be repaired by choosing REML instead of ML.

For a GWAS null model, \(X\) includes the intercept and the chosen non-genetic or ancestry-related covariates, but not the marker being tested. The covariance parameters \(\theta=(\sigma_g^2,\sigma_e^2)\) are estimated from this model. Once fitted, the null model supplies a covariance estimate for marker association testing. The precise association procedure is not the subject here; the key point is that null-model fitting estimates covariance parameters conditional on a particular fixed-effect design.

What ML and REML Maximize

For a proposed covariance matrix \(V\), ML estimates \(\beta\) by generalized least squares:

Define

Then \(y'Py\) is the generalized residual quadratic form after fitting the fixed effects. The ML profile likelihood for the covariance parameters uses this residual fit and the determinant of \(V\):

REML instead bases covariance estimation on contrasts that remove the fixed effects. Let \(L\) be an \(n\times(n-p)\) matrix of full column rank satisfying \(L'X=0\). Then \(L'y\) has mean zero and covariance \(L'VL\). REML maximizes the likelihood of these error contrasts. Equivalently, for a fixed design \(X\), its log likelihood can be written up to constants independent of \(\theta\) as

The determinant term involving \(X\) is part of the adjustment for estimating fixed effects. REML does not simply replace \(n\) by \(n-p\) in every mixed model; the displayed expression shows that its correction also depends on \(V\) and the design matrix. The likelihood is defined relative to a fixed set of contrasts, so its value is not directly comparable across different fixed-effect designs.

A Short Derivation: Why the Simple Variance Estimates Differ

The degrees-of-freedom difference can be seen exactly in the special case \(V=\sigma^2I\), with \(\sigma^2\) the only covariance parameter. For any fixed design \(X\), let

The ML profile likelihood is proportional to \((\sigma^2)^{-n/2}\exp\{-\mathrm{SSE}/(2\sigma^2)\}\). Taking its log and differentiating with respect to \(\sigma^2\) gives

For REML, the contrast vector has dimension \(n-p\), and its covariance is \(\sigma^2L'L\). Its log likelihood, as a function of \(\sigma^2\), contains \(-\frac{n-p}{2}\log(\sigma^2)-\frac{\mathrm{SSE}}{2\sigma^2}\). Differentiating yields

Under this independent, equal-variance model, \(\mathrm{SSE}/\sigma^2\) has a chi-squared distribution with \(n-p\) degrees of freedom, so this REML estimate is unbiased for \(\sigma^2\). The ML estimate is smaller on average because it divides by \(n\) even though \(p\) fixed effects were fitted. This exact result explains the basic adjustment, but it does not imply that every REML variance-component estimate is unbiased in a general mixed model.

Worked GWAS-Style Null-Model Calculations

Implications for Fitting a GWAS Null Model

When the null fixed-effect design is specified in advance and remains the same across a scan, REML is commonly used to estimate mixed-model covariance parameters. It accounts for the dimensions consumed by the fixed effects, which can matter in finite samples. ML is also a valid likelihood method under the model, but its variance-component estimates can be downward-biased in settings where many fixed effects have been estimated. In a general mixed model, the size and direction of the difference depend on the covariance structure and data; the simple denominator correction is not a universal formula.

The choice should follow the inferential task. For estimating covariance parameters under one fixed-effect design, REML is often a sensible default. If the scientific comparison changes the fixed-effect design, REML likelihood values from those designs are not directly comparable; ML is generally used for likelihood-based fixed-effect comparisons. This is a statement about likelihood comparison, not a claim that an ML comparison automatically solves every issue of model selection or test calibration.

A practical null-model fit should specify the fixed effects before estimating variance components. Adding or removing covariates changes the residual space and may change the variance estimates. When comparing fits for a GWAS scan, confirm that the phenotype, sample, fixed-effect design, and covariance structure are consistent. If marker-specific missingness changes the analyzed individuals, the null model and covariance matrix may no longer correspond to the same observations; simply reusing a fit without checking this alignment can be inappropriate.

Diagnostics, Alternatives, and Limitations

REML does not guarantee a well-behaved fit. Variance components can be estimated at a boundary, such as \(\widehat{\sigma}_g^2=0\), or the likelihood can be nearly flat in a direction that makes components difficult to distinguish. A nearly singular \(K\), a poorly conditioned fixed-effect design, or strong overlap between modeled covariance components can make optimization unstable. Inspect convergence information, parameter boundaries, and whether the fitted covariance is positive definite. Check that the covariate design has full column rank and that the intended observations are in the same order in \(y\), \(X\), and \(K\).

The model assumptions also deserve attention. Residual patterns inconsistent with the Gaussian mean or covariance model, influential observations, or a misspecified \(K\) can undermine both ML and REML. These methods estimate parameters within the model supplied; neither corrects confounding or makes an unsuitable covariance structure appropriate. If the phenotype is binary, the Gaussian model described here may be inappropriate, and methods for non-Gaussian mixed models require their own likelihood and estimation considerations.

A frequent mistake is to say, “REML is always unbiased, so it is always the right choice.” The correction is narrower: REML removes fixed-effect directions from the likelihood, and in the simple independent Gaussian model it gives the unbiased residual variance estimate. General mixed-model REML estimates are not guaranteed to be unbiased, and the choice still depends on whether the goal is covariance estimation or comparison of fixed-effect structures.