Wilks' Theorem for Regular Fixed-Effect Tests

Why a Likelihood-Ratio Statistic Has a Chi-Square Limit

A likelihood-ratio test compares how well two nested models describe the same data. Its statistic measures the improvement in maximized log likelihood when the less restricted model is used. For a regular fixed-effect test, Wilks’ theorem explains why this improvement, after doubling it, has an asymptotic chi-square distribution when the null hypothesis is true. The result is not a property of every likelihood comparison: it depends on the model being identifiable and sufficiently far from irregularities such as parameter-space boundaries.

In a generalized linear mixed model, the likelihood used for this result is ordinarily the marginal likelihood: the random effects are integrated out, rather than treated as freely estimated parameters. The theorem describes the large-sample behavior of a test based on that likelihood when the number of independent clusters grows and the model satisfies regularity conditions.

Definition: Let $\ell_F$ be the maximized marginal log likelihood under a full GLMM and $\ell_R$ the maximized marginal log likelihood under a nested reduced GLMM. The likelihood-ratio statistic is $$ D = 2(\ell_F-\ell_R). $$ For a regular test of $q$ independent fixed-effect restrictions, Wilks’ theorem states that, under the null hypothesis, $D$ converges in distribution to a chi-square random variable with $q$ degrees of freedom as the number of independent clusters increases.

For example, a test may compare a model allowing a treatment coefficient to vary freely with one that fixes that coefficient at zero. More generally, the null may impose several independent restrictions on a vector of fixed effects. The models must be nested: every parameter value permitted by the reduced model must also be permitted by the full model. Both fits must use the same observed data and the same marginal-likelihood basis for comparison.

Conditions for the Result

The phrase “under regularity conditions” carries important content. A useful GLMM setting has independent clusters indexed by $i=1,\ldots,m$. Within cluster $i$, responses are conditionally independent given its random effects, and their conditional distributions follow the specified GLM family and link. Random effects have a specified distribution, commonly multivariate normal, and are integrated out to obtain a cluster marginal density $f_i(y_i;\theta)$. The marginal log likelihood is $$ \ell_m(\theta)=\sum_{i=1}^{m}\log f_i(y_i;\theta), $$ where $\theta$ includes the fixed effects and nuisance parameters, such as random-effect covariance parameters and any dispersion parameter.

For the usual Wilks result, the null parameter value must be an interior point of a correctly specified, identifiable model. The likelihood must be sufficiently smooth near that value; its information matrix must be finite and nonsingular in the relevant directions; and the information accumulated across clusters must grow in a stable, nondegenerate way. The tested restrictions must be smooth and locally independent. The random-effects distribution and covariance parameters also need to be identifiable under the null. In particular, the standard result does not automatically apply when a nuisance parameter is on the boundary of its allowable space.

These conditions are asymptotic, not a claim that a finite data set has exactly a chi-square test statistic. The standard setting is an increasing number of independent clusters, with no single cluster dominating the information. A very small number of clusters, severe sparsity, separation, or weak identification can make the approximation unreliable even if a fitting algorithm returns estimates.

Why Maximizing Over Nuisance Parameters Still Gives Chi-Square

The intuition is that, near the true parameter, a smooth log likelihood is approximately quadratic. After nuisance parameters are allowed to adjust, the part of the likelihood improvement attributable to the tested restrictions is a squared, standardized score. A sum of $q$ squared independent standard-normal quantities has a chi-square distribution with $q$ degrees of freedom.

Here is the local argument. Write the parameter in local coordinates as $(\psi,\lambda)$, where $\psi$ is the $q$-dimensional component constrained by the null and $\lambda$ denotes nuisance parameters. At the null, let the score and information be partitioned into corresponding blocks. Removing the contribution explained by nuisance-parameter adjustment gives the efficient score and its information:

