Linear Mixed Models Method Card Template (Statistics, R and Python)

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

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.

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.

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.

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

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$):

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:

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.

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).

# 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

(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.

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}")

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.

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.

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")

(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

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])

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.

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

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].

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

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