Playground: Chapter 10

Mixed-Effects Models

This page practises the ideas of Chapter 10 in four parts. Practise the chapter works through the book’s exercises, with starter code, hints, and solutions. Go further goes beyond the book, with new questions and real longitudinal datasets. Check your understanding asks short questions whose answers open with a click. Do it yourself gives open tasks with an empty code space and no answers: you choose and check the models yourself, as you will in your own thesis.

The exercises use statsmodels.formula.api.mixedlm. The panel DataFrame panel is prepared as in the chapter: one row per student per semester, with the background variables and time counted from 0 (semester 1) to 3 (semester 4). In the exercises, replace each ______ with your own code and press Run Code.

Practise the chapter

These are the exercises at the end of Chapter 10, with the same numbers.

Exercise 1: Sleep over the semesters

Fit a random intercept model for sleep_hours over time. Calculate the ICC, and describe whether sleep changes over the semesters.

Solution

sleep_ri = smf.mixedlm("sleep_hours ~ time", panel, groups=panel["student_id"]).fit()
var_between = sleep_ri.cov_re.iloc[0, 0]
var_within = sleep_ri.scale
icc = var_between / (var_between + var_within)
print(f"ICC: {icc:.3f}")
print("95% CI for slope:\n", sleep_ri.conf_int().loc["time"])

The ICC is 0.83: 83% of the variation in sleep is between students, so students differ a great deal in how much they sleep, but each sleeps about the same from one semester to the next. The slope of time is practically zero (0.0005 hours per semester), and its 95% confidence interval ([-0.015, 0.016]) is narrow and centred on zero: sleep does not change over the two years, and the interval rules out any change larger than about a minute per semester.

Exercise 2: Is a random slope worth keeping?

Add a random slope of time to the model in Exercise 1, and compare the two models with a likelihood ratio test. Decide whether the random slope is worth keeping.

Solution

sleep_rs = smf.mixedlm("sleep_hours ~ time", panel, groups=panel["student_id"], re_formula="~time").fit()
lr_stat = -2 * (sleep_ri.llf - sleep_rs.llf)
p_val = stats.chi2.sf(lr_stat, df=2)
print(f"LR Chi2(2): {lr_stat:.2f}, p-value: {p_val:.4f}")
print("Random effects cov:\n", sleep_rs.cov_re)

The likelihood ratio test favours the random slope (\(\chi^2 \approx 7.8\), 2 df, \(p \approx 0.02\)). But sleep_rs.cov_re shows that the slopes hardly vary: their standard deviation is \(\sqrt{0.003} \approx 0.05\) hours (3 minutes) per semester, against \(\sqrt{0.90} \approx 0.95\) hours between students’ overall levels. Statistically detectable, practically unimportant. Either model is defensible; for a thesis, the simpler random intercept model with a sentence noting the test would be a reasonable choice.

Exercise 3: Do part-time students change differently?

Fit a model of wellbeing over time with study mode and their interaction, with random intercepts and slopes, and describe whether part-time students’ wellbeing changes at a different rate.

Solution

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

Full-time students start at 62.8 and decline by 0.47 points per semester. Part-time students start 3.9 points lower (95% CI [-5.9, -1.9]), but their slope differs by only 0.14 points per semester, with an interval from -0.34 to 0.61 that includes zero. Part-time students’ wellbeing is lower from the start, and the gap stays about the same over the two years: there is no evidence that it changes at a different rate.

Exercise 4: Supervisors and study hours

Fit a model with nested random intercepts for supervisors and students for study_hours, and calculate how much of its variation lies between supervisors.

Solution

# The fitted model directly estimates the supervisor variance component:
var_sup = study_levels.cov_re.iloc[0, 0]
var_stu = study_levels.vcomp[0]
var_res = study_levels.scale
var_tot = var_sup + var_stu + var_res

The share between supervisors is 0 (the supervisor variance is estimated as 0.000): the model found no differences between supervisors at all, beyond what differences between their students explain. About 91% of the variation lies between students and 9% within students from semester to semester. Supervisors matter for wellbeing, as the chapter showed, but not for how much students study. The zero variance is the model’s way of saying that the supervisor level can be dropped.

Exercise 5: Supervisors who differ more

Change the simulation at the start of the chapter so that supervisors differ more (style_sd = 10), and describe what happens to the false alarm rate of the ordinary regression.

Solution

runs = [simulate_clustered(style_sd=10) for _ in range(100)]

The ordinary regression now declares the useless course significant in about 46% of the studies, against about 27% in the chapter; the mixed model stays near 5–8%. When supervisors differ more, students of the same supervisor resemble each other more (the ICC rises), so the 200 students carry even less independent information: the study is really about 20 supervisors. The ordinary regression still counts 200 independent students, so its standard error is too small by more, and the difference between two sets of 10 supervisors looks like an effect of the course nearly half the time.

Exercise 6: Explaining it to a fellow student

In your own words, explain to a fellow student why an ordinary regression of all the semester records would be wrong. Write your answer first, then open the model answer.

Model answer

“Each student appears up to four times, and a student’s four records are much more alike than records of four different students: someone who sleeps little in semester 1 usually sleeps little in semester 4 too. An ordinary regression treats the 2,326 records as 2,326 independent pieces of evidence, but they are really about 600 people, and together they carry only about as much information as 700 or so independent records. So the ordinary regression thinks it knows more than it does: its standard errors are too small, and its p-values too small, and it will find ‘effects’ that are just differences between a few students. A mixed model knows which records belong to the same student and gets the uncertainty right.”

Go further

These exercises go beyond the book.

Exercise 7: How much information is in the records?

