Why Compare Two Mixed Models?
A mixed-model association test asks whether a genetic variant explains additional variation in a phenotype after accounting for covariates and relatedness. The likelihood-ratio test (LRT) answers this by comparing how well two nested models fit: a null model without the variant effect and a full model that includes it. It is a direct comparison of fitted likelihoods, rather than a test based on one model coefficient alone.
This comparison is useful when the variant effect is the scientific target and the fitted model is well specified. Its main practical cost in a genome-wide scan is that the full model may need to be fitted separately for each variant. Unlike a score test, which can use a fitted null model to screen variants efficiently, an LRT generally requires fitting the alternative model as well.
Models and the Likelihood-Ratio Statistic
For a continuous phenotype measured on \(n\) individuals, write \(y\) for the phenotype vector, \(X\) for the matrix of nuisance fixed-effect covariates, \(g\) for a vector of variant dosages, and \(K\) for a specified genomic relationship matrix. A common linear mixed model is
Here \(\beta\) contains the covariate effects, \(\gamma\) is the variant effect per dosage unit, and \(\theta\) denotes the covariance parameters, such as \(\tau^2\) and \(\sigma^2\). The random genetic effect represented through \(K\) accounts for covariance among individuals; it is not the variant effect being tested. The null model sets \(\gamma=0\), while the full model estimates \(\gamma\) along with the nuisance parameters.
For a single additive dosage effect, the null hypothesis is \(H_0:\gamma=0\) and the alternative is \(H_a:\gamma\ne 0\). The full model adds one parameter, so the usual reference distribution is \(\chi^2_1\). If the variant is represented by two independent effect terms, such as a two-degree-of-freedom genotype coding, the joint null sets both terms to zero and the usual reference distribution is \(\chi^2_2\).
The ordinary Gaussian log likelihood at parameters \((\beta,\gamma,\theta)\) is
The maximized likelihood includes the fitted fixed effects and covariance parameters. For a conventional LRT, maximize this likelihood separately under the null and full models, using the same individuals, phenotype, covariate specification, and covariance-model family. If the full model contains the null as a special case and both optimizations find the global maxima, then \(\ell_F\geq\ell_0\) and \(T\geq0\). A substantially negative computed value signals an optimization or model-comparison problem.
Calibration and Its Justification
The chi-square approximation is an asymptotic result, not a rule that applies automatically to every pair of mixed models. It relies on a correctly specified likelihood, nested models, identifiable parameters, and regular estimation behavior near the null. In particular, the tested effect should be an ordinary, freely varying fixed-effect parameter under the alternative, and nuisance parameters should not be in problematic boundary or non-identifiable configurations. With a single variant effect and adequately large sample size, the usual \(\chi^2_1\) approximation is common, but small samples, rare variants, or irregular variance-component fits can make it inaccurate.
The basic justification follows from a local quadratic approximation to the log likelihood. Let \(\eta\) denote the \(q\) added parameters, and let \(\lambda\) collect nuisance parameters. Near the null, the profile log likelihood for \(\eta\) is approximately a quadratic function. After scaling by its information matrix, this can be written locally as a constant plus a linear term minus one-half a squared distance. Maximizing over the \(q\) unrestricted coordinates gains approximately one-half the squared length of a \(q\)-dimensional standard normal vector, compared with fixing those coordinates at zero. Thus twice the gain is approximately the sum of \(q\) squared independent standard normal variables, which has a \(\chi^2_q\) distribution. This is the likelihood-ratio form of the standard asymptotic result for nested regular models.
This argument also explains the degrees of freedom: they count the independent restrictions released in the full model, not the number of observations or the number of covariance parameters. If a variance component is itself tested at zero, it lies on the boundary of its nonnegative parameter space; the regular argument above does not apply, and the reference distribution may be a mixture rather than an ordinary chi-square. A variant fixed-effect test is a different hypothesis, but boundary fits in nuisance variance components can still make routine calibration less dependable.
Worked Example: One Additive Variant Effect
Suppose a hypothetical study analyzes a continuous biomarker measured in 1,200 individuals. The fitted null and full models use the same individuals, covariates, and relationship matrix. The Gaussian covariance model is judged adequate, the dosage effect is an unconstrained fixed effect, and both models converge to their maximum-likelihood fits. The maximized log likelihoods are \(\ell_0=-128.40\) and \(\ell_F=-123.795\). The estimated dosage coefficient is \(0.18\) biomarker units per additional effect allele.
The full model adds one parameter, so the test has one degree of freedom. The statistic is
The upper-tail probability under \(\chi^2_1\) is approximately \(0.0024\). Under the stated model and approximation, this is evidence against \(H_0:\gamma=0\). The coefficient estimates an increase of \(0.18\) biomarker units per additional effect allele; the LRT p-value is evidence about the null, not the probability that the null is true. The assumptions are supported here by the stated common sample and model structure, convergence, and regular fixed-effect test, though model adequacy and the chi-square approximation still require scrutiny in a real analysis.
Worked Comparisons and Model Choices
Worked Example: A Two-Parameter Variant Test
In a hypothetical study of 800 individuals, a variant is represented by two genotype-effect terms rather than a single additive dosage term. The null model has neither term; the full model estimates both. Assume the same observations and nuisance structure are used, the likelihood is correctly specified, and both fits converge to their maximum-likelihood values. Suppose \(\ell_0=-90.000\) and \(\ell_F=-86.0925\).
The statistic is
The full model adds two parameters, so compare \(T\) with \(\chi^2_2\). Its upper-tail probability is approximately \(0.0201\). This joint test is evidence that at least one of the two modeled genotype effects is nonzero; it does not identify which term is responsible. That would require an additional analysis. The degrees of freedom are two because the null imposes two restrictions, not because there are two genotype groups or because the sample has 800 individuals.
Worked Example: Why REML Likelihoods Cannot Be Compared This Way
Consider a hypothetical single-variant test in which the full model adds genotype to the fixed-effect design. Restricted maximum likelihood (REML) removes fixed effects using a construction that depends on that design. Suppose software reports restricted log likelihoods of \(-53.0\) for the null design and \(-50.2\) for the full design.
A tempting calculation is \(2\{-50.2-(-53.0)\}=5.6\). But this does not provide the ordinary LRT: the two restricted likelihoods are based on different fixed-effect designs, so their difference is not a valid likelihood ratio for testing the added fixed effect. The fact that the number is positive does not repair the comparison.
Instead, fit both models by maximum likelihood. Suppose the corresponding maximized log likelihoods are \(-61.0\) and \(-58.9\). Then
For one added parameter, the \(\chi^2_1\) upper-tail probability is approximately \(0.0404\). This is the conventional LRT calculation, provided the regularity and model assumptions hold. REML remains useful for estimating covariance parameters, but comparing REML fits with different fixed-effect structures is not the way to test the added variant effect.
Computational and Inferential Trade-Offs
The three familiar tests use different parts of the model fit. A Wald test assesses the fitted variant coefficient relative to its standard error. An LRT compares the maximized fit with and without the variant. A score test evaluates evidence for a departure from the null using the null fit. For a regular, well-behaved large-sample problem, these tests often give similar results, but they are not numerically identical in every sample.
| Test | What it uses | Typical computational feature | Main caution |
|---|---|---|---|
| Wald | Full-model estimate and standard error | Requires fitting the alternative model | Can be sensitive to weak identification or an inaccurate local quadratic approximation |
| Likelihood ratio | Maximized null and full likelihoods | Usually requires fitting the alternative for each variant | Needs comparable likelihoods and appropriate calibration |
| Score | Evidence evaluated under the null model | Can be efficient for screening many variants from a fitted null | Tests locally at the null and can differ from full-model evidence in finite samples |
For a scan over millions of variants, repeatedly optimizing the covariance parameters in each full model may be prohibitively expensive. Implementations can reduce the burden through reuse of null-model quantities or specialized algorithms, but the conceptual LRT still depends on the alternative-model likelihood. A score test is often preferred for rapid genome-wide screening; an LRT can be useful for follow-up or when direct comparison of fitted models is important. Neither method protects against an incorrect covariance model, genotype coding error, or unmodeled confounding.
In practice, inspect convergence messages, fitted variance components, and whether the full likelihood is at least as large as the null likelihood. Check that the same individuals and covariates are used in both fits, and that the variant is not nearly constant or otherwise uninformative. A coefficient estimate and its uncertainty provide context for the test statistic, while residual and model diagnostics help assess whether the assumed phenotype distribution and covariance structure are plausible. A tiny p-value cannot compensate for a poorly specified model.
Common Mistakes and Limitations
- Using the wrong sign or omitting the factor of two. The statistic is \(2(\ell_F-\ell_0)\), with full minus null. Reversing the subtraction can produce a negative value despite correct nested fits.
- Comparing REML likelihoods after changing fixed effects. REML values depend on the fixed-effect design. Use maximum likelihood for the ordinary fixed-effect LRT comparison.
- Choosing degrees of freedom from the sample size or number of genotypes. Count the independent parameters added by the full model.
- Treating the chi-square approximation as exact. It is an asymptotic reference under regularity conditions; rare variants, small samples, failed fits, and boundary complications can undermine calibration.
- Interpreting a p-value as an effect size. The p-value measures compatibility of the data with the null under the model. Report the estimated effect and its units separately.
The LRT is therefore not automatically more accurate than a Wald or score test. Its advantage is a direct comparison of nested likelihoods; its disadvantages include alternative-model fitting costs and dependence on likelihood comparability and calibration. The test is most informative when the model assumptions are defensible, the optimization is reliable, and the reported result is paired with an interpretable effect estimate.
Check Your Understanding
Try each question, then compare your reasoning with the solutions below.
- A full model adds one dosage coefficient to a null mixed model. The maximized log likelihoods are \(-45.0\) and \(-42.8\), respectively. Calculate the LRT statistic and degrees of freedom.
- Why is a difference of REML log likelihoods generally unsuitable when the full model adds a fixed variant effect?
- A full model adds three independent variant-effect parameters. What chi-square reference distribution is used under regular conditions?
- In a genome-wide scan, why might a score test be used for screening and an LRT for selected follow-up variants?
Solutions. (1) \(T=2\{-42.8-(-45.0)\}=2(2.2)=4.4\), with \(df=1\). (2) The restricted likelihood depends on the fixed-effect design, so changing that design makes the REML likelihoods non-comparable for this ordinary fixed-effect LRT; fit the models by maximum likelihood instead. (3) Use \(\chi^2_3\), because three independent restrictions are released, subject to regularity conditions. (4) A score test can reuse the fitted null model and is often cheaper across many variants; an LRT directly compares null and full fits but usually requires alternative-model fitting, which may be practical for a smaller follow-up set.