# Linear Mixed Models: [Your study question or course] > A method card for linear mixed models (multilevel or hierarchical models) for clustered and longitudinal data: > students within tutorial groups, repeated measurements within people. It covers why ordinary least squares (OLS) > standard errors are wrong for such data, random intercepts and random slopes, fixed versus random effects, the > intraclass correlation (ICC), REML versus ML and when to use each, the same model fitted with `lme4::lmer` in R and > `MixedLM` in statsmodels, the degrees-of-freedom debate, convergence and singular fits, centring predictors, and > reporting. It is for students meeting mixed models in a course, thesis writers with multi-wave or nested data, and > analysts who need correct uncertainty for grouped data. It speeds up learning with one fixed modelling sequence, > a full worked example, a half-solved example, an exercise and a self-test with a separate answer key. The worked > examples use fictional data; replace them with your own. > Source basis: standard mixed-model methods (Laird and Ware 1982; Patterson and Thompson 1971; Bates et al. 2015; > Satterthwaite 1946; Kenward and Roger 1997) and the NIST/SEMATECH e-Handbook, written in original words. ## Document control | Field | Value | |---|---| | Course or project | [Course code and name, or thesis chapter] | | Learner or analyst | [Your name] | | Data set | [Name, source, number of clusters, observations per cluster] | | Software and versions | R 4.5.1, lme4 2.1.0 (optional: lmerTest 3.2.1, pbkrtest 0.5.5); Python 3.14, numpy 2.5.4, pandas 3.1.0, scipy 1.18.1, statsmodels 0.15.0 | | Estimation | [REML for final model; ML for comparing fixed effects] | | Date | YYYY-MM-DD | | Version | [0.1] | ## When to use it **One-sentence rule:** use a linear mixed model when a continuous outcome is measured on units that are grouped (clusters) or measured repeatedly (subjects), so that observations in the same group are more alike than observations in different groups, and you want correct standard errors and, optionally, group-specific intercepts or slopes. Think in **levels**. Level 1 is the observation (one quiz score, one student). Level 2 is the group it belongs to (the student whose scores they are, the tutorial group the student attends). A predictor lives at one level: week of term varies within a student (level 1); teaching format is the same for everyone in a tutorial group (level 2). The level of a predictor decides how much information you have about it, and this is exactly what OLS gets wrong. | Use when | Do not use when | |---|---| | Repeated measurements per subject, possibly with missing time points | Each unit is measured once and units are independent (use OLS) | | Units nested in clusters (students in classes, loaves in batches) | Very few clusters (fewer than about 5) for estimating a variance; use fixed cluster effects or cluster-robust SEs | | Time points unequal across subjects, or time is continuous | The outcome is binary or a count (use a generalised linear mixed model, `glmer`) | | You want individual trajectories (random slopes) | You care only about population-average effects and the clusters are a nuisance (GEE is an alternative) | | Data are unbalanced and some measurements are missing at random | Missingness depends on the unobserved outcome itself (MNAR); sensitivity analysis needed | **Alternatives:** repeated-measures ANOVA (complete balanced data only, see stats-repeated-measures); OLS with cluster-robust standard errors (many clusters, interest in averages only); generalised estimating equations (GEE); fixed-effects (within) regression with a dummy per cluster (controls all cluster-level confounding but cannot estimate cluster-level predictors); aggregating to cluster means (simple and valid for a cluster-level predictor in balanced data, see Example B). ## The question and the data **Research question (Example A):** how fast do quiz scores improve over the first weeks in the Lakemoor University statistics learning lab, and do students differ in their starting level and their rate of improvement? Eight students at Lakemoor University (fictional) each took a short practice quiz (0 to 50) in weeks 0, 1, 2, 3 and 4 of the learning lab: 40 observations nested in 8 students. | Variable | Level | Role | Values | |---|---|---|---| | id | 2 (student) | Grouping factor | S1 to S8 | | week | 1 (observation) | Predictor, time since start | 0 to 4 | | score | 1 (observation) | Outcome | 0 to 50 | | id | week 0 | week 1 | week 2 | week 3 | week 4 | |---|---|---|---|---|---| | *Example:* S1 | 29 | 29 | 31 | 35 | 39 | | S2 | 26 | 31 | 34 | 38 | 41 | | S3 | 26 | 29 | 35 | 35 | 38 | | S4 | 36 | 38 | 39 | 45 | 46 | | S5 | 33 | 39 | 41 | 44 | 44 | | S6 | 32 | 35 | 36 | 38 | 39 | | S7 | 30 | 27 | 28 | 33 | 37 | | S8 | 27 | 27 | 26 | 30 | 29 | The table is shown wide for reading; the code stores it long (40 rows of id, week, score), which is the layout every mixed-model function needs. ```vega-lite { "$schema": "https://vega.github.io/schema/vega-lite/v5.json", "description": "Quiz score by week for eight Lakemoor students, one line per student", "width": 360, "height": 240, "data": {"values": [ {"id": "S1", "week": 0, "score": 29}, {"id": "S1", "week": 1, "score": 29}, {"id": "S1", "week": 2, "score": 31}, {"id": "S1", "week": 3, "score": 35}, {"id": "S1", "week": 4, "score": 39}, {"id": "S2", "week": 0, "score": 26}, {"id": "S2", "week": 1, "score": 31}, {"id": "S2", "week": 2, "score": 34}, {"id": "S2", "week": 3, "score": 38}, {"id": "S2", "week": 4, "score": 41}, {"id": "S3", "week": 0, "score": 26}, {"id": "S3", "week": 1, "score": 29}, {"id": "S3", "week": 2, "score": 35}, {"id": "S3", "week": 3, "score": 35}, {"id": "S3", "week": 4, "score": 38}, {"id": "S4", "week": 0, "score": 36}, {"id": "S4", "week": 1, "score": 38}, {"id": "S4", "week": 2, "score": 39}, {"id": "S4", "week": 3, "score": 45}, {"id": "S4", "week": 4, "score": 46}, {"id": "S5", "week": 0, "score": 33}, {"id": "S5", "week": 1, "score": 39}, {"id": "S5", "week": 2, "score": 41}, {"id": "S5", "week": 3, "score": 44}, {"id": "S5", "week": 4, "score": 44}, {"id": "S6", "week": 0, "score": 32}, {"id": "S6", "week": 1, "score": 35}, {"id": "S6", "week": 2, "score": 36}, {"id": "S6", "week": 3, "score": 38}, {"id": "S6", "week": 4, "score": 39}, {"id": "S7", "week": 0, "score": 30}, {"id": "S7", "week": 1, "score": 27}, {"id": "S7", "week": 2, "score": 28}, {"id": "S7", "week": 3, "score": 33}, {"id": "S7", "week": 4, "score": 37}, {"id": "S8", "week": 0, "score": 27}, {"id": "S8", "week": 1, "score": 27}, {"id": "S8", "week": 2, "score": 26}, {"id": "S8", "week": 3, "score": 30}, {"id": "S8", "week": 4, "score": 29} ]}, "mark": {"type": "line", "point": true}, "encoding": { "x": {"field": "week", "type": "quantitative", "title": "Week in learning lab", "axis": {"tickCount": 5}}, "y": {"field": "score", "type": "quantitative", "title": "Quiz score (0 to 50)", "scale": {"domain": [20, 50]}}, "color": {"field": "id", "type": "nominal", "title": "Student"} } } ``` Read the chart. Students start at different levels (S4 at 36, S2 and S3 at 26), which suggests a random intercept. Their lines also climb at different rates (S2 gains 15 points, S8 only 2), which suggests a random slope. A single OLS line through all 40 points would ignore both facts. ## Assumptions and how to check each | Assumption | How to check | R | Python | If violated | |---|---|---|---|---| | Correct grouping structure (nesting, crossing) | Draw the hierarchy; check IDs are unique within the right level | `table(d$id)` | `d.groupby("id").size()` | Re-specify grouping, e.g. `(1 \| school/class)` | | Linearity of fixed effects | Plot outcome against predictors; residuals against fitted | `plot(m2)` | `plt.scatter(fit.fittedvalues, fit.resid)` | Add polynomial or spline terms | | Level-1 residuals roughly normal, constant variance | Q-Q plot and residual plot | `qqnorm(resid(m2))` | `stats.probplot(fit.resid, plot=None)` | Transform; robust SEs; check outliers | | Random effects roughly normal | Q-Q plot of predicted random effects (BLUPs) | `qqnorm(ranef(m2)$id[, 1])` | `fit.random_effects` | Usually robust for fixed effects; report | | Random effects independent of predictors | Compare within and between effects of level-1 predictors | Add cluster mean of predictor | Same | Include cluster-mean (contextual) term or use within-centring | | Missing data at most MAR | Describe dropout; compare completers | Design | Design | Sensitivity analysis | | Enough clusters for the variance you estimate | Count clusters | `ngrps(m2)` | `fit.model.n_groups` | Simplify random effects; fixed cluster effects | ## How it works **Plain language.** OLS assumes every observation is independent. Five quiz scores from the same student are not: if S4 is strong in week 0, S4 is likely to be strong in week 3 as well. Treating 40 correlated scores as 40 independent pieces of information gives the wrong standard errors. For predictors that are constant within a group, OLS acts as if it had more information than it has, so its SEs are too small and p-values too optimistic. For predictors that vary within groups, OLS lumps the stable between-student differences into the error term, so its SEs can instead be too large. A mixed model fixes both problems by modelling where the variation comes from. A mixed model has two parts. The **fixed effects** describe the average relationship in the population: the typical starting score and the typical weekly gain. The **random effects** describe how individual groups deviate from that average: each student gets their own intercept (random intercept) and, if you include it, their own slope (random slope). The model does not estimate a separate free parameter for each student; it estimates the **variance** of those deviations. Each student's own line is then a compromise between their own data and the average line, pulled toward the average more strongly when the student has few or noisy observations. This "shrinkage" (partial pooling) is a key benefit. Two-level growth model (observation $t$ within student $i$): $$ \begin{aligned} y_{ti} &= \beta_{0i} + \beta_{1i}\,\text{week}_{ti} + e_{ti}, & e_{ti} &\sim N(0, \sigma^2) \\ \beta_{0i} &= \gamma_{00} + u_{0i}, \quad \beta_{1i} = \gamma_{10} + u_{1i}, & \begin{pmatrix} u_{0i} \\ u_{1i} \end{pmatrix} &\sim N\!\left(\mathbf{0}, \begin{pmatrix} \tau_{00} & \tau_{01} \\ \tau_{01} & \tau_{11} \end{pmatrix}\right) \end{aligned} $$ In matrix form for all data, $\mathbf{y} = X\boldsymbol{\beta} + Z\mathbf{u} + \mathbf{e}$ (Laird and Ware 1982). **Intraclass correlation.** From the random-intercept-only (null) model, the share of total variance that lies between groups, which is also the expected correlation between two observations from the same group: $$ \text{ICC} = \frac{\tau_{00}}{\tau_{00} + \sigma^2} $$ **Design effect.** With $m$ observations per cluster, the variance of a cluster-level comparison is inflated by about $\text{DEFF} = 1 + (m - 1)\,\text{ICC}$ relative to independent data. With ICC = 0.47 and $m = 5$, DEFF = 2.9: 40 correlated scores carry the information of about 14 independent ones for between-student questions. **REML versus ML.** Maximum likelihood (ML) estimates variances as if the fixed effects were known, so it underestimates them in small samples (like dividing by $n$ instead of $n - 1$). Restricted maximum likelihood (REML; Patterson and Thompson 1971) corrects this and is the default for final estimates. Rules: (1) report REML estimates; (2) to compare models that differ in **fixed effects**, refit with ML, because REML likelihoods of models with different fixed effects are not comparable; (3) to compare models that differ only in **random effects**, REML likelihoods may be compared. | Symbol | Meaning | Example A (model m2, REML) | |---|---|---| | $\gamma_{00}$ | Average score at week 0 (fixed intercept) | 29.60 | | $\gamma_{10}$ | Average weekly gain (fixed slope) | 2.39 | | $\tau_{00}$ | Variance of student intercepts | 14.03 (SD 3.75) | | $\tau_{11}$ | Variance of student slopes | 0.56 (SD 0.75) | | $\tau_{01}$ | Intercept-slope covariance | 0.76 (correlation 0.27) | | $\sigma^2$ | Residual variance within student | 2.70 (SD 1.64) | | ICC | Between-student share, null model | 0.473 | ```mermaid flowchart TD A["Draw the hierarchy: what is nested in what?"] --> B["Null model: random intercept only"] B --> C["ICC: is there clustering worth modelling?"] C --> D["Add level-1 and level-2 fixed effects"] D --> E["Add random slopes justified by design"] E --> F["Converged and not singular?"] F -->|No| G["Simplify random part, rescale or centre, change optimiser"] G --> E F -->|Yes| H["Compare random parts with REML LRT"] H --> I["Compare fixed parts with ML LRT or Wald tests"] I --> J["Refit final model with REML"] J --> K["Check residuals and random effects"] K --> L["Report fixed effects with CIs, variance components, ICC"] ``` ## Worked example A (full) **Questions.** (1) What is the average weekly gain? $H_0: \gamma_{10} = 0$. (2) Do students differ in their weekly gain? $H_0: \tau_{11} = 0$. α = 0.05. **R code (lme4).** ```r # Lakemoor University learning lab: weekly quiz score (0-50), 8 students, weeks 0-4 d <- data.frame( id = factor(rep(paste0("S", 1:8), each = 5)), week = rep(0:4, times = 8), score = c(29, 29, 31, 35, 39, 26, 31, 34, 38, 41, 26, 29, 35, 35, 38, 36, 38, 39, 45, 46, 33, 39, 41, 44, 44, 32, 35, 36, 38, 39, 30, 27, 28, 33, 37, 27, 27, 26, 30, 29) ) library(lme4) # install.packages("lme4") once # 1. Null model and ICC m0 <- lmer(score ~ 1 + (1 | id), data = d) vc <- as.data.frame(VarCorr(m0)) icc <- vc$vcov[1] / sum(vc$vcov) print(VarCorr(m0)); cat("ICC =", round(icc, 3), "\n") # 2. Random intercept, then random intercept and slope (REML) m1 <- lmer(score ~ week + (1 | id), data = d) m2 <- lmer(score ~ week + (week | id), data = d) summary(m2) confint(m2, parm = "beta_", method = "Wald") # 3. Compare random-effect structures (same fixed effects): REML fits, LRT anova(m1, m2, refit = FALSE) # 4. Compare fixed effects: must use ML (anova() refits with ML by default) anova(m0, m1) # 5. Naive OLS for contrast summary(lm(score ~ week, data = d))$coefficients ``` ```text (output abridged) Groups Name Std.Dev. id (Intercept) 4.0283 Residual 4.2529 ICC = 0.473 Linear mixed model fit by REML ['lmerMod'] Formula: score ~ week + (week | id) REML criterion at convergence: 185 Random effects: Groups Name Variance Std.Dev. Corr id (Intercept) 14.0293 3.7456 week 0.5602 0.7485 0.27 Residual 2.6958 1.6419 Number of obs: 40, groups: id, 8 Fixed effects: Estimate Std. Error t value (Intercept) 29.6000 1.3985 21.165 week 2.3875 0.3221 7.413 2.5 % 97.5 % (Intercept) 26.858955 32.341045 week 1.756273 3.018727 m1: score ~ week + (1 | id) m2: score ~ week + (week | id) npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq) m1 4 198.76 205.52 -95.381 190.76 m2 6 197.01 207.15 -92.507 185.01 5.7477 2 0.05648 . refitting model(s) with ML (instead of REML) m0: score ~ 1 + (1 | id) m1: score ~ week + (1 | id) npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq) m0 3 247.87 252.94 -120.937 241.87 m1 4 200.26 207.01 -96.128 192.26 49.616 1 1.87e-12 *** Estimate Std. Error t value Pr(>|t|) (Intercept) 29.6000 1.270124 23.304808 4.149942e-24 week 2.3875 0.518526 4.604398 4.530660e-05 ``` **Python code (statsmodels MixedLM) on the same data.** ```python import numpy as np import pandas as pd import statsmodels.formula.api as smf scores = [29, 29, 31, 35, 39, 26, 31, 34, 38, 41, 26, 29, 35, 35, 38, 36, 38, 39, 45, 46, 33, 39, 41, 44, 44, 32, 35, 36, 38, 39, 30, 27, 28, 33, 37, 27, 27, 26, 30, 29] d = pd.DataFrame({"id": np.repeat([f"S{i}" for i in range(1, 9)], 5), "week": np.tile(range(5), 8), "score": scores}) # 1. Null model and ICC m0 = smf.mixedlm("score ~ 1", d, groups="id").fit(reml=True) tau2, sigma2 = float(m0.cov_re.iloc[0, 0]), m0.scale print(f"ICC = {tau2 / (tau2 + sigma2):.3f} (tau2 = {tau2:.3f}, sigma2 = {sigma2:.3f})") # 2. Random intercept and slope (REML) m2 = smf.mixedlm("score ~ week", d, groups="id", re_formula="~week").fit(reml=True) print(m2.summary()) # 3. ML fits for comparing fixed effects (likelihood-ratio test by hand) from scipy import stats m0ml = smf.mixedlm("score ~ 1", d, groups="id").fit(reml=False) m1ml = smf.mixedlm("score ~ week", d, groups="id").fit(reml=False) lr = 2 * (m1ml.llf - m0ml.llf) print(f"ML logLik: m0 = {m0ml.llf:.3f}, m1 = {m1ml.llf:.3f}; LR = {lr:.3f}, p = {stats.chi2.sf(lr, 1):.2e}") ``` ```text ICC = 0.473 (tau2 = 16.229, sigma2 = 18.087) Mixed Linear Model Regression Results ======================================================== Model: MixedLM Dependent Variable: score No. Observations: 40 Method: REML No. Groups: 8 Scale: 2.6958 Min. group size: 5 Log-Likelihood: -92.5074 Max. group size: 5 Converged: Yes Mean group size: 5.0 -------------------------------------------------------- Coef. Std.Err. z P>|z| [0.025 0.975] -------------------------------------------------------- Intercept 29.600 1.398 21.166 0.000 26.859 32.341 week 2.387 0.322 7.413 0.000 1.756 3.019 id Var 14.028 5.789 id x week Cov 0.759 0.832 week Var 0.560 0.307 ======================================================== ML logLik: m0 = -120.937, m1 = -96.128; LR = 49.616, p = 1.87e-12 ``` **Where R and Python agree and differ.** | Quantity | lme4 | statsmodels | Comment | |---|---|---|---| | Fixed intercept, slope | 29.6000, 2.3875 | 29.600, 2.387 | Identical (display rounding) | | Fixed-effect SEs | 1.3985, 0.3221 | 1.398, 0.322 | Identical | | Intercept variance | 14.0293 | 14.028 | Tiny optimiser difference | | Intercept-slope link | correlation 0.27 | covariance 0.759 | Same thing: $0.27 \times 3.7456 \times 0.7485 \approx 0.76$ | | Residual variance | 2.6958 | Scale 2.6958 | Identical | | REML log-likelihood | -92.507 | -92.5074 | Identical | | p-values for fixed effects | None printed | z test, normal reference | The degrees-of-freedom debate, below | | SEs for variances | Not printed (use profile CIs) | Printed | Wald SEs for variances are unreliable; do not test variances with them | **Calculation by hand for the key statistics.** - ICC from the null model: $\tau_{00} = 4.0283^2 = 16.23$, $\sigma^2 = 4.2529^2 = 18.09$, ICC $= 16.23 / (16.23 + 18.09) = 0.473$. Nearly half the variation in scores is between students. - Fixed slope Wald test: $t = 2.3875 / 0.3221 = 7.41$; 95% Wald CI $2.3875 \pm 1.96 \times 0.3221 = [1.76, 3.02]$. - Fixed-effect LRT (ML): $2 \times (-96.128 - (-120.937)) = 49.62$ on 1 df, *p* < .001. - Random-slope LRT (REML): $185.01$ subtracted from $190.76$ gives 5.75 on 2 df, *p* = .056. Because a variance cannot be negative, the null value sits on the boundary and the plain chi-square p-value is conservative; a common correction (Stram and Lee 1994) uses a 50:50 mixture of chi-square with 1 and 2 df, which gives *p* of about .036 here. Either way the evidence that students differ in their gain is modest; keep the random slope because the design (growth over time) makes it plausible and Barr et al. (2013) argue for keeping slopes the design justifies. **OLS versus mixed model.** OLS gives the same fixed estimates here (balanced data) but different SEs: intercept SE 1.27 (mixed 1.40, so OLS is too small for the between-student quantity) and slope SE 0.52 (mixed 0.32, so OLS is too large for this within-student quantity, because it leaves the stable student differences in the error). OLS SEs are wrong in both directions, depending on the level of the predictor. **Centring.** Week is coded 0 to 4, so the intercept (29.6) is the expected score at the first quiz, which is meaningful. If week were coded as calendar weeks 10 to 14, the intercept would be the predicted score at week 0 of the calendar, far outside the data, and its variance and its correlation with the slope would be hard to interpret. Centring week at its mean (week minus 2) would move the intercept to $29.6 + 2 \times 2.3875 = 34.4$, the expected score at mid-period. Choose the zero point that answers your question, and say what it is. **Interpretation.** On average, students scored 29.6 at the first quiz and gained 2.39 points per week (95% CI [1.76, 3.02]). Students differed substantially in starting level (SD of intercepts 3.75 points) and somewhat in weekly gain (SD of slopes 0.75 points: about two thirds of students gain between 1.6 and 3.1 points per week if the slope distribution is normal). Students who started higher tended to gain slightly faster (correlation 0.27), but with 8 students this correlation is very imprecise. **APA-style report.** A linear mixed model with random intercepts and random slopes for week by student (REML, lme4) showed that quiz scores increased by 2.39 points per week, *b* = 2.39, *SE* = 0.32, 95% CI [1.76, 3.02], from an average of 29.60 at week 0. The intraclass correlation from the null model was .47. Students varied in their starting level (*SD* = 3.75) and in their weekly gain (*SD* = 0.75); residual *SD* = 1.64. Adding week improved fit, χ²(1) = 49.62, *p* < .001 (ML). ## Worked example B (half solved) **Setting.** Lakemoor University ran STAT 101 in six tutorial groups of five students. Three groups (T1, T3, T5) used a flipped format, three (T2, T4, T6) a standard lecture format. Outcome: exam score (0 to 100). Format is a level-2 (tutorial-level) predictor. | Tutorial | Format | Exam scores | Mean | |---|---|---|---| | T1 | flipped | 71, 67, 64, 66, 63 | 66.2 | | T2 | lecture | 73, 62, 62, 68, 67 | 66.4 | | T3 | flipped | 64, 69, 74, 61, 67 | 67.0 | | T4 | lecture | 48, 44, 54, 48, 57 | 50.2 | | T5 | flipped | 69, 55, 67, 69, 70 | 66.0 | | T6 | lecture | 54, 51, 52, 64, 53 | 54.8 | ```r library(lme4) b <- data.frame( tutorial = factor(rep(paste0("T", 1:6), each = 5)), format = factor(rep(c("flipped", "lecture"), times = 3, each = 5), levels = c("lecture", "flipped")), exam = c(71, 67, 64, 66, 63, 73, 62, 62, 68, 67, 64, 69, 74, 61, 67, 48, 44, 54, 48, 57, 69, 55, 67, 69, 70, 54, 51, 52, 64, 53) ) summary(lm(exam ~ format, data = b))$coefficients # naive OLS mb <- lmer(exam ~ format + (1 | tutorial), data = b) summary(mb) confint(mb, parm = "beta_", method = "Wald") ``` ```text (output abridged) Estimate Std. Error t value Pr(>|t|) (Intercept) 57.133333 1.757524 32.50786 8.816877e-24 formatflipped 9.266667 2.485514 3.72827 8.663933e-04 Random effects: Groups Name Variance Std.Dev. tutorial (Intercept) 30.01 5.478 Residual 24.90 4.990 Number of obs: 30, groups: tutorial, 6 Fixed effects: Estimate Std. Error t value (Intercept) 57.133 3.415 16.730 formatflipped 9.267 4.830 1.919 2.5 % 97.5 % (Intercept) 50.4400568 63.82661 formatflipped -0.1990557 18.73239 ``` ```python import numpy as np import pandas as pd import statsmodels.formula.api as smf b = pd.DataFrame({ "tutorial": np.repeat([f"T{i}" for i in range(1, 7)], 5), "flipped": np.repeat([1, 0] * 3, 5), "exam": [71, 67, 64, 66, 63, 73, 62, 62, 68, 67, 64, 69, 74, 61, 67, 48, 44, 54, 48, 57, 69, 55, 67, 69, 70, 54, 51, 52, 64, 53], }) ols = smf.ols("exam ~ flipped", data=b).fit() print("OLS:", ols.params.round(3).to_dict(), "SE:", ols.bse.round(3).to_dict()) mb = smf.mixedlm("exam ~ flipped", b, groups="tutorial").fit(reml=True) print(mb.summary().tables[1]) ``` ```text OLS: {'Intercept': 57.133, 'flipped': 9.267} SE: {'Intercept': 1.758, 'flipped': 2.486} Coef. Std.Err. z P>|z| [0.025 0.975] Intercept 57.133 3.415 16.730 0.000 50.440 63.827 flipped 9.267 4.830 1.919 0.055 -0.199 18.732 tutorial Var 30.007 5.355 ``` 1. OLS says the flipped format raises scores by 9.27 points with SE = [ ] and *p* = [ ]. 2. The mixed model gives the same estimate with SE = [ ]. The SE grew by a factor of [ ]. 3. Conditional ICC (with format in the model) = 30.01 / (30.01 + [ ]) = [ ]. 4. Format varies only between tutorials. How many independent units inform the format comparison: 30 or [ ]? If you compare the six tutorial means with a pooled t-test (three per format), you get *t*(4) = 1.92, *p* = .127. Compare this with the statsmodels z-test p-value of [ ]. Which p-value is more trustworthy with six clusters, and why? [ ] 5. Report sentence: [ ]. ## Exercise C (unsolved) Harbour City Library (fictional) recorded the number of reference questions answered on five weekdays at each of four branches. | Branch | Mon | Tue | Wed | Thu | Fri | |---|---|---|---|---|---| | North | 20 | 24 | 24 | 18 | 19 | | South | 18 | 22 | 20 | 22 | 14 | | East | 25 | 20 | 22 | 20 | 19 | | West | 21 | 22 | 19 | 20 | 22 | Tasks: (a) fit a random-intercept model `q ~ 1 + (1 | branch)` in R and the equivalent `MixedLM` in Python; (b) report the branch variance, the residual variance and the ICC; (c) explain the message both programs print; (d) decide what to do and how to report it; (e) say why four branches is a weak basis for estimating a between-branch variance in any case. (Counts are treated as continuous here for practice; for real count data see stats-count-regression.) ## Reading software output | Output line | What it means | What to report | |---|---|---| | `fit by REML` / `Method: REML` | Estimation method | State it; ML for fixed-effect comparisons | | `REML criterion` / `Log-Likelihood` | Fit measure; minus two times the log-likelihood | Only for model comparison | | Random effects: Variance, Std.Dev. | Variance components for each random term | Variances or SDs for every component | | `Corr` (lme4) / `Cov` (statsmodels) | Association between random intercepts and slopes | Correlation, with caution in small samples | | Residual / Scale | Level-1 residual variance | Residual SD or variance | | `Number of obs`, `groups` | Observations and clusters | Both, always | | Fixed effects Estimate, Std. Error | Population-average coefficients | Estimate with 95% CI | | `t value` (lme4) | Wald ratio; no p-value by design | With p from Satterthwaite or Kenward-Roger df | | `P>\|z\|` (statsmodels) | Wald z-test, normal reference, infinite df | Treat as approximate; optimistic with few clusters | | `boundary (singular) fit` / `SingularMatrixWarning` | A variance estimated at zero or a correlation at plus or minus 1 | Simplify the random part; report it | | `Model failed to converge` / `Converged: No` | Optimiser did not settle | Rescale, centre, simplify, change optimiser | | `anova(m0, m1)` Chisq | Likelihood-ratio test between nested models | χ²(df), p and which likelihood (ML or REML) | ## Reporting template **Fill-in sentence.** A linear mixed model with [fixed effects] and random [intercepts / intercepts and slopes for predictor] by [grouping factor] was fitted with [REML] in [lme4 2.1.0]. [Outcome] [increased / decreased] by [b] per [unit] (*SE* = [value], 95% CI [[lower], [upper]], *t*([Satterthwaite df]) = [value], *p* = [value]). The intraclass correlation from the null model was [ICC]. Random intercepts had *SD* = [value] [and random slopes *SD* = [value], correlation = [value]]; residual *SD* = [value]. Data: [n observations] in [k clusters]. | Fixed effect | Estimate | SE | 95% CI | df method | p | |---|---|---|---|---|---| | *Example:* Intercept (week 0) | 29.60 | 1.40 | [26.86, 32.34] | Wald | [value] | | *Example:* Week | 2.39 | 0.32 | [1.76, 3.02] | Wald | [value] | | Random effect | Variance | SD | Correlation | |---|---|---|---| | *Example:* Student intercept | 14.03 | 3.75 | | | *Example:* Student slope (week) | 0.56 | 0.75 | 0.27 | | *Example:* Residual | 2.70 | 1.64 | | **Results checklist** - [ ] Hierarchy described: observations, clusters, number of each, observations per cluster. - [ ] Full model formula, including the random part and the coding or centring of predictors. - [ ] Estimation method (REML or ML) and software with version. - [ ] Fixed effects with SE, 95% CI and the df method behind any p-value. - [ ] Variance components (and correlations) and ICC. - [ ] Model comparisons with the correct likelihood (ML for fixed effects). - [ ] Convergence and singular-fit status; any simplification of the random part and why. - [ ] Residual and random-effect diagnostics summarised. - [ ] Missing data: how much, and the MAR assumption. ## Common mistakes and misconceptions 1. **Fitting OLS to clustered data.** The estimates may be fine but the SEs are wrong: too small for cluster-level predictors (Example B: SE 2.49 instead of 4.83), possibly too large for within-cluster predictors. 2. **Comparing models with different fixed effects using REML likelihoods.** Refit with ML (R's `anova()` does it automatically; in statsmodels pass `reml=False`). 3. **"Random effects are for variables I don't care about."** Random versus fixed is about whether the levels are a sample from a larger population whose variance you want to estimate, not about interest. 4. **A grouping factor with very few levels as a random effect.** With four branches (Exercise C) the variance is barely estimable. Use fixed effects for the groups, or report the limitation. 5. **Ignoring a singular-fit message.** A zero variance or a correlation of plus or minus 1 means the random part is too complex for the data. Simplify and say so; do not report the correlation as a finding. 6. **Trusting z-test p-values with few clusters.** With six clusters a z-test treats the df as infinite. Use Satterthwaite or Kenward-Roger df (Example B: *p* = .055 with z, .127 with 4 df). 7. **Not centring or rescaling.** Predictors on large scales (year 2026, income in USD) make intercepts meaningless and can cause convergence warnings. Centre at a meaningful value; rescale to sensible units. 8. **Reading the ICC from a model with predictors as "the" ICC.** The null-model ICC describes clustering in the raw outcome; the conditional ICC describes what is left after the fixed effects. Say which one you report. 9. **Calling a random slope "non-significant, therefore students do not differ".** Variance tests have low power and a boundary problem; absence of evidence is not evidence of absence. ## The degrees-of-freedom debate, convergence and singular fits **Why lme4 prints no p-values.** In a mixed model the reference distribution of the Wald t ratio is not exactly a t distribution with a known df, because the SE depends on estimated variance components. The lme4 authors (Bates et al. 2015) leave p-values out on purpose. Options, from most to least common: - **Satterthwaite approximation** (Satterthwaite 1946), via the `lmerTest` package: `library(lmerTest)` before fitting makes `summary()` show df and p. It is the usual default in psychology and education. - **Kenward-Roger** (Kenward and Roger 1997) adjusts both the SE and the df; available through `pbkrtest` and `emmeans`. It is often preferred with few clusters. - **Likelihood-ratio tests** with ML fits (`anova(m0, m1)`), which work well with many clusters. - **Parametric bootstrap** (`pbkrtest`, or simulate and refit), robust but slow. - **Wald z** (statsmodels default): acceptable with many clusters (roughly 40 or more), optimistic with few. In balanced designs the answer for a cluster-level predictor is often clear: Example B has 6 clusters and 2 fixed parameters at that level, so the natural df is 6 minus 2 = 4, which matches the pooled t-test on cluster means. Refitting Example B after `library(lmerTest)` (version 3.2.1) confirms it: Satterthwaite df = 4.00, *t* = 1.92, *p* = .127, and `summary(mb, ddf = "Kenward-Roger")` gives the same df and p here. Whatever you choose, name the method. **Convergence warnings.** Messages such as "Model failed to converge" mean the optimiser stopped near, but not exactly at, the optimum. Work through: (1) check the scale of predictors and rescale (divide by 10 or 100); (2) centre predictors; (3) try another optimiser (`lmerControl(optimizer = "bobyqa")`) or `allFit()` to see whether estimates agree across optimisers; (4) simplify the random part. If all optimisers give the same estimates to reported precision, the warning was a false alarm; say so. **Singular fit.** A singular (boundary) fit means some variance is estimated as zero or some correlation as plus or minus 1. It is common with few clusters or with complex random-slope structures. The fit is not wrong as such, but the random part is more complex than the data can support. Simplify (drop the correlation with `(1 + week || id)`, then the slope), report the simplification, and do not interpret a correlation of 1 as a real finding. **Fixed versus random, once more.** Treat a grouping factor as random when its levels are a sample from a wider set (students, classes, batches, branches), there are enough of them (often 5 to 10 at the very least), and you want to generalise beyond them. Treat it as fixed when the levels are the complete set of interest (three named ovens) or there are very few. ## Self-test questions 1. Explain in two sentences why OLS standard errors are wrong for clustered data, and why they can be wrong in either direction. 2. A null model gives a between-cluster variance of 9 and a residual variance of 27. Compute the ICC and the design effect for clusters of 11. 3. You want to test whether adding a fixed effect of gender improves the model. Should the two models be fitted with REML or ML? Why? 4. What does `boundary (singular) fit` tell you, and name two remedies. 5. Why does statsmodels give *p* = .055 for the flipped format while a cluster-level t-test gives *p* = .127? 6. Your time variable is the calendar year (2020 to 2026). What problem does this create for the random intercept and its correlation with the slope, and how do you fix it? ## My notes and next steps | Item | My entry | |---|---| | Confidence (1 to 5) for drawing the hierarchy and choosing the grouping factor | [ ] | | Confidence (1 to 5) for choosing random intercepts and slopes | [ ] | | Confidence (1 to 5) for REML vs ML and model comparison | [ ] | | Confidence (1 to 5) for writing the report | [ ] | | Open questions | [What I still do not understand] | | Mistakes logged in stats-error-log | [Link or entry numbers] | | Next review date (spacing: 2 days, 1 week, 3 weeks) | YYYY-MM-DD | | Next card to study | [stats-logistic-regression / stats-regression-model-report] | ## Sources and further reading Based on standard mixed-model methods described in the sources below; all wording, data and examples on this card are original. - Patterson, H. D., and Thompson, R. (1971). Recovery of inter-block information when block sizes are unequal. *Biometrika*, 58(3), 545–554. - Laird, N. M., and Ware, J. H. (1982). Random-effects models for longitudinal data. *Biometrics*, 38(4), 963–974. - Satterthwaite, F. E. (1946). An approximate distribution of estimates of variance components. *Biometrics Bulletin*, 2(6), 110–114. - Kenward, M. G., and Roger, J. H. (1997). Small sample inference for fixed effects from restricted maximum likelihood. *Biometrics*, 53(3), 983–997. - Bates, D., Mächler, M., Bolker, B., and Walker, S. (2015). Fitting linear mixed-effects models using lme4. *Journal of Statistical Software*, 67(1), 1–48. - Kuznetsova, A., Brockhoff, P. B., and Christensen, R. H. B. (2017). lmerTest package: tests in linear mixed effects models. *Journal of Statistical Software*, 82(13), 1–26. - Stram, D. O., and Lee, J. W. (1994). Variance components testing in the longitudinal mixed effects model. *Biometrics*, 50(4), 1171–1177. - Barr, D. J., Levy, R., Scheepers, C., and Tily, H. J. (2013). Random effects structure for confirmatory hypothesis testing: keep it maximal. *Journal of Memory and Language*, 68(3), 255–278. - Seabold, S., and Perktold, J. (2010). statsmodels: econometric and statistical modeling with Python. *Proceedings of the 9th Python in Science Conference*, 92–96. - NIST/SEMATECH e-Handbook of Statistical Methods, section 7.4.4 What are variance components?, https://www.itl.nist.gov/div898/handbook/ (accessed 2026-10-25). Public domain. - Further reading (by name only): de Smith, M. J. (2024). *Statistical Analysis Handbook*, 2024 edition. StatsRef, https://www.statsref.com/; OpenIntro Statistics; APA *Publication Manual* (7th ed.) and APA JARS-Quant for reporting standards. ## Answer key **Example B.** 1. OLS: SE = 2.49, *p* < .001 (8.66e-04). 2. Mixed model: SE = 4.83; the SE grew by a factor of about 1.94 (4.83 / 2.49). 3. Conditional ICC = 30.01 / (30.01 + 24.90) = 0.55. 4. Six independent units (tutorials), not 30. The statsmodels z-test gives *p* = .055; the cluster-level *t*(4) gives *p* = .127. The t-based value is more trustworthy: with six clusters the variance between tutorials is estimated from very little information, and a z-test ignores that uncertainty (it uses infinite df). Satterthwaite and Kenward-Roger df (lmerTest) both give 4.00 here, *p* = .127. 5. "A linear mixed model with a random intercept for tutorial group (REML) estimated that the flipped format was associated with exam scores 9.27 points higher than the lecture format, *SE* = 4.83, 95% CI [-0.20, 18.73], *t*(4) = 1.92, *p* = .127. Tutorial groups differed considerably (intercept *SD* = 5.48; residual *SD* = 4.99; conditional ICC = .55). With only six groups the comparison is imprecise." Note also that format was not randomised to tutorials here, so a causal claim is not justified. **Exercise C.** R prints `boundary (singular) fit: see help('isSingular')`; the branch variance is 0.000 and the residual variance 6.261 (*SD* 2.50), so the ICC is 0. The fixed intercept is 20.55 (*SE* = 0.56), identical to the OLS mean and SE. Python prints `SingularMatrixWarning` and `ConvergenceWarning: The MLE may be on the boundary of the parameter space` with `branch Var 0.000` and Scale 6.2605. The branch means (North 21.0, South 19.2, East 21.2, West 20.8) vary no more than expected from day-to-day noise alone, so the estimate of between-branch variance is zero. What to do: report that the between-branch variance was estimated at zero, that the model reduces to OLS, and either keep the random intercept (the estimates are the same) or drop it and say why. Four clusters give a variance estimate based on effectively three df, so even a non-zero value would be very imprecise; with the branches being the complete set of interest, fixed branch effects would be a reasonable alternative. **Self-test.** 1. Observations in the same cluster are correlated, so they carry less independent information than OLS assumes. For cluster-level predictors OLS overstates the information (SEs too small); for within-cluster predictors it leaves stable cluster differences in the error term (SEs can be too large). 2. ICC = 9 / (9 + 27) = 0.25. DEFF = 1 + (11 - 1) × 0.25 = 3.5. 3. ML. REML maximises a likelihood of residuals that depends on the fixed-effects design, so REML likelihoods for models with different fixed effects are not comparable. 4. A variance component was estimated at zero (or a correlation at plus or minus 1): the random part is too complex for the data. Remedies: drop the random-effect correlation, drop the random slope, use fixed effects for few groups, or collect more clusters. 5. statsmodels uses a z-test (infinite df). The cluster-level test uses 4 df, reflecting that only six tutorials inform the comparison; the t distribution with 4 df has much heavier tails. 6. The intercept refers to year 0, far outside the data, so the intercept variance and its correlation with the slope describe an extrapolation and are often near plus or minus 1; the model may also fail to converge. Centre time, e.g. year minus 2020 (intercept at the first year) or minus 2023 (the middle). ## Revision history | Version | Date | Author | Change | |---|---|---|---| | 0.1 | YYYY-MM-DD | [Name] | First version from template | | [0.2] | YYYY-MM-DD | [Name] | [Replaced example data with own data] |