Cross-Level Interactions in Multilevel Models

Why Cross-Level Interactions Matter

An individual-level association need not be the same in every group. For example, the association between a patient’s treatment exposure and an outcome may differ across clinics with different staffing levels. A cross-level interaction represents this possibility directly: a predictor measured for individuals is interacted with a predictor measured for their group.

The interaction is not merely a product term added for convenience. It expresses a specific model claim: the association between the individual-level predictor and the outcome changes systematically with the group-level predictor. In a generalized linear mixed model (GLMM), that claim is first expressed on the model’s link scale. Its implications for the response scale depend on the link and outcome distribution.

Model Specification and Assumptions

Let $j=1,\ldots,J$ index groups and $i=1,\ldots,n_j$ index observations within group $j$. Let $Y_{ij}$ be the response, $x_{ij}$ an individual-level predictor, and $z_j$ a group-level predictor. For interpretation, suppose $x_{ij}$ has been centered within groups and $z_j$ has been centered around a meaningful overall reference. Other predictors can be included, but are omitted here to make the interaction clear.

Definition: A cross-level interaction occurs when a predictor measured at one hierarchy level modifies the association of a predictor measured at another level. In the model below, $x_{ij}z_j$ is the product of an individual-level predictor and a group-level predictor: $$ g\!\left(E[Y_{ij}\mid u_{0j},u_{1j}]\right) = \beta_0+\beta_x x_{ij}+\beta_z z_j+\beta_{xz}x_{ij}z_j+u_{0j}+u_{1j}x_{ij}. $$ Here $g$ is the link function, $\beta_{xz}$ is the fixed cross-level interaction coefficient, and $u_{0j}$ and $u_{1j}$ are group-specific random intercept and slope deviations. The random effects are commonly modeled jointly as mean-zero multivariate normal variables with covariance matrix $\Sigma$, independent across groups and independent of the predictors. Conditional on the random effects and predictors, responses are assumed independent and to follow the specified response distribution.

The response distribution and link must suit the outcome: for example, a normal distribution with identity link for a continuous outcome, a Bernoulli distribution with logit link for a binary outcome, or a Poisson distribution with log link for a count. The linear predictor specifies the conditional mean, given the group’s random effects. A model may omit $u_{1j}$ if random slopes are not included, but doing so assumes there is no residual group-to-group variation in the $x$ slope beyond the fixed effects.

The group-level predictor must genuinely vary between groups, and the individual-level predictor must vary sufficiently within groups to estimate its within-group association. To estimate moderation, there must also be adequate information across the range of $z_j$; many observations within only a few groups do not substitute for information from more groups. As in the earlier material on centering predictors within nested groups, centering helps make the lower-order coefficients correspond to clear reference values. The conditional distribution, random-effects structure, and link should be defensible for the data at hand.

What the Interaction Coefficient Means

For a fixed group-level value $z_j=z$, the part of the linear predictor that changes with $x_{ij}$ has slope $\beta_x+\beta_{xz}z+u_{1j}$. Thus, $\beta_{xz}$ is the change in the fixed $x$ slope for a one-unit increase in $z$, measured on the link scale. It is not, by itself, the slope of $x$, nor is it generally the change in the response mean per unit of $z$.

\begin{aligned} \eta_{ij}(x+1,z)-\eta_{ij}(x,z) &= \beta_x+\beta_{xz}z+u_{1j},\\ \{\eta_{ij}(x+1,z+1)-\eta_{ij}(x,z+1)\} &-\{\eta_{ij}(x+1,z)-\eta_{ij}(x,z)\} = \beta_{xz}, \end{aligned}

The second line is a direct justification: compare the one-unit $x$ contrast at $z+1$ with that at $z$. The intercept terms and random intercept cancel, as does the random slope; the remaining difference is $\beta_{xz}$. This is moderation of the fixed link-scale slope. For a particular group, the slope also includes its random deviation $u_{1j}$.

When $z=0$, the fixed $x$ slope is $\beta_x$; when $x=0$, the fixed group-level contrast per unit of $z$ is $\beta_z$. These lower-order coefficients are conditional on the other predictor being at its reference value. Centering changes the reference values and therefore the interpretation of the lower-order coefficients, even though a consistent linear recentering does not change the fitted interaction coefficient.

Worked Example: A Continuous Outcome in Clinics

