By the end of today you can

• recognise repeated-measures and nested data;
• explain, with a simulation, why ordinary regression gets their uncertainty wrong;
• fit random intercept and random slope models with statsmodels;
• interpret fixed effects, random effects, and the intraclass correlation;
• test whether an effect changes over time, and compare models;
• model a count outcome with a generalised estimating equation;
• report a mixed-effects model in a thesis.

Two kinds of structure

Observations within a group are more alike than observations from different groups.

Less information than the rows suggest

A student who scores high this semester is likely to score high next semester.

Students of one supervisor share that supervisor’s style.

Methods that ignore this get the uncertainty wrong, sometimes badly.

A simulated training course with no effect

def simulate_clustered():
    supervisor = np.repeat(np.arange(20), 10)
    trained = np.repeat(np.tile([0, 1], 10), 10)   # no real effect
    style = np.random.normal(scale=5, size=20)[supervisor] # shared style
    wellbeing = 60 + style + np.random.normal(scale=10, size=200)
    # ... test "trained" with OLS and with MixedLM

20 supervisors, 10 students each, half the supervisors trained. Repeated 500 times.

Ignoring the grouping creates false alarms

The course went to 20 supervisors, not 200 independent students.

The mixed model stays close to 5–8%, recognizing that there are only 20 supervisors.

A study of sleep deprivation

sleepstudy.head()

Reaction times of 18 people over 10 days of sleep restriction, one measurement per person per day.

180 rows, but only 18 people.

Everyone gets slower, but differently

People start at different levels (intercepts) and slow at different rates (slopes).

Random intercepts

sleep_ri = smf.mixedlm("Reaction ~ Days", sleepstudy, groups=sleepstudy["Subject"]).fit()
print(sleep_ri.summary())

groups=sleepstudy["Subject"]: one average line, but a separate starting level for each person.

Fixed and random effects

The intraclass correlation

\[ \text{ICC} = \frac{\text{variance between groups}}{\text{variance between groups} + \text{variance within groups}} \]

var_between = sleep_ri.cov_re.iloc[0, 0]
var_within = sleep_ri.scale
icc = var_between / (var_between + var_within)

About 59% of the variation is stable differences between people.

ICC 0: the grouping does not matter. The higher it is, the more wrong an ordinary regression becomes.

Random slopes

sleep_rs = smf.mixedlm("Reaction ~ Days", sleepstudy, groups=sleepstudy["Subject"], re_formula="~Days").fit()

re_formula="~Days": a separate intercept and slope for each person.

Slopes vary with SD ≈ 5.9 ms per day. Ignoring it made the average too confident.

Is the random slope worth it?

lr_stat = -2 * (sleep_ri.llf - sleep_rs.llf)
p_val = stats.chi2.sf(lr_stat, df=2)

A likelihood ratio test compares the two models. The p-value is tiny (\(p < 0.001\)): random slopes improve the model clearly.

AIC agrees: it balances fit against complexity, and lower is better.

The wellbeing panel

panel = (
    semesters
    .merge(students, on="student_id")
    .assign(time=lambda d: d["semester"] - 1)
)

Up to four records per student, 2326 in all. time is 0 in semester 1.

Students differ far more than they change

Thirty random students: each stays roughly at their own level.

The ICC for wellbeing

wb_null = smf.mixedlm("wellbeing ~ 1", panel, groups=panel["student_id"]).fit()
icc = wb_null.cov_re.iloc[0, 0] / (wb_null.cov_re.iloc[0, 0] + wb_null.scale)

About 78% of all the variation in wellbeing lies between students.

The design effect

\[ \text{design effect} = 1 + (m - 1) \times \text{ICC} \]

2326 records carry about as much information as 716 independent observations.

Does wellbeing decline? The wrong way

ols_fit = smf.ols("wellbeing ~ time", panel).fit()

Slope -0.38, standard error 0.22: not significant (\(p = 0.55\)).

The right way: a growth model

wb_growth = smf.mixedlm("wellbeing ~ time", panel, groups=panel["student_id"], re_formula="~time").fit()

The ordinary regression puts the large differences between students into its error, drowning the small change within each student.

The error can go either way

The direction depends on the data. The error is always there.

Confidence intervals and inference

The time interval runs from -0.65 to -0.22: entirely below zero. Wellbeing declines, as RQ8 predicted.

Summary and parameters

print(wb_growth.summary().tables[1])

The summary table gives the parameter estimates, standard errors, z-scores, p-values, and 95% confidence intervals.

Random-effects covariances (wb_growth.cov_re) report intercept and slope variations.

Did the workshop effect last?

wb_workshop = smf.mixedlm(
    "wellbeing ~ C(semester) * workshop",
    panel,
    groups=panel["student_id"],
    re_formula="~time"
).fit()

C(semester) gives each semester its own average; the interaction lets the gap differ by semester.