U_{\psi\cdot\lambda} = U_\psi-I_{\psi\lambda}I_{\lambda\lambda}^{-1}U_\lambda, \qquad I_{\psi\cdot\lambda} = I_{\psi\psi}-I_{\psi\lambda}I_{\lambda\lambda}^{-1}I_{\lambda\psi}.

Under the null and the regularity conditions, the efficient score, standardized by its information, has an asymptotic multivariate standard-normal distribution:

I_{\psi\cdot\lambda}^{-1/2}U_{\psi\cdot\lambda} \ \xrightarrow{d}\ N_q(0,I_q).

The local quadratic approximation to the profile log likelihood says that the improvement from freeing $\psi$, while re-optimizing $\lambda$, is approximately one half of the squared length of this standardized score. Since the likelihood-ratio statistic doubles that improvement,

D \ \approx\ U_{\psi\cdot\lambda}^{\mathsf T} I_{\psi\cdot\lambda}^{-1} U_{\psi\cdot\lambda} = \left\|I_{\psi\cdot\lambda}^{-1/2}U_{\psi\cdot\lambda}\right\|^2 \ \xrightarrow{d}\ \chi^2_q.

The last equality in distribution follows because the squared length of a $q$-dimensional standard-normal vector is the sum of $q$ squared standard-normal variables. This also explains why nuisance parameters do not simply add degrees of freedom to the reference distribution: they are estimated in both models, and their local effect is accounted for through the efficient score and information. The result concerns the number of independent restrictions being tested, not the total number of parameters fitted.

Worked Example: A One-Restriction Treatment Test

Suppose hypothetical data come from 80 independent clinics. A binary-outcome GLMM includes a clinic random intercept and a treatment fixed effect $\beta_T$. The null hypothesis is $H_0:\beta_T=0$; the full model estimates $\beta_T$, while the reduced model fixes it at zero. Assume both models converged, the random-intercept variance estimate is positive and interior, the treatment effect is identifiable, and the marginal likelihood is well behaved. These facts make the regular chi-square approximation plausible for this illustration.

Suppose the maximized marginal log likelihoods are $\ell_F=-128.40$ and $\ell_R=-131.00$. The statistic is $$ D=2(\ell_F-\ell_R) =2[-128.40-(-131.00)] =2(2.60)=5.20. $$

There is one independent restriction, so the reference distribution is $\chi^2_1$. The upper-tail p-value is $P(\chi^2_1\geq 5.20)\approx 0.0226$. At significance level $\alpha=0.05$, this is evidence against the null of no treatment fixed effect, conditional on the model and its assumptions. It does not establish that treatment causes an outcome change; interpretation also depends on study design and the specified model.

Worked Example: A Joint Test of Two Restrictions

Consider a hypothetical GLMM for a count outcome, with a random intercept for each of 60 independent sites. A full model includes two prespecified fixed-effect coefficients, $\beta_1$ and $\beta_2$, while a reduced model fixes both at zero. Assume these are two identifiable, independent restrictions, and the nuisance estimates are interior. Suppose the full and reduced maximized marginal log likelihoods are $-205.70$ and $-209.80$, respectively.

The likelihood-ratio statistic is $$ D=2[-205.70-(-209.80)] =2(4.10)=8.20. $$

For two restrictions, the reference distribution is $\chi^2_2$. Its upper-tail probability at 8.20 is $e^{-8.20/2}=e^{-4.10}\approx 0.0166$. Thus, at the 0.05 level, the joint null hypothesis is rejected. This is evidence that the two coefficients are not both zero under the specified model; it does not, by itself, identify which coefficient contributes most to the evidence.

A Local Quadratic Illustration

A one-dimensional version makes the approximation concrete. Suppose a marginal log likelihood near the null value $\psi=0$ is locally $\ell(\psi)-\ell(0)\approx U\psi-\tfrac12 I\psi^2$, where $U$ is the score at zero and $I>0$ is the information. Maximizing this quadratic gives $\hat{\psi}=U/I$. Substitution yields a maximized log-likelihood gain of $U^2/(2I)$, so the likelihood-ratio statistic is $U^2/I$. Under the null, $U/\sqrt{I}$ is approximately standard normal; consequently $U^2/I$ is approximately $\chi^2_1$.

