Random Slopes Guide¶

Random slopes let each group have its own relationship between a predictor and the outcome — rather than a single population slope shared across all groups. This page explains when to use them, how to specify them, and how to interpret the results.

When to add a random slope¶

Start with random intercepts only (groups=). Consider adding a random slope when:

Theory suggests it: subjects or groups plausibly differ in how much a predictor affects them (not just their baseline)
Model comparison favours it: a likelihood ratio test or AIC comparison of intercept-only vs slope model is significant
Residual plots show group-by-predictor interaction: conditional residuals fan out systematically with the predictor value

Rule of thumb: if you have fewer than ~10 observations per group for the predictor of interest, a random slope may be poorly identified. Check that the model converges and that the slope variance is not estimated near zero.

Syntax¶

Correlated intercept + slope¶

The default: the intercept and slope for each group are modelled as jointly Normal, allowing them to correlate (e.g. groups with a high baseline also tend to have a steeper slope).

import interlace

result = interlace.fit(
    "rt ~ condition",
    data=df,
    random=["(1 + condition | subject)"],
)

Equivalent lme4 R syntax: rt ~ condition + (1 + condition | subject)

Independent (uncorrelated) parameterisation¶

Use || to force the intercept and slope variance to be estimated independently, with covariance fixed to zero. This is more parsimonious and useful when correlation cannot be estimated reliably (few groups, sparse data per group).

result = interlace.fit(
    "rt ~ condition",
    data=df,
    random=["(1 + condition || subject)"],
)

Equivalent lme4 R syntax: rt ~ condition + (1 + condition || subject)

Mixing slope and intercept-only terms¶

Use random= for all terms when combining random slopes for one factor with random intercepts for another:

result = interlace.fit(
    "rt ~ condition",
    data=df,
    random=[
        "(1 + condition | subject)",  # by-subject slope
        "(1 | item)",                 # item intercept only
    ],
)

Interpreting results¶

BLUPs (random_effects)¶

When a model includes random slopes, result.random_effects[group] returns a DataFrame rather than a Series — one column per random effect term:

print(result.random_effects["subject"])
#             (Intercept)  condition
# subject_01       -12.3       0.42
# subject_02         8.7      -0.31
# subject_03         2.1       0.05
# ...

The (Intercept) column is the by-subject deviation from the grand mean. The condition column is the by-subject deviation from the population slope for condition. A subject with condition = 0.42 responds more strongly to condition than average.

Variance components and covariance matrix¶

# Variance of random intercepts and slopes
print(result.variance_components)
# {'subject': {'(Intercept)': 45.2, 'condition': 3.1}, 'residual': 12.8}

# Full covariance matrix (intercept-slope correlation included)
print(result.varcov)
#                 (Intercept)  condition
# (Intercept)         45.2       -8.4
# condition            -8.4        3.1

The off-diagonal of varcov gives the intercept-slope covariance. A negative value means groups with a higher baseline tend to show a smaller (or negative) slope.

Comparing to the intercept-only model¶

Fit both models with method="ML" and compare:

m_intercept = interlace.fit("rt ~ condition", data=df,
                             groups=["subject", "item"], method="ML")

m_slopes = interlace.fit("rt ~ condition", data=df,
                          random=["(1 + condition | subject)", "(1 | item)"],
                          method="ML")

import scipy.stats
# Correlated slope adds 2 parameters (slope variance + covariance)
lrt = 2 * (m_slopes.llf - m_intercept.llf)
p   = scipy.stats.chi2.sf(lrt, df=2)
print(f"LRT χ²(2) = {lrt:.2f}, p = {p:.4f}")

If p < 0.05, the slope variance is statistically meaningful. See the Model Comparison Guide for the full LRT workflow.

Common issues¶

Singular fit / convergence warning¶

If the optimizer warns about a singular fit or near-zero slope variance, the random slope may be over-parameterised for your data:

Try the independent (||) parameterisation first
Check that you have enough observations per group (≥ 5–10 per group for the predictor)
Compare AIC with the intercept-only model; if they are similar, prefer the simpler model

Switching from `groups=` to `random=`¶

# Before
result = interlace.fit("y ~ x", data=df, groups=["g1", "g2"])

# After — equivalent intercept-only model using random=
result = interlace.fit(
    "y ~ x",
    data=df,
    random=["(1 | g1)", "(1 | g2)"],
)

Both produce identical results. Use whichever is clearer for your use case.

Uncertainty in BLUPs¶

BLUPs are point estimates — each group’s deviation from the population mean. The random_effects_se property and random_effects_ci() method expose the posterior standard errors and normal-approximation confidence intervals for those estimates.

Standard errors¶

se = result.random_effects_se
# Intercept-only model → pd.Series indexed by group
print(se["subject"])
# subject_01    4.32
# subject_02    3.87
# subject_03    5.01
# ...

# Random-slope model → pd.DataFrame, one column per term
print(se["subject"])
#             (Intercept)  condition
# subject_01         4.32       0.61
# subject_02         3.87       0.54

Confidence intervals¶

# 95 % CIs (default)
ci = result.random_effects_ci()
print(ci["subject"])
#             lower    upper
# subject_01  -20.8     -3.8
# subject_02   -2.1     17.3
# ...

# 90 % CIs
ci_90 = result.random_effects_ci(level=0.90)

For random-slope models random_effects_ci() returns a DataFrame with a MultiIndex column: (term, "lower") and (term, "upper") pairs.

Caterpillar plot¶

A caterpillar plot orders groups by their BLUP and overlays the CI — a quick visual check for which groups stand out from the population mean:

import pandas as pd
from plotnine import ggplot, aes, geom_point, geom_errorbar, geom_hline, coord_flip, theme_bw

blups = result.random_effects["subject"]
ci    = result.random_effects_ci()["subject"]

caterpillar = pd.DataFrame({
    "group":  blups.index,
    "blup":   blups.values,
    "lower":  ci["lower"].values,
    "upper":  ci["upper"].values,
}).sort_values("blup").assign(rank=lambda d: range(len(d)))

(
    ggplot(caterpillar, aes(x="rank", y="blup"))
    + geom_errorbar(aes(ymin="lower", ymax="upper"), width=0.3, alpha=0.5)
    + geom_point(size=2)
    + geom_hline(yintercept=0, linetype="dashed", color="grey")
    + coord_flip()
    + theme_bw()
)

Groups whose CI excludes zero differ reliably from the population mean. Note that these are normal-approximation CIs — treat them as indicative rather than exact for small group counts or boundary variance estimates.