Random Slopes at Multiple Nested Levels

Why Allow Slopes to Vary at More Than One Level?

In a nested study, observations may be grouped within units that are themselves grouped within larger units. For example, repeated measurements may be collected from patients within clinics, with clinics nested within regions. A predictor’s association with the outcome may differ from one clinic to another, and clinics in the same region may also share a broader pattern of variation. A model with random slopes at both levels represents these two sources of unit-to-unit variation in the predictor effect.

The essential distinction is between a fixed slope, which describes the average link-scale association across the population of units, and random slopes, which describe deviations from that average at particular units. In a GLMM, these combine on the linear predictor scale. When the link is nonlinear, the corresponding response-scale change also depends on the intercept and on the predictor value.

Model and Assumptions

Index observations by \(i\), the level-2 units within a level-3 unit by \(j\), and the level-3 units by \(k\). For example, \(i\) could index patients, \(j\) clinics, and \(k\) regions. Let \(x_{ijk}\) be a predictor whose effect is allowed to vary at both levels. A random-intercept-and-slope model at each level is

g\{\mu_{ijk}\} = \beta_0+\beta_1x_{ijk} +b_{0k}+b_{1k}x_{ijk} +u_{0jk}+u_{1jk}x_{ijk}, \qquad \mu_{ijk}=E(Y_{ijk}\mid b_{0k},b_{1k},u_{0jk},u_{1jk}).

Here \(g\) is the link function, \(\beta_0\) and \(\beta_1\) are fixed effects, \(b_{0k}\) and \(b_{1k}\) are the level-3 random intercept and slope, and \(u_{0jk}\) and \(u_{1jk}\) are the level-2 random intercept and slope. Because the level-2 units are nested, each \(j\) refers to a unit within a particular \(k\). For a Bernoulli outcome, for example, \(g\) may be the logit link; for a normally distributed outcome with an identity link, the linear predictor is the conditional mean.

Definition: A random slope at a level is a random coefficient multiplying a predictor at that level. In this model, the slope for an observation in level-2 unit \(j\) within level-3 unit \(k\), on the link scale, is \(\beta_1+b_{1k}+u_{1jk}\). The random intercepts shift the linear predictor but do not directly change this link-scale slope.

A standard specification assumes that random-effect vectors at different level-3 units are independent and that, conditional on their level-3 unit, the level-2 random-effect vectors are independent across level-2 units. The level-2 and level-3 random-effect vectors are assumed independent of each other in this formulation. Within a level, an intercept and a slope may be correlated:

\begin{pmatrix}b_{0k}\\b_{1k}\end{pmatrix} \sim N\left[ \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix} \tau_{00}^{(3)} & \tau_{01}^{(3)}\\ \tau_{01}^{(3)} & \tau_{11}^{(3)} \end{pmatrix} \right], \qquad \begin{pmatrix}u_{0jk}\\u_{1jk}\end{pmatrix} \sim N\left[ \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix} \tau_{00}^{(2)} & \tau_{01}^{(2)}\\ \tau_{01}^{(2)} & \tau_{11}^{(2)} \end{pmatrix} \right].

The covariance matrices must be positive semidefinite. The conditional response distribution must also be appropriate for the outcome, and observations are assumed conditionally independent given the random effects and fixed effects, unless the response model explicitly includes further dependence. The grouping structure must genuinely be nested, and there must be enough information within units to estimate slope variation. In particular, a predictor that never varies within a level-2 unit cannot identify a level-2 random slope from within-unit comparisons.

The value \(x=0\) defines the reference point for the random intercepts. Centering \(x\) at a substantively meaningful value can make intercepts easier to interpret and may improve numerical stability, but it does not remove the need for adequate predictor variation at each level.

Worked Example: Combining Two Random Slopes

Suppose a hypothetical continuous-outcome model uses the identity link, and the fitted fixed slope for a treatment-intensity predictor is \(\beta_1=2.0\) outcome units per intensity unit. For one region, the estimated random slope is \(b_{1k}=-0.4\); for one clinic within it, the estimated random slope is \(u_{1jk}=0.7\). Find the clinic’s conditional slope.

The model’s link-scale slope in this clinic is