Worked Example: Turning a Local Likelihood Gain into a Test

For a hypothetical fitted GLMM, suppose the local quadratic approximation for a single tested coefficient has score $U=2.40$ and information $I=4.00$. The approximate maximizing value is $$ \hat{\psi}=\frac{U}{I}=\frac{2.40}{4.00}=0.600. $$

The approximate log-likelihood gain is $U^2/(2I)=2.40^2/(2\times4.00)=5.76/8.00=0.720$. Doubling it gives $D\approx1.44$, equivalently $U^2/I=5.76/4.00=1.44$. Using a $\chi^2_1$ reference distribution gives an upper-tail p-value of approximately $0.230$. This local illustration does not replace fitting the two models; it shows how the quadratic likelihood approximation produces the squared-normal form underlying Wilks’ theorem.

Using the Approximation and Recognizing Problems

In practice, the theorem justifies comparing the observed statistic with its chi-square reference distribution only after checking that the comparison is meaningful. The full and reduced models should represent the same response data, clusters, random-effects structure, and marginal-likelihood framework, differing only in the fixed-effect restrictions under test. The optimization should reach stable solutions, with no unexplained warnings, extreme estimates, or large discrepancies across reasonable starting values. A numerical maximization failure can distort either log likelihood and therefore the difference $D$.

The principal alternatives answer related questions but have different finite-sample behavior. A Wald test evaluates an estimate relative to its estimated standard error; a score test evaluates the likelihood slope at the null fit. Both can be convenient, but neither is the same finite-sample calculation as comparing the two maximized likelihoods. When the cluster count is small or the regularity approximation is suspect, a parametric bootstrap can instead simulate data under the fitted null model, refit both models for each simulated data set, and compare the observed statistic with the simulated null distribution. That procedure requires its own modeling and computational checks; it does not make a misspecified model valid.

Common pitfall: “The statistic is a likelihood ratio, so it must have a chi-square distribution.” This is wrong: the chi-square law is an asymptotic result under regularity conditions and a true null. In particular, an estimated variance component on the boundary, nonidentifiability, very few independent clusters, or unstable fits can invalidate the ordinary Wilks approximation. Check the model and fitting behavior before relying on the reference distribution.

Wilks’ theorem is valuable because it provides a general calibration for likelihood improvement without requiring a separate derivation for every fixed-effect model. Its limitation is equally important: the theorem guarantees an asymptotic reference law, not exact small-sample calibration. A small p-value is evidence against the specified null under the assumed likelihood and regularity conditions; it is not the probability that the null is true.

Check Your Understanding

Answer each question using the conditions and calculations in this tutorial.

A full model has maximized log likelihood $-74.2$, and its reduced nested model has maximized log likelihood $-76.0$. Calculate $D$.
For the statistic in question 1, suppose one independent restriction is tested. Is the 0.05-level test significant? Give an approximate p-value.
Why are the random-effect covariance parameters not counted as tested degrees of freedom when they are nuisance parameters estimated in both models?
Name two conditions that can make the ordinary Wilks chi-square approximation unreliable in a GLMM.
Under the null, what is the interpretation of a p-value from a Wilks-based test?

Solutions. (1) $D=2[-74.2-(-76.0)]=2(1.8)=3.60$. (2) With one restriction, compare with $\chi^2_1$: the upper-tail p-value is approximately $0.058$, so the test is not significant at 0.05. (3) Wilks’ result concerns the restrictions imposed by the null; nuisance parameters are re-optimized in both models, and their local adjustment is incorporated into the efficient information. (4) Examples include a covariance parameter on the boundary, nonidentifiability, very few independent clusters, or unstable likelihood maximization. (5) It is the probability, assuming the null and model conditions hold, of obtaining a test statistic at least as large as the observed one; it is not the probability that the null is true.