Every record is used, including the first year of students who later left.

The effect fades

Predictions for the average student

pred_grid = pd.MultiIndex.from_product(
    [[1, 2, 3, 4], ["Invited", "Not invited"]],
    names=["semester", "workshop"]
).to_frame().reset_index(drop=True)
pred_grid["time"] = pred_grid["semester"] - 1
pred_grid["predicted_wellbeing"] = wb_workshop.predict(pred_grid)

predict() evaluates the population-level fixed effects.

Every gap after semester 1 differs from the baseline gap. A yearly follow-up session might keep the benefit going.

A singular fit / boundary convergence

gpa_growth = smf.mixedlm("gpa ~ time", panel_gpa, groups=panel_gpa["student_id"], re_formula="~time").fit()
# ConvergenceWarning: The MLE may be on the boundary of the parameter space

The model is more complex than the data supports: GPA slopes hardly vary across students (variance ≈ 0).

Simplify to a random intercept only.

Grades stay stable

gpa_growth = smf.mixedlm("gpa ~ time", panel_gpa, groups=panel_gpa["student_id"]).fit()

The slope of time is -0.000: almost exactly zero (\(p = 0.95\)).

Wellbeing falls while grades hold. A non-change is a finding too.

Students within supervisors

wb_levels = smf.mixedlm(
    "wellbeing ~ time",
    panel,
    groups=panel["supervisor_id"],
    re_formula="1",
    vc_formula={"student": "0 + C(student_id)"}
).fit()

vc_formula specifies nested student random intercepts. Who supervises a student makes a measurable difference (13% of variance).

Other outcomes: counts and binary

Supervisor meetings per semester are a count outcome.

A count model for meetings

meetings_gee = smf.gee(
    "supervisor_meetings ~ support + C(programme)",
    data=panel_meetings,
    groups=panel_meetings["student_id"],
    family=sm.families.Poisson()
).fit()
np.exp(meetings_gee.params)

Each extra point of support: about 49% more meetings. PhD students: about 23% more.

The direction is uncertain: more meetings may also make students feel supported.

Mixed models and missing data

Students who left after year 1 were not a random group (Chapter 6).

Mixed models use every record each student gave, instead of dropping them.

The results hold as long as leaving depends on things the model can see: missing at random.

Reporting a mixed model

• the fixed effects, with confidence intervals;
• the random effects, as SDs or variance shares;
• the number of observations and of groups;
• which random effects were included.

Writing it up

A linear mixed-effects model with random intercepts and slopes for students (2326 observations of 600 students) showed that wellbeing declined by 0.43 points per semester (95% CI [0.22, 0.65]). Students differed considerably in their overall wellbeing (SD of intercepts = 10.7); 78% of the variation was between students.

In your field: agriculture

chick_model = smf.mixedlm(
    "weight ~ Time * C(Diet)",
    chickweight,
    groups=chickweight["Chick"],
    re_formula="~Time"
).fit()

50 chicks weighed from birth to 21 days on four diets: repeated measurements in groups.

Diet 1: about 6.3 g a day. Diet 3 adds 5.1 g a day, 82% faster.

Practical lab: the Chapter 10 playground

Work through the playground exercises in your browser, with hints and solutions.

statsmodels runs right in your browser through Pyodide; every exercise also runs locally in VS Code or Jupyter.

Practical exercises 1–3: growth models

1. A random intercept model for sleep over time: ICC and change.
2. Add a random slope and compare with likelihood ratio test.
3. Does part-time students’ wellbeing change at a different rate?

Practical exercises 4–6: levels and simulation

4. How much of the variation in study hours lies between supervisors?
5. With supervisors who differ more, what happens to the false alarm rate?
6. Explain to a fellow student why a regression of all semester records is wrong.

Try this yourself

Look at your own data for structure.

• is anyone measured more than once?
• do participants share a group: class, clinic, supervisor, site?
• estimate the ICC with a null model;
• calculate the design effect;
• decide which random effects your model needs.

Troubleshooting guide (Part 1)

Troubleshooting guide (Part 2)

Completion checklist

• the grouping structure is identified;
• the ICC and design effect are reported;
• random intercepts are included for each grouping;
• random slopes are tested with likelihood ratio tests;
• singular fits / boundary variances are simplified;
• fixed effects come with confidence intervals;
• numbers of observations and groups are reported.

Misconceptions to leave behind (Part 1)

Misconceptions to leave behind (Part 2)

The chapter in one sentence

When observations share a person or a group, model that structure: mixed models separate differences between groups from change within them and get the uncertainty right.

Next: Chapter 11

The next chapter turns from explanation to prediction:

• machine learning and generalisation;
• overfitting;
• a fair test with training and test data;
• cross-validation and tuning;
• the final test, and predictions about people.

Questions

What is the grouping structure in your own data?

How many truly independent units does your study have?