Tutorials › Generalized Linear Mixed Models (GLMM) › Testing Fixed-Effect Interactions in Clustered Data

Testing Fixed and Random Effects · Tutorial 365 of 698

Testing Fixed-Effect Interactions in Clustered Data

Learn to define the interaction null hypothesis and choose a joint test that matches the interaction’s degrees of freedom and model structure.

Advanced 11 min read

What You'll Learn

  • Write interaction terms in a GLMM for categorical and continuous predictors.
  • State the joint null hypothesis for all coefficients in an interaction block.
  • Distinguish a joint interaction test from tests of individual coefficients.
  • Calculate a joint Wald statistic using the covariance among interaction estimates.
  • Compare nested models with a likelihood-ratio test for an interaction.
  • Interpret interaction effects on the model’s link scale and check key assumptions.

Why Test an Interaction as a Block?

An interaction asks whether the association between one predictor and the outcome changes across values or levels of another predictor. In clustered data, that question is answered within a generalized linear mixed model (GLMM), which represents dependence through random effects while modelling the outcome conditionally on those effects. The central testing issue is that an interaction can involve several coefficients. Testing one coefficient at a time does not, in general, test whether the interaction as a whole is present.

For example, if a three-level group factor interacts with time, the interaction typically has two coefficients when one group is the reference. The joint null hypothesis sets both coefficients to zero. A test of just one coefficient asks whether one particular contrast is zero; it leaves the other contrast unrestricted. These are different hypotheses and can lead to different conclusions.

Definition: A fixed-effect interaction test evaluates whether a specified block of interaction coefficients is jointly zero, conditional on the other fixed effects and the random-effects structure in the model. If the block is \(\boldsymbol{\beta}_{I}\), the joint null is \(H_0:\boldsymbol{\beta}_{I}=\boldsymbol{0}\). Its degrees of freedom are the number of independent restrictions being tested, under the regularity conditions for the chosen test.

Formulating the Interaction in a GLMM

Let \(Y_{ij}\) denote an outcome for observation \(i\) in cluster \(j\), and let \(b_j\) denote the cluster’s random effect or random-effects vector. A GLMM specifies the conditional distribution of \(Y_{ij}\) given \(b_j\), and relates its conditional mean \(\mu_{ij}=E(Y_{ij}\mid b_j)\) to predictors through a link function:

$$ g(\mu_{ij})=\boldsymbol{x}_{ij}^{\mathsf T}\boldsymbol{\beta} +\boldsymbol{z}_{ij}^{\mathsf T}b_j. $$

Here, \(\boldsymbol{\beta}\) contains fixed-effect coefficients, \(\boldsymbol{x}_{ij}\) contains the fixed-effect predictors (including interaction terms), and \(\boldsymbol{z}_{ij}\) defines the random-effects contribution. The particular response distribution and link depend on the outcome. For a binary outcome with a logit link, fixed-effect coefficients describe changes in conditional log odds, given the random effects.

Suppose \(A\) is a categorical predictor with levels \(0,1,2\), where level 0 is the reference, and \(T\) is a continuous predictor. With treatment coding, an interaction model includes the terms \(A_1,A_2,T,A_1T,A_2T\), in addition to an intercept and any other specified predictors. Its interaction block is \((\beta_{1T},\beta_{2T})\), and the joint null is \(H_0:\beta_{1T}=\beta_{2T}=0\). Under this null, the slope for \(T\) is the same across all three levels of \(A\). The alternative allows at least one of the two slopes to differ.

The interaction’s interpretation depends on the link and coding. With a logit link, the interaction coefficients represent differences in slopes on the conditional log-odds scale. They are not, in general, differences in probabilities. Also, when an interaction is included, a main-effect coefficient is conditional on the other interacting predictor being at its reference level or numerical value zero. Centering a continuous predictor can make that reference value more meaningful, but does not change the underlying fitted relationship.

Worked Example: Interpreting a One-Coefficient Interaction

Suppose a hypothetical study records a binary outcome for patients in 18 clinics. Each patient has one record, and the model includes a clinic random intercept. A continuous exposure \(T\), measured in hours, interacts with a binary group indicator \(A\) (0 = reference group, 1 = comparison group). A fitted conditional logit model is

$$ \operatorname{logit}\{P(Y_{ij}=1\mid b_j)\} =\beta_0+b_j+\beta_A A_{ij}+\beta_T T_{ij} +\beta_{AT}A_{ij}T_{ij}. $$

Assume the estimated interaction coefficient is \(\hat{\beta}_{AT}=-0.60\). The conditional log-odds slope for \(T\) is \(\beta_T\) in the reference group and \(\beta_T+\beta_{AT}\) in the comparison group. Their difference is \(-0.60\) log-odds units per hour. Exponentiating gives \(\exp(-0.60)=0.5488\): the conditional odds ratio associated with a one-hour increase in \(T\) in the comparison group is about 0.549 times the corresponding conditional odds ratio in the reference group. This describes a ratio of odds ratios, not a probability difference.

Because this interaction has one coefficient, the individual test of \(H_0:\beta_{AT}=0\) is also the joint test of the entire interaction block. With a multi-coefficient interaction, the two tests would not be interchangeable.