Hypothetical setting. Suppose a Gaussian mixed model relates a continuous symptom score to each patient’s weekly treatment hours $x_{ij}$, centered at that patient’s clinic mean, and clinic staffing quality $z_j$, measured in standard-deviation units and centered at its overall mean. Lower scores indicate fewer symptoms. Assume the fitted fixed coefficients are $\hat\beta_x=-0.40$, $\hat\beta_z=0.30$, and $\hat\beta_{xz}=0.25$, with an identity link. For illustration, consider a clinic whose random slope deviation is $u_{1j}=-0.05$.

At average staffing quality, $z_j=0$, the group-specific conditional slope for treatment hours is $\hat\beta_x+\hat\beta_{xz}(0)+u_{1j}=-0.40+0-0.05=-0.45$ symptom-score points per additional weekly treatment hour. At staffing quality one standard deviation above average, it is $-0.40+(0.25)(1)-0.05=-0.20$ points per hour. At one standard deviation below average, it is $-0.40+(0.25)(-1)-0.05=-0.70$ points per hour.

The fixed interaction estimate says that each one-standard-deviation increase in staffing quality makes the fixed treatment-hours slope 0.25 score points per hour less negative. The random slope deviation shifts all three slopes for this particular clinic by the same $-0.05$. This model describes conditional associations; these hypothetical estimates alone do not establish that increasing treatment hours or staffing quality causes symptom changes.

Interpreting Interactions for Non-Gaussian Outcomes

With a nonlinear link, the interaction has a link-scale interpretation, not necessarily a constant difference in response-scale slopes. For a logit link, a one-unit increase in $x$ multiplies the conditional odds by $\exp(\beta_x+\beta_{xz}z+u_{1j})$. For a log link, it multiplies the conditional mean or rate by that same exponential expression. These are conditional interpretations, given the random effects. They should not be confused with marginal effects averaged over the random-effects distribution.

Worked Example: A Binary Outcome and Odds Ratios

Hypothetical setting. Consider a Bernoulli GLMM for whether a patient adheres to a care plan. Let $x_{ij}$ be counseling minutes in ten-minute units, centered within clinic, and let $z_j$ be a centered clinic-level support score in standard-deviation units. Suppose the fitted logit model has $\hat\beta_0=-1.00$, $\hat\beta_x=-0.35$, $\hat\beta_z=0.40$, and $\hat\beta_{xz}=0.50$. For a clear fixed-effect calculation, set the random effects to zero.

At $z= -1$, the log-odds slope for a ten-minute increase in counseling is $-0.35+(0.50)(-1)=-0.85$. The conditional odds ratio is $\exp(-0.85)=0.427$. At $z=1$, the slope is $-0.35+(0.50)(1)=0.15$, giving odds ratio $\exp(0.15)=1.162$. Thus the estimated counseling association differs across the two support levels on the odds scale.

To see why these are not fixed probability changes, take $x=0$. At $z=-1$, the linear predictor is $-1.00+0.40(-1)=-1.40$, so the conditional probability is $1/(1+\exp(1.40))=0.198$. At $x=1$, it is $-1.40-0.85=-2.25$, giving probability $1/(1+\exp(2.25))=0.095$. At $z=1$, the corresponding linear predictors are $-1.00+0.40=-0.60$ at $x=0$ and $-0.60+0.15=-0.45$ at $x=1$, with probabilities $0.354$ and $0.389$. The probabilities change by different amounts because the inverse-logit transformation is nonlinear. The calculations condition on zero random effects and do not represent probabilities averaged over clinics.

Worked Example: A Count Outcome and Rate Ratios

Hypothetical setting. Suppose a Poisson GLMM models the number of follow-up calls received by patients, using the log of each patient’s observation time as an offset. Let $x_{ij}$ be the number of prior contacts, centered within practice, and let $z_j$ be a practice-level workload measure in units of ten additional patients per staff member, centered at its mean. Suppose the fixed estimates are $\hat\beta_x=0.08$ and $\hat\beta_{xz}=-0.12$, with a log link. Set random effects to zero for the fixed-effect comparison.