The chapter calculated the design effect for wellbeing. Calculate it for GPA and for sleep, and the number of independent records each set of records is worth.

Solution

design_effect = 1 + (records_per_student - 1) * icc_gpa

With 3.9 records per student, GPA has an ICC of 0.71 and a design effect of 3.0: its 2,325 records are worth about 767 independent observations. Sleep has an ICC of 0.83 and a design effect of 3.4, so its records are worth only about 679. The more stable a variable is within each student, the more its repeated records repeat each other, and the less each new semester adds.

Exercise 8: Averaging first, or a mixed model

One simple way to avoid repeated records is to average each student’s records and run an ordinary regression on the 600 averages. Compare this with a mixed model, and with the wrong analysis of all records, for the difference in wellbeing between part-time and full-time students.

Solution

# The three models illustrate the point:
print("Averages OLS:\n", m_avg.summary().tables[1])
print("Mixed Model:\n", m_mix.summary().tables[1])
print("Naive OLS:\n", m_naive.summary().tables[1])

The averages and the mixed model agree almost perfectly: part-time students are ~3.7 points lower, with a standard error of ~0.98 in both. The ordinary regression of all 2,326 records gives a similar estimate (-3.8) but an artificially small standard error (0.54), because it counts every repeated record as independent. For a variable that does not change within a student, such as study mode, averaging is a sound shortcut. However, it cannot answer questions about change over time, and it weights a student with one record equally with one with four; that is why mixed models are essential.

Exercise 9: How long a study needs to run

The sleepstudy dataset follows 18 subjects over 10 days of sleep restriction. Fit the random slope model to the first five days only (days 0 to 4) and compare it with the model for all ten days.

Solution

early_data = sleepstudy.query("Days <= 4")
early = smf.mixedlm("Reaction ~ Days", early_data, groups=early_data["Subject"], re_formula="~Days").fit()

With five days, the estimated average slope is 8.2 ms per day (SE 2.4), compared to 10.5 ms per day (SE 1.5) with all ten days: a milder, less precise estimate from half the observation window. The residual variance is also smaller early on, showing that subjects become more erratic as cumulative sleep debt builds up.

Exercise 10: Growth curves for chicks

The chickweight dataset records the weight of 50 chicks measured over 21 days. Fit a model where each chick has its own initial weight and growth rate over Time.

Solution

chick_model = smf.mixedlm("weight ~ Time", chickweight, groups=chickweight["Chick"], re_formula="~Time").fit()
print(chick_model.summary())

The average chick gains approximately 8.5 grams per day. The standard deviation of chick growth rates is about 3.8 grams per day: chicks vary substantially in growth trajectory.

Check your understanding

Answer each question in your own words first, then click to see a model answer.

1. Why are repeated records of the same student not independent?

Answer

Because each student brings stable characteristics to every record: their habits, circumstances, personality, and the way they answer questions. Two records of the same student are therefore more alike than records of two different students, and the second record partly repeats the information in the first.

2. What does the intraclass correlation (ICC) measure?

Answer

The share of all the variation that lies between groups, such as students or supervisors, rather than within them. An ICC of 0 means that the grouping does not matter; an ICC near 1 means that members of a group are almost identical. The higher the ICC, the more wrong an analysis that ignores the grouping would be.

3. What is the difference between a fixed effect and a random effect?

Answer

A fixed effect is a quantity of direct interest, estimated as one number: the average effect of time, or the difference between study modes. A random effect describes how groups vary around the average, estimated as a standard deviation: how much students’ starting levels or growth rates differ. The individual students are not of interest in themselves; they are a sample from a population whose variation the random effect describes.

4. What does a singular fit / zero variance estimate mean, and what should you do about it?

Answer

That the model’s random effects are more complex than the data can support: typically, a variance component estimated as zero or on the boundary. It does not mean that the model is wrong in its fixed effects. The usual response is to simplify the random effects, for example by dropping a random slope or a grouping level that shows zero variation, and to report the simplification transparently.

5. In the chapter’s simulation, why did the ordinary regression produce so many false alarms?

Answer

The course was given to supervisors, and there were only 20 of them. Students of the same supervisor shared the supervisor’s style, so the 200 students were far from 200 independent pieces of evidence. The ordinary regression counted them as independent, underestimated the standard error, and mistook chance differences between two small groups of supervisors for an effect of the course.

Do it yourself

These tasks have no starter code and no answers. Each code space below is empty and ready to run: write your own code, as you would for a thesis. The datasets (students, semesters, panel, sleepstudy, chickweight) and standard packages (numpy, pandas, statsmodels) are already imported.

1. Model the change in exercise_days over the four semesters. Choose the random effects with a likelihood ratio test, calculate the ICC, and write two sentences reporting the result.

2. Test whether the effect of the workshop invitation on wellbeing over time differs between Master’s and PhD students, with a three-way interaction of time, workshop, and programme. Describe the result with predicted wellbeing for the four groups in semester 2.

3. Write the paragraph reporting a mixed model for GPA over the semesters, like the chapter’s reporting paragraph for wellbeing, including the random effects and the number of observations and students, with every number calculated by code.

4. Explore chick growth curves with dietary interactions using chickweight. Fit a model with Time * C(Diet) and random intercepts and slopes per chick, and interpret whether diets alter growth rates.

Work on your own computer

Download the Chapter 10 project

The project contains the data and all four parts of this page as a clean Python project: the exercises with blanks, the questions to check your understanding (answers in solutions.py), and the open tasks, each with space to write your code.

Download chapter10.zip, unzip it, and open the folder in VS Code or your favourite Python editor.
Install the required packages:

pip install pandas numpy scipy statsmodels matplotlib seaborn