\beta_1+b_{1k}+u_{1jk} =2.0-0.4+0.7=2.3.

With an identity link, this is also the change in the conditional mean outcome for a one-unit increase in intensity: 2.3 outcome units per intensity unit for this clinic, holding the model’s other terms fixed. The region’s negative deviation reduces the population slope, while the clinic’s positive deviation more than offsets that reduction. These conditional estimates describe variation represented by the model; they do not establish a causal effect.

Interpreting Slopes in a Nonlinear GLMM

For a nonlinear link, the combined coefficient remains a slope on the link scale, but it is not generally a constant change in the response mean. The response mean is \(h(\eta_{ijk})\), where \(h=g^{-1}\) and \(\eta_{ijk}\) is the linear predictor. For a continuous predictor, the conditional response-scale derivative is

\frac{\partial \mu_{ijk}}{\partial x_{ijk}} = h'(\eta_{ijk})\bigl(\beta_1+b_{1k}+u_{1jk}\bigr).

Thus the random slopes govern unit-specific changes on the link scale, while the response-scale effect also depends on the current linear predictor. For a binary outcome with a logit link, the link-scale slope is a change in log odds per predictor unit. The corresponding change in probability depends on both the slope and the starting probability.

Worked Example: A Clinic-Specific Log-Odds Slope

Consider a hypothetical Bernoulli GLMM with a logit link. Let \(\beta_0=-1.0\), \(\beta_1=0.6\), and, for a particular region and clinic, let \(b_{0k}=0.2\), \(b_{1k}=-0.3\), \(u_{0jk}=-0.1\), and \(u_{1jk}=0.4\). Calculate the conditional probabilities at \(x=0\) and \(x=1\), and compare them.

The clinic-specific intercept and slope on the log-odds scale are

\beta_0+b_{0k}+u_{0jk}=-1.0+0.2-0.1=-0.9, \qquad \beta_1+b_{1k}+u_{1jk}=0.6-0.3+0.4=0.7.

At \(x=0\), the linear predictor is \(-0.9\), so the conditional probability is

P(Y=1\mid x=0,\text{ random effects}) = \frac{1}{1+\exp(0.9)} \approx 0.2891.

At \(x=1\), the linear predictor is \(-0.9+0.7=-0.2\), giving

P(Y=1\mid x=1,\text{ random effects}) = \frac{1}{1+\exp(0.2)} \approx 0.4502.

Using the unrounded probabilities, the estimated conditional probability difference is approximately \(0.450166-0.289050=0.1611\), or 16.11 percentage points for this one-unit predictor change in this clinic. The random effects yield a clinic-specific log-odds slope of 0.7, not a probability increase of 0.7. The probability difference applies to these specified covariates and random effects; it need not be the same at other starting values or in other clinics.

Distribution of the Combined Slope

The random slopes at the two levels add for a particular level-2 unit. Under the stated independence between levels, their sum has variance equal to the sum of their variances. This gives a direct interpretation of how much slope variation the model allows across the lowest-level units.

Justification: The conditional link-scale slope is \(\beta_1+b_{1k}+u_{1jk}\). If \(b_{1k}\) and \(u_{1jk}\) are independent, mean-zero normal variables with variances \(\tau_{11}^{(3)}\) and \(\tau_{11}^{(2)}\), then their sum is normal, with mean zero and variance \(\tau_{11}^{(3)}+\tau_{11}^{(2)}\). Therefore, the combined slope has mean \(\beta_1\) and variance \(\tau_{11}^{(3)}+\tau_{11}^{(2)}\). This statement concerns the distribution across units, not uncertainty in an estimated slope for one particular unit.

Worked Example: Interpreting Overall Slope Variation

Suppose a hypothetical fitted model has fixed slope \(\beta_1=0.80\), level-3 random-slope standard deviation \(0.25\), and level-2 random-slope standard deviation \(0.40\). Assume the model’s random slopes are independent across levels and normally distributed as specified above. Find the standard deviation of the combined slope and the approximate probability that a randomly selected level-2 unit has a negative link-scale slope.

The two random-slope variances are \(0.25^2=0.0625\) and \(0.40^2=0.1600\). By independence, the combined slope variance and standard deviation are

0.0625+0.1600=0.2225, \qquad \sqrt{0.2225}\approx 0.4717.