At average workload, $z=0$, a one-unit increase in prior contacts changes the log rate by $0.08$, so the conditional rate ratio is $\exp(0.08)=1.083$, an estimated 8.3% higher rate. At workload two units above average, the slope is $0.08+(-0.12)(2)=-0.16$, and the rate ratio is $\exp(-0.16)=0.852$, or an estimated 14.8% lower rate per additional prior contact. The interaction is a change of $-0.12$ in the log-rate slope per ten-patient workload unit. It does not say that workload changes the count by a fixed number, and the offset means these interpretations concern rates per unit of observation time.

Assessing Evidence and Diagnosing Problems

A common inferential target is $H_0:\beta_{xz}=0$, against an alternative that the interaction is nonzero. A Wald test divides the estimated interaction by its standard error and compares the resulting statistic with an approximate standard normal reference distribution. A likelihood-ratio comparison can instead compare models with and without the fixed interaction, retaining the same response data and random-effects structure; under regularity conditions, twice the log-likelihood difference is approximately chi-squared with one degree of freedom. These are large-sample approximations, so uncertainty can be substantial when the number of groups is small or the predictor range is limited. The interaction estimate and its confidence interval should be reported in context, not replaced by a significance label.

Before interpreting a fitted interaction, inspect whether there is meaningful variation in $x$ within groups and in $z$ across groups, and whether groups cover overlapping ranges of the predictors. If high values of $z$ occur only in groups with unusually high or low $x$, estimates may rely heavily on extrapolation. Examine fitted response patterns over the observed $x,z$ combinations, residual behavior appropriate to the outcome, convergence messages, and the estimated random-effects variances and correlations. A singular or unstable random-effects fit, sparse outcome information in some groups, or strong sensitivity to a few groups can make interaction estimates unreliable.

The main alternatives answer different questions. Omitting $x z$ imposes a common fixed $x$ slope across values of $z$; it is more parsimonious but rules out systematic moderation. Including a random slope without $x z$ allows slopes to vary across groups without explaining that variation through $z$. Including both allows a systematic average relationship with $z$ and residual group-to-group slope variation. Fitting separate models for subgroups may seem intuitive, but it can discard information, complicate comparisons, and fail to represent a continuous moderator smoothly. These alternatives are not interchangeable: choose based on the scientific question, design, and available information.

Common Mistakes and Limitations

Reading $\beta_x$ as the slope everywhere. With an interaction, the fixed slope at $z$ is $\beta_x+\beta_{xz}z$, with a further group-specific deviation if a random slope is present.
Calling a link-scale interaction a constant probability difference. In a logistic model, the interaction changes the log-odds slope; the resulting probability difference depends on the starting probability.
Confusing systematic moderation with random slope variation. A random slope describes residual variation among groups; the product term describes how the average slope changes with the measured group-level predictor.
Interpreting the product without its lower-order terms. The model should ordinarily include $x$ and $z$ along with $xz$, so the interaction has a coherent conditional interpretation.

Even a well-fitting interaction model does not guarantee that the chosen functional form is adequate. A linear change in slope over $z$ may be too restrictive; sparse groups or limited overlap can make the interaction imprecise; and unmeasured group-level differences can complicate substantive interpretation. Model-based conditional associations should not be presented as causal effects without a design and assumptions that justify causal claims.

Key takeaway: In a cross-level interaction, $\beta_{xz}$ measures how the fixed association of an individual-level predictor changes with a group-level predictor on the GLMM’s link scale. Translate that result through the link before describing response-scale implications.

Check Your Understanding

Try each question, then check the brief solutions below.

In a log-link model, $\hat\beta_x=0.20$ and $\hat\beta_{xz}=-0.10$. What is the fixed rate ratio for a one-unit increase in $x$ when $z=2$?
For a logistic model, explain what $\beta_{xz}=0$ says and what it does not say about probabilities.
Why does a random slope alone not answer the same question as a cross-level interaction?
Name two design features that can make a cross-level interaction difficult to estimate reliably.

Solutions. (1) The log-rate slope is $0.20+(-0.10)(2)=0$, so the rate ratio is $\exp(0)=1$. (2) It says the fixed $x$ slope on the log-odds scale does not change with $z$; it does not imply a constant probability difference or that probabilities are equal across groups. (3) A random slope permits unexplained group-to-group slope variation, whereas the interaction models systematic average slope change associated with $z$. Both terms can be included. (4) Examples include too few groups, little variation in $z$, little within-group variation in $x$, and poor overlap in predictor ranges across groups.