Joint Wald Tests

A joint Wald test uses the estimated interaction coefficients and their estimated covariance matrix. Let \(\hat{\boldsymbol{\beta}}_I\) be the vector of \(q\) interaction estimates and let \(\widehat{\operatorname{Var}}(\hat{\boldsymbol{\beta}}_I)\) be its estimated covariance matrix. Under \(H_0:\boldsymbol{\beta}_I=\boldsymbol{0}\), and when the usual large-sample approximation is adequate, the Wald statistic is approximately chi-square with \(q\) degrees of freedom.

$$ W= \hat{\boldsymbol{\beta}}_I^{\mathsf T} \left\{\widehat{\operatorname{Var}}(\hat{\boldsymbol{\beta}}_I)\right\}^{-1} \hat{\boldsymbol{\beta}}_I \ \overset{\text{approx.}}{\sim}\ \chi^2_q. $$

The covariance matrix matters because interaction estimates can be correlated. The joint statistic measures the estimated vector’s distance from zero while accounting for its uncertainty and dependence among components. Replacing the covariance matrix with separate standard errors and adding squared individual statistics would generally omit covariance terms and give the wrong joint test.

Justification: For a regular fixed-effect model, the estimated coefficient vector is approximately multivariate normal: \(\hat{\boldsymbol{\beta}}_I\) has approximate mean \(\boldsymbol{\beta}_I\) and covariance matrix \(V\). Under the joint null, \(V^{-1/2}\hat{\boldsymbol{\beta}}_I\) is approximately a vector of \(q\) independent standard normal variables. The sum of their squared values is therefore approximately \(\chi^2_q\), and that sum is exactly the quadratic form \(W=\hat{\boldsymbol{\beta}}_I^{\mathsf T}V^{-1}\hat{\boldsymbol{\beta}}_I\). In practice \(V\) is estimated, so this reference law is an approximation.

Worked Example: A Two-Degree-of-Freedom Wald Test

In a hypothetical binary-outcome GLMM, a two-degree-of-freedom interaction has estimates \(\hat{\beta}_1=0.50\) and \(\hat{\beta}_2=-0.30\). Suppose their estimated covariance matrix is

$$ \widehat{\operatorname{Var}}(\hat{\boldsymbol{\beta}}_I)= \begin{pmatrix} 0.16&0.04\\ 0.04&0.09 \end{pmatrix}. $$

The joint null is \(H_0:\beta_1=\beta_2=0\). The determinant of the covariance matrix is \(0.16(0.09)-0.04^2=0.0128\), so its inverse is

$$ \begin{pmatrix} 0.16&0.04\\ 0.04&0.09 \end{pmatrix}^{-1} = \begin{pmatrix} 7.03125&-3.125\\ -3.125&12.5 \end{pmatrix}. $$

Substitution into the joint Wald statistic gives

$$ \begin{aligned} W &= \begin{pmatrix}0.50&-0.30\end{pmatrix} \begin{pmatrix}7.03125&-3.125\\-3.125&12.5\end{pmatrix} \begin{pmatrix}0.50\\-0.30\end{pmatrix}\\ &=7.03125(0.50)^2 +2(-3.125)(0.50)(-0.30) +12.5(-0.30)^2\\ &=1.7578125+0.9375+1.125=3.8203125. \end{aligned} $$

Using a \(\chi^2_2\) reference distribution, the upper-tail probability is \(P(\chi^2_2\geq 3.8203)=\exp(-3.8203/2)\approx0.1481\). At \(\alpha=0.05\), this test does not reject the joint null. This result is not evidence that both interaction coefficients are exactly zero; it indicates that this test does not find sufficient evidence against their joint null at the stated level.

For comparison, the individual Wald statistics are \(0.50^2/0.16=1.5625\) and \((-0.30)^2/0.09=1.000\). Their approximate two-sided \(z\)-test p-values are about 0.211 and 0.317, respectively. The joint test is a distinct calculation: it uses the off-diagonal covariance of 0.04 and tests both restrictions together.

Likelihood-Ratio Tests for an Interaction

A likelihood-ratio test compares a model containing the interaction with a nested model in which the entire interaction block is removed. The models must describe the same response data, and the reduced model must be obtained by imposing the interaction null while retaining the other intended terms. For the standard large-sample test of fixed effects, the models are fitted by maximum likelihood and the random-effects structure is held the same. The test statistic is twice the difference in maximized log likelihoods; under regularity conditions and the joint null, it is approximately chi-square with degrees of freedom equal to the number of interaction restrictions.

$$ G^2=2\{\ell_{\mathrm{full}}-\ell_{\mathrm{reduced}}\} \ \overset{\text{approx.}}{\sim}\ \chi^2_q. $$

For a \(q\)-coefficient interaction block, the likelihood-ratio test and the joint Wald test address the same null, but use different approximations to the evidence. They need not produce identical statistics or p-values, particularly when information is limited or estimates are imprecise. The degrees of freedom come from the number of restrictions, not the number of clusters or observations.