The combined slope is therefore normally distributed with mean \(0.80\) and standard deviation \(0.4717\). Standardizing zero gives

P(\text{combined slope}<0) = P\left(Z<\frac{0-0.80}{0.4717}\right) = P(Z<-1.696) \approx 0.045.

Under this fitted distributional model, about 4.5% of level-2 units are predicted to have a negative link-scale slope. This is a model-based population proportion, not a probability that a fixed, known clinic’s true slope is negative. The calculation also relies on the independence and normality assumptions; a different covariance structure would change the variance calculation.

Diagnostics, Alternatives, and Limitations

Multiple-level random slopes add parameters and require information at each level. Check whether the predictor varies within enough units at both levels, whether there are enough level-3 units and enough level-2 units per level-3 unit, and whether the fitted model is numerically stable. These are design and estimation questions, not merely software warnings.

Convergence warnings or a singular fit. A variance estimated at or near zero, or a covariance matrix on its boundary, may indicate that the data do not support the full random-effects structure. It can also signal scaling or optimization problems. Inspect the estimates and sensitivity to reasonable predictor scaling rather than treating a warning as proof that the model is wrong.
Weak within-unit predictor variation. If few units contain observations at multiple predictor values, the corresponding slope variances can be poorly determined. Random-intercept variation cannot substitute for missing information about within-unit slopes.
Distributional mismatch. Severe outliers or patterns in estimated random effects may raise concerns about the Gaussian random-effects assumption. Such diagnostics are imperfect, especially with few higher-level units, so they should be interpreted cautiously.
Conditional response-model problems. Residual patterns or poor predictive calibration may indicate that the response distribution, link, fixed-effects structure, or conditional independence assumption is inadequate. Random slopes do not automatically fix these issues.

The main alternatives answer different questions. A random-intercept-only model allows units to differ in baseline level but constrains the predictor’s link-scale slope to be common. A model with a random slope at only one nested level permits slope variation there but not at the other. Treating every unit’s slope as an unrelated fixed effect may use many parameters and does not provide the same partial pooling as a random-effects model. The multiple-level random-slope model is appropriate when the hierarchy and scientific question support variation at both levels and the data contain enough information to estimate it.

A common mistake is to add the two random slopes but interpret the result as a response-scale probability change. The sum \(\beta_1+b_{1k}+u_{1jk}\) is a link-scale slope. For a logit model it is a log-odds change per predictor unit; the probability change must be calculated from the linear predictor at specified covariate values, as in the worked example. A second pitfall is to infer that the random slope at a higher level is identified simply because there are many observations overall. Information about higher-level variation depends importantly on the number of higher-level units and the distribution of predictor values across them.

Key takeaway: In a nested GLMM with random slopes at two levels, the unit-specific link-scale slope is the fixed slope plus one random slope from each level. Interpret that sum on the link scale; translate to the response scale only with the inverse link and the unit’s full linear predictor.

Check Your Understanding

Answer each question using the model and assumptions stated in this tutorial. Brief solutions follow.

If \(\beta_1=1.2\), \(b_{1k}=0.3\), and \(u_{1jk}=-0.5\), what is the unit-specific link-scale slope?
In a logit model, what does a unit-specific slope of \(-0.4\) mean?
Under independence across levels, what is the variance of the combined slope if the level-3 and level-2 random-slope variances are \(0.09\) and \(0.16\)?
Why can a level-2 random slope be difficult to estimate if the predictor is constant within each level-2 unit?
Does a link-scale slope of \(0.7\) in a logistic model imply a probability increase of 0.7 for a one-unit predictor increase? Explain.

Solutions. (1) \(1.2+0.3-0.5=1.0\) link-scale unit per predictor unit. (2) For that unit, a one-unit predictor increase is associated with a decrease of 0.4 in log odds, conditional on the other model terms; the corresponding probability change depends on the starting linear predictor. (3) Independence gives \(0.09+0.16=0.25\). (4) Without within-unit predictor variation, the data provide no within-unit contrasts to distinguish slope variation at that level from other sources of variation. (5) No. The value 0.7 is a log-odds change; the probability change must be computed through the inverse-logit at specified covariate and random-effect values.