TopolearnStudents in schools, patients over time and users across sessions are not independent rows.
When data are grouped, observations within a group share something unobserved. A linear mixed model represents that with random effects: group-specific deviations drawn from a normal distribution, estimated alongside the usual fixed effects.
smf.mixedlm("score ~ hours", data=df, groups=df["school"]).fit() gives each school its own random intercept. Add random slopes with re_formula="~hours", letting the effect of hours vary by school too.
Reading the results:
- fixed effects (
res.fe_params) are population-level effects, interpreted like OLS coefficients, Group Varin the summary is the variance of the random intercepts; compare it with the residual variance (res.scale) to see how much variation is between groups (the intraclass correlation),res.random_effectsgives each group's estimated deviation, shrunk towards zero for small groups (partial pooling).
Why not dummies per group? Mixed models generalise to new groups, handle many small groups gracefully, and let you include group-level predictors that would be collinear with group dummies.
Fitting uses REML by default (reml=False for comparing fixed-effect structures by likelihood). Convergence warnings are common with complex random effects; simplify the random structure or rescale predictors. For non-normal outcomes, BinomialBayesMixedGLM exists, and GEE is an alternative when you want population-averaged effects.
Example
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
rng = np.random.default_rng(11)
school = np.repeat(np.arange(30), 20) # 30 schools x 20 students
hours = rng.uniform(0, 10, 600)
school_int = rng.normal(0, 4, 30)[school] # schools differ in baseline
df = pd.DataFrame({"school": school, "hours": hours,
"score": 50 + 2 * hours + school_int + rng.normal(0, 3, 600)})
res = smf.mixedlm("score ~ hours", data=df, groups=df["school"]).fit()
print(res.fe_params.round(2).to_dict()) # {'Intercept': 50.6, 'hours': 1.99}
group_var = res.cov_re.iloc[0, 0]
print(round(group_var, 1), round(res.scale, 1)) # 29.5 9.7
print(round(group_var / (group_var + res.scale), 2)) # 0.75 (intraclass correlation)
print(round(res.random_effects[0].iloc[0], 2)) # 2.33 (school 0 above average)Remember
- Random effects model group-level variation; fixed effects are population effects
- Random intercepts via
groups=, random slopes viare_formula= - Group estimates are shrunk towards the mean (partial pooling)
Watch out for
- Ignoring grouping and treating repeated measurements as independent rows.