Worked Example: Testing a Group-by-Time Interaction by Likelihood Ratio

A hypothetical dataset contains binary outcomes for 600 patients in 30 clinics. Each patient belongs to one of three groups and has a measurement at one of two time points. A clinic random intercept is included in both models. With group 0 as the reference, the full model has two group-by-time interaction coefficients. Suppose the fitted full model has maximized log likelihood \(-121.4\), while the nested model without either interaction term has maximized log likelihood \(-124.0\).

The null hypothesis is that both interaction coefficients are zero. The likelihood-ratio statistic is

$$ G^2=2\{-121.4-(-124.0)\}=2(2.6)=5.2. $$

There are two restrictions, so the reference distribution is \(\chi^2_2\). Its upper-tail probability is \(\exp(-5.2/2)=\exp(-2.6)\approx0.07427\). At \(\alpha=0.05\), the result does not reject the joint null. The data do not provide sufficient evidence at that level that the group differences change over time on the model’s conditional log-odds scale. This conclusion concerns the interaction block; it does not establish that group or time has no association with the outcome.

Assumptions, Diagnostics, and Common Pitfalls

Neither a Wald test nor a likelihood-ratio test repairs a poorly specified GLMM. Interpret the reference p-value only after checking that the model and comparison are appropriate. Important conditions include:

  • Correct interaction specification. Include the lower-order terms associated with the interaction. For example, a model with \(A\times T\) ordinarily includes both \(A\) and \(T\), so the interaction represents a change in the association rather than replacing the main effects.
  • Appropriate conditional response model. The chosen response distribution and link should suit the outcome. The random-effects structure should represent the intended clustering, and observations should be conditionally independent given the specified random effects and fixed predictors.
  • Identifiable, estimable coefficients. Sparse cells, few clusters, near-separation, or strong predictor collinearity can make interaction estimates unstable. Check convergence messages, standard errors, estimated correlations, and whether fitted values or estimates change substantially under reasonable model checks.
  • Suitable reference approximation. The chi-square approximations are large-sample results under regularity conditions. A small number of clusters, sparse outcomes, or estimates near a boundary can make calibration poor. A non-significant result in such settings may reflect limited information rather than a well-supported absence of interaction.
  • Valid nested comparison. For a likelihood-ratio test, the two models must be fitted to the same observations and use the same response specification and random-effects structure. A change in likelihood caused by changing multiple model components cannot be attributed solely to the interaction.

A common mistake is to test the interaction by looking at each coefficient’s p-value and declaring the interaction absent if neither is below 0.05. That reasoning tests separate restrictions, not the joint hypothesis. Conversely, finding one small individual p-value does not by itself provide the specified joint test of a multi-parameter interaction. The correction is to state the full null and use a joint Wald test or a nested-model likelihood-ratio test with the appropriate degrees of freedom.

Another frequent error is to interpret a non-significant interaction as proof that effects are identical across groups. A test that does not reject \(H_0\) does not establish exact equality. If the scientific question concerns whether differences are small enough to be practically negligible, the analysis needs a prespecified equivalence margin and an appropriate inferential procedure; a conventional non-rejection alone does not answer that question.

The alternatives have different strengths. A joint Wald test uses one fitted model and the estimated covariance matrix, making it convenient for testing a specified coefficient block. A likelihood-ratio test compares nested fits and evaluates the loss of fit when the whole block is removed. Wald tests can be sensitive to the local curvature and uncertainty approximation; likelihood-ratio tests can also be poorly calibrated in small or sparse samples. Both usually rely on asymptotic chi-square calibration. Neither test tells the full pattern of group-specific slopes by itself, so estimated coefficients and their uncertainty should accompany the test when describing the fitted interaction.

Key takeaway: First define the interaction block and its joint null. Then use a test whose degrees of freedom equal the number of restrictions, check that the model and approximation are credible, and interpret the result on the link scale with the predictor units and coding made explicit.

Check Your Understanding

Answer each question, then compare your reasoning with the brief solutions.

  1. A four-level factor interacts with a continuous predictor using reference-level indicator coding. How many interaction coefficients are tested jointly?
  2. In a logit GLMM, why is an interaction coefficient not generally a probability difference?
  3. A joint Wald test gives \(W=6.0\) for a three-coefficient interaction. At \(\alpha=0.05\), should the null be rejected? Use the approximate \(\chi^2_3\) critical value 7.815.
  4. For a two-coefficient interaction, the full and reduced maximized log likelihoods are \(-80.0\) and \(-83.5\). Compute the likelihood-ratio statistic and state its degrees of freedom.

Solutions. (1) There are \(4-1=3\) interaction coefficients, so the joint test has 3 degrees of freedom. (2) The coefficient describes a change in the conditional log-odds slope (or contrast) under the logit link. Probabilities are a nonlinear transformation of log odds and also depend on the other predictors and random effects. (3) Since \(6.0<7.815\), do not reject the joint null at 0.05. This is insufficient evidence against the specified null, not proof that the coefficients are zero. (4) \(G^2=2\{-80.0-(-83.5)\}=7.0\), with 2 degrees of freedom because two interaction restrictions are tested.

Clinical Biostats