The V Lab
AudienceLearners in medicine, public health, sociology, psychology, and health data science
Study timeApproximately 180–240 minutes
PrerequisitesLinear and logistic regression, confidence intervals, variance, and basic R

About the data in this tutorial Every patient, clinic, hospital, student, school, therapist, and follow-up record was simulated with a fixed random seed and contains no real personal information. The simulation mechanisms validate code and interpretation; they provide no evidence about real treatments, institutions, or social policies.

How to use this tutorial

A multilevel analysis should not begin with a complicated lmer() formula. Work through the material in this order:

Identify analysis units and dependence → state the target relationship → decompose variance → separate within- and between-cluster effects → choose a random structure → fit and diagnose → predict, conduct sensitivity analyses, and report

Four reproducible simulations run through this page: a continuous outcome for patients nested in clinics, a random-slope model for students nested in schools, a three-level model with therapy sessions nested in people nested in therapists, and a binary outcome for patients nested in hospitals. All primary code relies only on base R, lme4, and knitr for rendering the page.

Learning objectives

After completing this tutorial, you should be able to:

  • identify observation, individual, and cluster levels, including nesting, cross-classification, and repeated measures;
  • explain fixed effects, random effects, variance components, and partial pooling without interpreting “random” as unimportant;
  • calculate and correctly interpret an ICC from a random-intercept model;
  • distinguish complete pooling, no pooling, and partial pooling;
  • use cluster-mean decomposition to separate within- and between-cluster relationships;
  • interpret random slopes, intercept–slope covariance, and cross-level interactions;
  • choose an LMM or GLMM for continuous, binary, and longitudinal outcomes;
  • distinguish cluster-specific conditional effects from population-marginal effects;
  • use ML and REML correctly and recognize convergence, singular fits, and influential clusters;
  • explain the limits of mixed models for missingness, serial correlation, and causal identification;
  • use lme4::lmer(), lme4::glmer(), lme4::VarCorr(), and prediction interfaces; and
  • write an auditable report covering sample levels, fixed effects, random structure, diagnostics, and limitations.

1 Why ordinary regression may not be enough

1.1 Independence follows from the data-generating process

Patients in the same clinic may share referral criteria, clinicians, equipment, and regional environments. Students in the same school share curricula and resources. Diary observations from the same person share stable traits and histories. Even when every row looks like a complete record, these shared causes generally make errors dependent.

Ignoring that dependence can:

  • distort standard errors and confidence intervals;
  • mix relationships from different levels;
  • compress cluster heterogeneity into a single residual term;
  • conflate prediction for a new person with prediction for an observed person; and
  • create a large row count without much independent information for higher-level variables.
design_map <- data.frame(
  Field = c(
    "Medicine: one outcome", "Medicine: repeated follow-up",
    "Sociology", "Psychology", "Medicine: binary outcome"
  ),
  `Level 1` = c("Patient", "Follow-up occasion", "Student", "Therapy session", "Patient"),
  `Level 2` = c("Clinic", "Patient", "School", "Person", "Hospital"),
  `Level 3 or crossed structure` = c(
    "—", "Hospital", "Community may cross schools", "Therapist", "—"
  ),
  `Primary question` = c(
    "Are patients in the same clinic correlated?",
    "Do trajectories vary by patient and hospital?",
    "Do individual SES and school-composition effects differ?",
    "Do daily and chronic stress have different relationships?",
    "What are hospital heterogeneity and the conditional OR?"
  ),
  check.names = FALSE
)

knitr::kable(design_map, caption = "Levels, units, and questions in three application fields")
Levels, units, and questions in three application fields
Field Level 1 Level 2 Level 3 or crossed structure Primary question
Medicine: one outcome Patient Clinic — Are patients in the same clinic correlated?
Medicine: repeated follow-up Follow-up occasion Patient Hospital Do trajectories vary by patient and hospital?
Sociology Student School Community may cross schools Do individual SES and school-composition effects differ?
Psychology Therapy session Person Therapist Do daily and chronic stress have different relationships?
Medicine: binary outcome Patient Hospital — What are hospital heterogeneity and the conditional OR?

An ID does not automatically require a random intercept A “level” is defined relative to the research question and dependence mechanism; it is not an intrinsic property of a variable. Explain what records share and whether inference concerns people, clusters, or both. Do not add (1 | hospital_id) mechanically just because a data set contains hospital_id.

1.2 Nesting, crossing, and multiple membership

  • Nesting: each patient belongs to one clinic, and each occasion belongs to one person.
  • Cross-classification: a student belongs to both a school and a residential community when schools do not belong to only one community.
  • Multiple membership: a patient is served jointly by several therapists during treatment.
  • Unequal cluster sizes: hospitals contain different numbers of patients, as is common in real data.
  • Informative cluster size: cluster size itself is associated with outcome risk and may require additional handling.
Check your understanding: If 5,000 patients come from eight hospitals, what is the hospital-level sample size? For a hospital-level exposure or policy effect, independent information comes primarily from eight hospitals, not 5,000 patients. More patients improve estimation of hospital means but do not create more independent hospitals.

2 The two-level random-intercept model

2.1 Common relationships and cluster deviations

Let individual ii be nested in cluster jj:

Yij=β0+β1Xij+β2Zj+u0j+εij, Y_{ij}=\beta_0+\beta_1X_{ij}+\beta_2Z_j+u_{0j}+\varepsilon_{ij}, u0j∼N(0,τ00),εij∼N(0,σ2). u_{0j}\sim N(0,\tau_{00}),\qquad \varepsilon_{ij}\sim N(0,\sigma^2).

XijX_{ij} is an individual-level variable and ZjZ_j a cluster-level variable. Fixed effects β\beta describe average conditional relationships shared across clusters. The random intercept u0ju_{0j} describes cluster jj’s deviation from the overall intercept, with distributional variance τ00\tau_{00}.

“Fixed” and “random” are technical terms A fixed effect does not mean that a variable never changes, and a random effect does not mean unimportant or necessarily sampled by simple random sampling. The model estimates common coefficients β\beta directly while using a distribution to describe many cluster deviations uju_j.

2.2 The ICC describes correlation, not a causal share

The intraclass correlation coefficient in an empty random-intercept model is:

ICC=ρ=τ00τ00+σ2. ICC=\rho=\frac{\tau_{00}}{\tau_{00}+\sigma^2}.

It is the model correlation between two randomly selected individuals from the same cluster. In an empty model, it can also describe the share of total variance at the between-cluster level. It cannot be interpreted as “the cluster caused this percentage of the outcome.” A conditional ICC after covariate adjustment answers a different variance-decomposition question.

A small ICC does not automatically justify ignoring clustering. Its consequences also depend on average cluster size, number of clusters, the level at which an exposure varies, and cluster-size imbalance.

2.3 Complete pooling, no pooling, and partial pooling

Strategy Approach Benefit Limitation
Complete pooling Ignore cluster differences Simple Compresses dependence and heterogeneity into residuals
No pooling Estimate every cluster separately Retains cluster differences Small clusters are unstable and many parameters are required
Partial pooling Shrink cluster estimates toward the overall distribution Stable while retaining heterogeneity Relies on the random-effect distribution and model

Under a simple empty model, a cluster deviation can be approximated by:

ûj≈λj(Y‾j−β̂0),λj=τ00τ00+σ2/nj. \widehat u_j\approx \lambda_j(\bar Y_j-\widehat\beta_0),\qquad \lambda_j=\frac{\tau_{00}}{\tau_{00}+\sigma^2/n_j}.

Small clusters contain less information, have smaller λj\lambda_j, and shrink more. Large clusters retain more of their own information. Partial pooling does not “average away institutional differences”; it makes a principled compromise between cluster data and the population distribution.

3 Medical application I: patients nested in clinics

3.1 Simulated cohort and time zero

Forty-eight clinics each enroll 22–36 patients. The intervention is assigned at the clinic level, the outcome is six-month systolic blood pressure, and pretreatment age and baseline blood pressure are recorded. Because treatment varies at the clinic level, its higher-level independent information comes from the number of clinics.

set.seed(20260821)
J <- 48L
n_j <- sample(22:36, J, replace = TRUE)
clinic <- sprintf("C%02d", seq_len(J))
treat_j <- sample(rep(0:1, each = J / 2))
u0 <- rnorm(J, 0, 4.5)

med <- do.call(rbind, lapply(seq_len(J), function(j) {
  n <- n_j[j]
  age <- pmin(pmax(rnorm(n, 60, 11), 30), 85)
  baseline <- rnorm(n, 138 + 0.08 * (age - 60), 12)
  sbp_6m <- 132 - 4.2 * treat_j[j] +
    0.55 * (baseline - 138) + 0.12 * (age - 60) +
    u0[j] + rnorm(n, 0, 9)
  data.frame(
    clinic = clinic[j],
    treatment = treat_j[j],
    age_c10 = (age - 60) / 10,
    baseline_c10 = (baseline - 138) / 10,
    sbp_6m = sbp_6m
  )
}))
med$treatment <- factor(
  med$treatment,
  0:1,
  c("Usual care", "Intervention")
)

med_audit <- data.frame(
  Patients = nrow(med),
  Clinics = length(unique(med$clinic)),
  `Smallest clinic` = min(table(med$clinic)),
  `Median clinic` = median(table(med$clinic)),
  `Largest clinic` = max(table(med$clinic)),
  check.names = FALSE
)

knitr::kable(med_audit, caption = "Audit of the simulated patient–clinic cohort")
Audit of the simulated patient–clinic cohort
Patients Clinics Smallest clinic Median clinic Largest clinic
1425 48 22 29.5 36

3.2 Empty model, ICC, and partial pooling

m_med_empty <- lme4::lmer(
  sbp_6m ~ 1 + (1 | clinic),
  data = med,
  REML = TRUE
)

med_empty_icc <- icc_random_intercept(m_med_empty, "clinic")
med_empty_vc <- as.data.frame(lme4::VarCorr(m_med_empty))

med_empty_variance <- data.frame(
  Component = c("Clinic random intercept", "Patient-level residual"),
  Variance = c(
    med_empty_vc$vcov[
      med_empty_vc$grp == "clinic" & is.na(med_empty_vc$var2)
    ],
    stats::sigma(m_med_empty)^2
  ),
  `Standard deviation` = sqrt(c(
    med_empty_vc$vcov[
      med_empty_vc$grp == "clinic" & is.na(med_empty_vc$var2)
    ],
    stats::sigma(m_med_empty)^2
  )),
  check.names = FALSE
)

knitr::kable(
  med_empty_variance,
  digits = 3,
  caption = "Variance components from the empty medical random-intercept model"
)
Variance components from the empty medical random-intercept model
Component Variance Standard deviation
Clinic random intercept 27.1 5.21
Patient-level residual 128.6 11.34

The empty-model ICC is 0.174: two randomly selected patients from the same clinic have a model-scale correlation of about 17.4% in six-month systolic blood pressure. This does not mean that clinics causally produced 17.4% of blood pressure.

clinic_raw <- aggregate(sbp_6m ~ clinic, data = med, FUN = mean)
clinic_raw$n <- as.numeric(table(med$clinic)[clinic_raw$clinic])
clinic_re <- lme4::ranef(m_med_empty)$clinic
clinic_partial <- data.frame(
  clinic = rownames(clinic_re),
  partial = unname(lme4::fixef(m_med_empty)[1] + clinic_re[, 1])
)
clinic_compare <- merge(clinic_raw, clinic_partial, by = "clinic")
clinic_compare <- clinic_compare[order(clinic_compare$sbp_6m), ]
plot_y <- seq_len(nrow(clinic_compare))

plot(
  clinic_compare$sbp_6m,
  plot_y,
  pch = 16,
  col = palette_ml["orange"],
  xlim = range(c(clinic_compare$sbp_6m, clinic_compare$partial)),
  yaxt = "n",
  xlab = "Mean six-month systolic blood pressure",
  ylab = "Clinics ordered by raw mean",
  main = "Raw means and partial pooling",
  las = 1
)
segments(
  clinic_compare$sbp_6m,
  plot_y,
  clinic_compare$partial,
  plot_y,
  col = palette_ml["gray"]
)
points(
  clinic_compare$partial,
  plot_y,
  pch = 17,
  col = palette_ml["teal"]
)
abline(v = lme4::fixef(m_med_empty)[1], lty = 2, col = palette_ml["navy"])
legend(
  "bottomright",
  legend = c("Raw clinic mean", "Partially pooled mean", "Overall intercept"),
  col = c(palette_ml["orange"], palette_ml["teal"], palette_ml["navy"]),
  pch = c(16, 17, NA),
  lty = c(NA, NA, 2),
  bty = "n"
)
Forty-eight clinics are ordered by raw mean systolic blood pressure. Orange points are raw means, teal triangles are partially pooled estimates, and most connecting segments point toward the overall mean. Smaller clinics generally move farther.

Raw clinic means and partially pooled means from the empty random-intercept model. Segments connect the two estimates for each clinic; smaller clinics generally shrink more toward the overall mean.

Conditional modes are often called BLUPs, but they are not error-free “clinic truths.” Any use for institutional ranking, sanctions, or resource allocation must display uncertainty, case mix, and the extent of shrinkage.

3.3 Adjusted clinic random-intercept model

m_med <- lme4::lmer(
  sbp_6m ~ treatment + baseline_c10 + age_c10 + (1 | clinic),
  data = med,
  REML = TRUE
)

med_labels <- c(
  `(Intercept)` = "Usual care with baseline covariates at their centers",
  treatmentIntervention = "Intervention vs usual care",
  baseline_c10 = "Per 10 mmHg higher baseline systolic pressure",
  age_c10 = "Per 10 years older"
)
med_fixed <- lmm_fixed_table(m_med, med_labels)
med_adjusted_icc <- icc_random_intercept(m_med, "clinic")

knitr::kable(
  med_fixed,
  digits = 3,
  caption = "Fixed effects from the adjusted patient–clinic linear mixed model"
)
Fixed effects from the adjusted patient–clinic linear mixed model
Term Estimate Standard error 95% CI lower 95% CI upper t value
(Intercept) Usual care with baseline covariates at their centers 132.219 1.031 130.199 134.24 128.27
treatmentIntervention Intervention vs usual care -4.113 1.457 -6.969 -1.26 -2.82
baseline_c10 Per 10 mmHg higher baseline systolic pressure 5.568 0.217 5.143 5.99 25.71
age_c10 Per 10 years older 0.573 0.228 0.127 1.02 2.52

Mean six-month systolic blood pressure is 4.11 mmHg lower under the intervention than usual care (Wald 95% CI -6.97 to -1.26). This is a conditional mean difference given baseline pressure, age, and the clinic random intercept. Because treatment is assigned at the clinic level, causal interpretation additionally requires clinic-level exchangeability, adequate overlap, and a clear assignment mechanism.

The adjusted ICC is 0.207, higher than the empty-model value of 0.174. This is not contradictory: removing explainable patient-level variance can change the composition of the remaining total variance. The ICC is not a “model quality score” that must decline whenever covariates are added.

med_residuals <- residuals(m_med, type = "pearson")
med_fitted <- fitted(m_med)
old_par <- par(mfrow = c(1, 2), mar = c(4.2, 4.2, 3, 1))

plot(
  med_fitted,
  med_residuals,
  pch = 16,
  cex = 0.55,
  col = grDevices::adjustcolor(palette_ml["teal"], alpha.f = 0.45),
  xlab = "Conditional fitted value",
  ylab = "Pearson residual",
  main = "Residuals versus fitted values",
  las = 1
)
abline(h = 0, lty = 2, col = palette_ml["vermillion"])

stats::qqnorm(
  as.numeric(scale(med_residuals)),
  pch = 16,
  cex = 0.55,
  col = palette_ml["blue"],
  main = "Residual Q–Q plot"
)
stats::qqline(as.numeric(scale(med_residuals)), col = palette_ml["vermillion"])
Two diagnostic panels. Residual points in the left panel scatter around a horizontal zero line. Standardized residual points in the right panel follow the normal reference line approximately with random tail departures.

Conditional residual diagnostics for the patient–clinic linear mixed model. The left panel compares fitted values with residuals, and the right panel is a normal Q–Q plot of standardized residuals. Diagnostics look for nonlinearity, heteroskedasticity, and tail departures. 

par(old_par)

Residual plots are not the only diagnostics. Also inspect clinic-level influence, clinic size, the random-effect distribution, the functional form of baseline blood pressure, and possible clinic-specific residual variances.

4 Within- and between-cluster effects

4.1 Why centering changes the research question

For an individual-level variable XijX_{ij} measured within clusters:

Xij=(Xij−X‾j)+X‾j, X_{ij}=(X_{ij}-\bar X_j)+\bar X_j, Yij=β0+βW(Xij−X‾j)+βBX‾j+u0j+εij. Y_{ij}=\beta_0+\beta_W(X_{ij}-\bar X_j)+ \beta_B\bar X_j+u_{0j}+\varepsilon_{ij}.

  • βW\beta_W: the relationship associated with a one-unit difference between two individuals in the same cluster;
  • βB\beta_B: the relationship associated with a one-unit difference between cluster means; and
  • βB−βW\beta_B-\beta_W: often called the contextual contrast, although causal interpretation still requires additional assumptions.

Grand-mean centering primarily changes the interpretation of zero and may improve numerical stability; it does not automatically separate within- and between-cluster relationships. Group-mean centering helps decompose the question, but it does not automatically remove all cluster-level confounding.

One coefficient can mix two questions that point in different directions Entering time-varying stress without decomposition mixes “this person is more stressed today than usual” with “this person is chronically more stressed than other people.” Likewise, a student’s household SES and the socioeconomic composition of the school should not be summarized by one undecomposed coefficient.

5 Sociology application: students nested in schools

5.1 Simulating SES, school support, and achievement

set.seed(20260822)
J <- 70L
n_j <- sample(24:36, J, replace = TRUE)
school <- sprintf("S%02d", seq_len(J))
school_ses <- rnorm(J, 0, 0.8)
school_ses_gc <- school_ses - mean(school_ses)
support_gc <- rnorm(J, 0, 1)
support_gc <- support_gc - mean(support_gc)
rho <- 0.25
z0 <- rnorm(J)
z1 <- rnorm(J)
b0 <- 5.0 * z0
b1 <- 1.5 * (rho * z0 + sqrt(1 - rho^2) * z1)

soc <- do.call(rbind, lapply(seq_len(J), function(j) {
  n <- n_j[j]
  ses_wc <- rnorm(n, 0, 0.9)
  ses_wc <- ses_wc - mean(ses_wc)
  ses <- school_ses[j] + ses_wc
  achievement <- 70 +
    2.8 * ses_wc + 5.0 * school_ses_gc[j] +
    2.0 * support_gc[j] - 1.2 * ses_wc * support_gc[j] +
    b0[j] + b1[j] * ses_wc + rnorm(n, 0, 7)
  data.frame(
    school = school[j],
    ses = ses,
    school_support_gc = support_gc[j],
    achievement = achievement
  )
}))
soc$school_ses <- ave(soc$ses, soc$school, FUN = mean)
soc$ses_wc <- soc$ses - soc$school_ses
soc$school_ses_gc <- soc$school_ses - mean(soc$school_ses)

soc_audit <- data.frame(
  Students = nrow(soc),
  Schools = length(unique(soc$school)),
  `Smallest school` = min(table(soc$school)),
  `Largest school` = max(table(soc$school)),
  `Largest absolute within-school mean of ses_wc` =
    max(abs(tapply(soc$ses_wc, soc$school, mean))),
  check.names = FALSE
)

knitr::kable(
  soc_audit,
  digits = 5,
  caption = "Audit of the student–school simulation and centering"
)
Audit of the student–school simulation and centering
Students Schools Smallest school Largest school Largest absolute within-school mean of ses_wc
2098 70 24 36 0

Within every school, ses_wc has mean zero and represents a student’s position relative to peers in that school. Here, school_ses_gc locates a school’s mean SES relative to an overall mean weighted by the number of students. If the research question requires every school to receive equal weight, first extract one unique mean per school before centering and document that choice in the analysis plan. School support is grand-mean centered at the school level with schools weighted equally.

school_means <- aggregate(
  cbind(achievement, school_ses) ~ school,
  data = soc,
  FUN = mean
)
old_par <- par(mfrow = c(1, 2), mar = c(4.2, 4.2, 3, 1))
plot(
  soc$ses_wc,
  soc$achievement,
  pch = 16,
  cex = 0.45,
  col = grDevices::adjustcolor(palette_ml["teal"], alpha.f = 0.30),
  xlab = "Student SES minus school mean SES",
  ylab = "Achievement",
  main = "Within-school relationship",
  las = 1
)
abline(stats::lm(achievement ~ ses_wc, data = soc), col = palette_ml["navy"], lwd = 2)
plot(
  school_means$school_ses,
  school_means$achievement,
  pch = 17,
  col = palette_ml["orange"],
  xlab = "School mean SES",
  ylab = "School mean achievement",
  main = "Between-school relationship",
  las = 1
)
abline(
  stats::lm(achievement ~ school_ses, data = school_means),
  col = palette_ml["vermillion"],
  lwd = 2
)
Two scatterplots. More than two thousand student points in the left panel show a positive within-school association between SES and achievement. Seventy school-mean points in the right panel show a steeper positive between-school association.

Within- and between-school relationships for student SES. The left panel relates a student’s SES relative to the school mean to achievement; the right panel relates school mean SES to school mean achievement. These are questions at different levels. 

par(old_par)

5.2 Random slopes and a cross-level interaction

A random-intercept model assumes that the within-school SES slope is identical in every school. Allowing that slope to vary gives:

Yij=β0+β1Xij+u0j+u1jXij+εij, Y_{ij}=\beta_0+\beta_1X_{ij}+u_{0j}+u_{1j}X_{ij}+\varepsilon_{ij}, (u0ju1j)∼N[(00),(τ00τ01τ01τ11)]. \begin{pmatrix}u_{0j}\\u_{1j}\end{pmatrix} \sim N\left[ \begin{pmatrix}0\\0\end{pmatrix}, \begin{pmatrix}\tau_{00}&\tau_{01}\\\tau_{01}&\tau_{11}\end{pmatrix} \right].

Here β1\beta_1 is the average within-school slope, τ11\tau_{11} describes between-school slope heterogeneity, and τ01\tau_{01} describes how intercepts and slopes vary together. After adding a cross-level interaction with school support ZjZ_j, the average SES slope is β1+β3Zj\beta_1+\beta_3Z_j.

m_soc_empty <- lme4::lmer(
  achievement ~ 1 + (1 | school),
  data = soc,
  REML = TRUE
)

m_soc <- lme4::lmer(
  achievement ~ ses_wc + school_ses_gc + school_support_gc +
    ses_wc:school_support_gc + (1 + ses_wc | school),
  data = soc,
  REML = TRUE,
  control = lme4::lmerControl(optimizer = "bobyqa")
)

soc_labels <- c(
  `(Intercept)` = "Intercept at average school SES and support",
  ses_wc = "Within-school effect of student SES",
  school_ses_gc = "Between-school effect of mean SES",
  school_support_gc = "Main effect of school support",
  `ses_wc:school_support_gc` = "Within-school SES × school support"
)
soc_fixed <- lmm_fixed_table(m_soc, soc_labels)
soc_empty_icc <- icc_random_intercept(m_soc_empty, "school")

soc_vc_raw <- as.data.frame(lme4::VarCorr(m_soc))
soc_vc_table <- data.frame(
  Group = c("School", "School", "School", "Residual"),
  Component = c(
    "Random intercept", "Random SES slope",
    "Intercept–slope correlation", "Student-level residual"
  ),
  Value = c(
    soc_vc_raw$sdcor[
      soc_vc_raw$grp == "school" & soc_vc_raw$var1 == "(Intercept)" &
        is.na(soc_vc_raw$var2)
    ],
    soc_vc_raw$sdcor[
      soc_vc_raw$grp == "school" & soc_vc_raw$var1 == "ses_wc" &
        is.na(soc_vc_raw$var2)
    ],
    soc_vc_raw$sdcor[
      soc_vc_raw$grp == "school" & !is.na(soc_vc_raw$var2)
    ],
    stats::sigma(m_soc)
  ),
  check.names = FALSE
)

knitr::kable(
  soc_fixed,
  digits = 3,
  caption = "Fixed effects from the student–school random-slope model"
)
Fixed effects from the student–school random-slope model
Term Estimate Standard error 95% CI lower 95% CI upper t value
(Intercept) Intercept at average school SES and support 70.88 0.644 69.62 72.143 110.06
ses_wc Within-school effect of student SES 2.77 0.260 2.26 3.274 10.64
school_ses_gc Between-school effect of mean SES 4.32 0.719 2.91 5.732 6.01
school_support_gc Main effect of school support 2.78 0.712 1.38 4.173 3.90
ses_wc:school_support_gc Within-school SES × school support -1.01 0.285 -1.57 -0.454 -3.55 
knitr::kable(
  soc_vc_table,
  digits = 3,
  caption = "Random-effect standard deviations, correlation, and residual from the student–school model"
)
Random-effect standard deviations, correlation, and residual from the student–school model
Group Component Value
School Random intercept 5.233
School Random SES slope 1.598
School Intercept–slope correlation 0.287
Residual Student-level residual 6.994

Within schools, a one-unit increase in SES is associated with an average achievement increase of 2.77 points; the between-school relationship for mean SES is 4.32 points. These are different quantities and cannot be replaced by one coefficient for raw SES. The cross-level interaction is -1.01: in this simulation, the positive within-school SES slope becomes weaker as school support increases. This describes effect modification; it does not prove that increasing support would necessarily alter the SES mechanism.

The school ICC from the empty model is 0.437. The random-slope standard deviation is 1.598, indicating between-school heterogeneity in the within-school SES relationship. The intercept–slope correlation is 0.287, and its interpretation depends on the zero point ses_wc=0.

soc_re <- lme4::ranef(m_soc)$school
selected_schools <- rownames(soc_re)[round(seq(1, nrow(soc_re), length.out = 12))]
school_profile <- unique(
  soc[c("school", "school_ses_gc", "school_support_gc")]
)
ses_grid <- seq(-2, 2, length.out = 80)
soc_beta <- lme4::fixef(m_soc)
school_predictions <- stats::setNames(
  lapply(selected_schools, function(s) {
    profile_s <- school_profile[school_profile$school == s, ]
    support_s <- profile_s$school_support_gc
    school_ses_s <- profile_s$school_ses_gc
    soc_beta["(Intercept)"] +
      soc_beta["school_ses_gc"] * school_ses_s +
      soc_beta["school_support_gc"] * support_s +
      soc_re[s, "(Intercept)"] +
      (soc_beta["ses_wc"] +
         soc_beta["ses_wc:school_support_gc"] * support_s +
         soc_re[s, "ses_wc"]) * ses_grid
  }),
  selected_schools
)
average_school_prediction <-
  soc_beta["(Intercept)"] + soc_beta["ses_wc"] * ses_grid
left_ylim <- grDevices::extendrange(
  c(unlist(school_predictions), average_school_prediction),
  f = 0.05
)
old_par <- par(mfrow = c(1, 2), mar = c(4.2, 4.2, 3, 1))

plot(
  NA,
  xlim = range(ses_grid),
  ylim = left_ylim,
  xlab = "Student SES minus school mean SES",
  ylab = "Conditional predicted achievement",
  main = "School-specific random slopes",
  las = 1
)
for (s in selected_schools) {
  lines(
    ses_grid,
    school_predictions[[s]],
    col = grDevices::adjustcolor(palette_ml["blue"], alpha.f = 0.45),
    lwd = 1.2
  )
}
lines(
  ses_grid,
  average_school_prediction,
  col = palette_ml["vermillion"],
  lwd = 2.6
)

support_sd <- stats::sd(
  soc$school_support_gc[!duplicated(soc$school)]
)
support_values <- c(-support_sd, 0, support_sd)
support_colors <- c(
  palette_ml["orange"], palette_ml["navy"], palette_ml["teal"]
)
support_predictions <- vapply(
  support_values,
  function(z) {
    soc_beta["(Intercept)"] + soc_beta["school_support_gc"] * z +
      (soc_beta["ses_wc"] +
         soc_beta["ses_wc:school_support_gc"] * z) * ses_grid
  },
  numeric(length(ses_grid))
)
right_ylim <- grDevices::extendrange(support_predictions, f = 0.05)
plot(
  NA,
  xlim = range(ses_grid),
  ylim = right_ylim,
  xlab = "Student SES minus school mean SES",
  ylab = "Fixed-part predicted achievement",
  main = "Cross-level modification by school support",
  las = 1
)
for (k in seq_along(support_values)) {
  lines(
    ses_grid,
    support_predictions[, k],
    col = support_colors[k],
    lwd = 2.4
  )
}
legend(
  "topleft",
  legend = c("Support = −1 SD", "Support = mean", "Support = +1 SD"),
  col = support_colors,
  lwd = 2.4,
  bty = "n"
)
Two panels. In the left panel, school-specific achievement lines have different intercepts and slopes. In the right panel, three fixed-prediction lines show a flatter SES slope at higher levels of school support.

Student–school random slopes and the cross-level interaction. The left panel shows conditional SES slopes for a subset of schools; the right panel shows fixed-part predicted relationships when school support is low, average, or high. 

par(old_par)

5.3 Correlation depends on covariate values under random slopes

When the centered SES values for two students in the same school are x1x_1 and x2x_2, the covariance induced by the random effects is:

Cov⁡(Y1,Y2)=τ00+(x1+x2)τ01+x1x2τ11. \operatorname{Cov}(Y_1,Y_2)= \tau_{00}+(x_1+x_2)\tau_{01}+x_1x_2\tau_{11}.

There is therefore no single ICC that is the same at every SES value.

G_soc <- as.matrix(lme4::VarCorr(m_soc)$school)
sigma2_soc <- stats::sigma(m_soc)^2

conditional_corr <- function(x1, x2, G, sigma2) {
  z1 <- c(1, x1)
  z2 <- c(1, x2)
  covariance <- drop(t(z1) %*% G %*% z2)
  variance1 <- drop(t(z1) %*% G %*% z1) + sigma2
  variance2 <- drop(t(z2) %*% G %*% z2) + sigma2
  covariance / sqrt(variance1 * variance2)
}

soc_corr_table <- data.frame(
  `Student 1 SES` = c(0, -1),
  `Student 2 SES` = c(0, 1),
  `Model-implied within-school correlation` = c(
    conditional_corr(0, 0, G_soc, sigma2_soc),
    conditional_corr(-1, 1, G_soc, sigma2_soc)
  ),
  check.names = FALSE
)

knitr::kable(
  soc_corr_table,
  digits = 3,
  caption = "Within-school correlations at selected SES combinations under the random-slope model"
)
Within-school correlations at selected SES combinations under the random-slope model
Student 1 SES Student 2 SES Model-implied within-school correlation
0 0 0.359
-1 1 0.315
Check your understanding: Why can a random intercept not automatically substitute for a random slope? A random intercept allows school mean levels to differ but still forces the SES slope to be identical across schools. If the scientific question and repeated information support slope heterogeneity, including only a random intercept incorrectly restricts the covariance structure and can understate uncertainty in the fixed slope.

6 Psychology application: occasions nested in people nested in therapists

6.1 Three-level longitudinal data

One hundred sixty participants are treated by 20 therapists, with seven therapy sessions recorded for each participant. Time begins at zero so the intercept represents baseline symptoms. Stress is decomposed into a within-person deviation, stress_wp, and a person’s mean stress, stress_pm_gc.

set.seed(20260823)
K <- 20L
therapists <- sprintf("T%02d", seq_len(K))
therapist_u <- rnorm(K, 0, 1.2)
persons_per_therapist <- 8L
session0 <- 0:6
psy_list <- list()
person_index <- 0L

for (k in seq_len(K)) {
  intervention <- sample(rep(0:1, each = persons_per_therapist / 2))
  for (p in seq_len(persons_per_therapist)) {
    person_index <- person_index + 1L
    z0p <- rnorm(1)
    z1p <- rnorm(1)
    p_u0 <- 3.0 * z0p
    p_u1 <- 0.45 * (-0.20 * z0p + sqrt(1 - 0.20^2) * z1p)
    stress_pm <- rnorm(1, 0, 0.9)
    stress_wp <- rnorm(length(session0), 0, 0.8)
    stress_wp <- stress_wp - mean(stress_wp)
    symptoms <- 22 - 0.90 * session0 -
      0.20 * intervention[p] - 0.35 * session0 * intervention[p] +
      1.20 * stress_wp + 2.00 * stress_pm +
      therapist_u[k] + p_u0 + p_u1 * session0 +
      rnorm(length(session0), 0, 2.5)
    psy_list[[person_index]] <- data.frame(
      therapist = therapists[k],
      person = sprintf("P%03d", person_index),
      session0 = session0,
      intervention = intervention[p],
      stress = stress_pm + stress_wp,
      symptoms = symptoms
    )
  }
}
psy <- do.call(rbind, psy_list)
psy$stress_pm <- ave(psy$stress, psy$person, FUN = mean)
psy$stress_wp <- psy$stress - psy$stress_pm
psy$stress_pm_gc <- psy$stress_pm - mean(psy$stress_pm)
psy$intervention <- factor(
  psy$intervention,
  0:1,
  c("Control", "Intervention")
)

psy_audit <- data.frame(
  Observations = nrow(psy),
  People = length(unique(psy$person)),
  Therapists = length(unique(psy$therapist)),
  `Occasions per person` = paste(range(table(psy$person)), collapse = "–"),
  `People per therapist` = paste(
    range(table(unique(psy[c("person", "therapist")])$therapist)),
    collapse = "–"
  ),
  check.names = FALSE
)

knitr::kable(psy_audit, caption = "Audit of the three-level psychology longitudinal data")
Audit of the three-level psychology longitudinal data
Observations People Therapists Occasions per person People per therapist
1120 160 20 7–7 8–8

6.2 Three-level random intercepts and a person-specific time slope

Ytij=β0+β1Timetij+β2Aij+β3TimetijAij+βWStresswithin+βBStressbetween+v0j+u0ij+u1ijTimetij+εtij. Y_{tij}=\beta_0+\beta_1Time_{tij}+\beta_2A_{ij}+ \beta_3Time_{tij}A_{ij}+\beta_WStress_{within}+ \beta_BStress_{between}+v_{0j}+u_{0ij}+u_{1ij}Time_{tij}+\varepsilon_{tij}. 

m_psy <- lme4::lmer(
  symptoms ~ session0 * intervention + stress_wp + stress_pm_gc +
    (1 + session0 | person) + (1 | therapist),
  data = psy,
  REML = TRUE,
  control = lme4::lmerControl(optimizer = "bobyqa")
)

psy_labels <- c(
  `(Intercept)` = "Baseline symptoms in the control group",
  session0 = "Average change per session in the control group",
  interventionIntervention = "Intervention-group difference at baseline",
  stress_wp = "Within-person effect of stress",
  stress_pm_gc = "Between-person effect of mean stress",
  `session0:interventionIntervention` = "Additional per-session change in the intervention group"
)
psy_fixed <- lmm_fixed_table(m_psy, psy_labels)

psy_vc_raw <- as.data.frame(lme4::VarCorr(m_psy))
psy_vc_table <- data.frame(
  Component = c(
    "Person random-intercept SD", "Person random time-slope SD",
    "Intercept–slope correlation", "Therapist intercept SD", "Residual SD"
  ),
  Value = c(
    psy_vc_raw$sdcor[
      psy_vc_raw$grp == "person" & psy_vc_raw$var1 == "(Intercept)" &
        is.na(psy_vc_raw$var2)
    ],
    psy_vc_raw$sdcor[
      psy_vc_raw$grp == "person" & psy_vc_raw$var1 == "session0" &
        is.na(psy_vc_raw$var2)
    ],
    psy_vc_raw$sdcor[
      psy_vc_raw$grp == "person" & !is.na(psy_vc_raw$var2)
    ],
    psy_vc_raw$sdcor[
      psy_vc_raw$grp == "therapist" & is.na(psy_vc_raw$var2)
    ],
    stats::sigma(m_psy)
  ),
  check.names = FALSE
)

knitr::kable(
  psy_fixed,
  digits = 3,
  caption = "Fixed effects from the three-level occasion–person–therapist model"
)
Fixed effects from the three-level occasion–person–therapist model
Term Estimate Standard error 95% CI lower 95% CI upper t value
(Intercept) Baseline symptoms in the control group 22.335 0.524 21.307 23.363 42.59
session0 Average change per session in the control group -0.857 0.062 -0.979 -0.735 -13.74
interventionIntervention Intervention-group difference at baseline -0.621 0.531 -1.661 0.419 -1.17
stress_wp Within-person effect of stress 1.326 0.104 1.123 1.529 12.78
stress_pm_gc Between-person effect of mean stress 1.654 0.289 1.088 2.220 5.73
session0:interventionIntervention Additional per-session change in the intervention group -0.403 0.088 -0.576 -0.230 -4.57 
knitr::kable(
  psy_vc_table,
  digits = 3,
  caption = "Random effects and residual from the three-level psychology model"
)
Random effects and residual from the three-level psychology model
Component Value
Person random-intercept SD 2.869
Person random time-slope SD 0.295
Intercept–slope correlation 0.043
Therapist intercept SD 1.642
Residual SD 2.508

Symptoms change by an average of -0.86 points per session in the control group; the intervention group has an additional per-session change of -0.4 points. On a day when stress is one unit above a person’s usual level, symptoms are higher by an average of 1.33 points. People whose long-term mean stress differs by one unit differ by 1.65 points on average. These two coefficients answer different questions.

Only 20 therapists are represented, so the therapist variance and any therapist ranking would have substantial uncertainty. A therapist random intercept represents dependence, but it does not show that between-therapist differences are caused by treatment quality.

selected_people <- unique(psy$person)[round(seq(1, length(unique(psy$person)), length.out = 12))]
psy_selected <- psy[psy$person %in% selected_people, ]
psy_beta <- lme4::fixef(m_psy)

plot(
  NA,
  xlim = range(psy$session0),
  ylim = range(psy_selected$symptoms),
  xlab = "Session (0 = baseline)",
  ylab = "Symptom score",
  main = "Individual and fixed-part mean trajectories",
  las = 1
)
for (id in selected_people) {
  d <- psy_selected[psy_selected$person == id, ]
  color <- if (d$intervention[1] == "Intervention") {
    grDevices::adjustcolor(palette_ml["teal"], alpha.f = 0.45)
  } else {
    grDevices::adjustcolor(palette_ml["orange"], alpha.f = 0.45)
  }
  lines(d$session0, d$symptoms, col = color, lwd = 1.2)
  points(d$session0, d$symptoms, col = color, pch = 16, cex = 0.45)
}
session_grid <- 0:6
fixed_control <- psy_beta["(Intercept)"] + psy_beta["session0"] * session_grid
fixed_intervention <- psy_beta["(Intercept)"] +
  psy_beta["interventionIntervention"] +
  (psy_beta["session0"] +
     psy_beta["session0:interventionIntervention"]) * session_grid
lines(session_grid, fixed_control, col = palette_ml["orange"], lwd = 3)
lines(session_grid, fixed_intervention, col = palette_ml["teal"], lwd = 3)
legend(
  "topright",
  legend = c("Control fixed trajectory", "Intervention fixed trajectory"),
  col = c(palette_ml["orange"], palette_ml["teal"]),
  lwd = 3,
  bty = "n"
)
Twelve individual symptom trajectories generally decline across seven sessions but have different slopes. Two thick lines show a faster average decline in the intervention group than in the control group.

Symptom trajectories for a subset of participants and the model’s fixed-part mean trajectories. Thin lines are observed individual trajectories; thick lines compare fixed-part change for control and intervention groups with stress held at its mean.

A random time slope allows people to change at different rates, but it does not automatically remove AR(1) correlation between residuals from adjacent sessions. If the residuals retain a clear time-series structure, use a method that supports the relevant correlation structure and conduct sensitivity analyses.

7 Medical application II: a binary outcome for patients nested in hospitals

7.1 Why a logistic GLMM is needed

For a binary infection outcome:

Yij∣uj∼Bernoulli(pij),logit⁡(pij)=XijTβ+uj. Y_{ij}\mid u_j\sim Bernoulli(p_{ij}),\qquad \operatorname{logit}(p_{ij})=X_{ij}^{T}\beta+u_j.

Here, exp⁡(β)\exp(\beta) is a cluster-specific conditional OR given the hospital random effect and model covariates. It generally differs from a marginal OR averaged over the distribution of hospitals.

set.seed(20260829)
J <- 56L
n_j <- sample(36:54, J, replace = TRUE)
hospital <- sprintf("H%02d", seq_len(J))
checklist_j <- sample(rep(0:1, each = J / 2))
h_u <- rnorm(J, 0, 0.65)

bin <- do.call(rbind, lapply(seq_len(J), function(j) {
  n <- n_j[j]
  age <- pmin(pmax(rnorm(n, 62, 12), 18), 90)
  emergency <- rbinom(n, 1, 0.28)
  eta <- -2.25 - 0.55 * checklist_j[j] +
    0.70 * emergency + 0.22 * ((age - 60) / 10) + h_u[j]
  data.frame(
    hospital = hospital[j],
    checklist = checklist_j[j],
    emergency = emergency,
    age_c10 = (age - 60) / 10,
    infection = rbinom(n, 1, plogis(eta))
  )
}))
bin$checklist <- factor(bin$checklist, 0:1, c("Usual care", "Checklist"))
bin$emergency <- factor(bin$emergency, 0:1, c("Elective", "Emergency"))

bin_audit <- data.frame(
  Patients = nrow(bin),
  Hospitals = length(unique(bin$hospital)),
  `Infection events` = sum(bin$infection),
  `Infection proportion` = pct(mean(bin$infection)),
  `Hospital size range` = paste(range(table(bin$hospital)), collapse = "–"),
  check.names = FALSE
)

knitr::kable(bin_audit, caption = "Audit of the simulated binary hospital-infection data")
Audit of the simulated binary hospital-infection data
Patients Hospitals Infection events Infection proportion Hospital size range
2501 56 298 11.9% 36–53 
hospital_rates <- aggregate(infection ~ hospital, data = bin, FUN = mean)
hospital_rates$n <- as.numeric(table(bin$hospital)[hospital_rates$hospital])
hospital_rates <- hospital_rates[order(hospital_rates$infection), ]

plot(
  hospital_rates$infection,
  seq_len(nrow(hospital_rates)),
  pch = 16,
  cex = 0.7 + 1.2 * sqrt(hospital_rates$n / max(hospital_rates$n)),
  col = palette_ml["blue"],
  yaxt = "n",
  xlab = "Observed infection proportion",
  ylab = "Hospitals ordered by infection proportion",
  main = "Raw hospital infection proportions",
  las = 1
)
abline(v = mean(bin$infection), lty = 2, col = palette_ml["vermillion"])
Points for fifty-six hospital infection proportions are ordered by value and distributed around a dashed overall infection-rate line. Several small hospitals appear near the extremes.

Observed infection proportions by hospital and the overall proportion. Hospital sample sizes and case mix differ, so raw proportions should not be interpreted directly as hospital quality rankings.

7.2 Conditional OR, latent ICC, and median odds ratio

m_bin <- lme4::glmer(
  infection ~ checklist + emergency + age_c10 + (1 | hospital),
  data = bin,
  family = stats::binomial,
  nAGQ = 1,
  control = lme4::glmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 2e5)
  )
)

bin_vc <- as.data.frame(lme4::VarCorr(m_bin))
tau2_bin <- bin_vc$vcov[
  bin_vc$grp == "hospital" &
    bin_vc$var1 == "(Intercept)" &
    is.na(bin_vc$var2)
]
stopifnot(length(tau2_bin) == 1L)
icc_latent <- tau2_bin / (tau2_bin + pi^2 / 3)
mor <- exp(stats::qnorm(0.75) * sqrt(2 * tau2_bin))
pearson <- residuals(m_bin, type = "pearson")
pearson_ratio <- sum(pearson^2) /
  (stats::nobs(m_bin) - length(lme4::fixef(m_bin)))

bin_coef <- coef(summary(m_bin))
bin_or <- wald_or(m_bin)
bin_labels <- c(
  `(Intercept)` = "Usual care, elective admission, age 60",
  checklistChecklist = "Infection checklist vs usual care",
  emergencyEmergency = "Emergency vs elective admission",
  age_c10 = "Per 10-year increase in age"
)
bin_results <- data.frame(
  Term = unname(bin_labels[rownames(bin_coef)]),
  `Conditional OR` = bin_or[, "OR"],
  `95% CI lower` = bin_or[, "lower"],
  `95% CI upper` = bin_or[, "upper"],
  `p value` = vapply(bin_coef[, "Pr(>|z|)"], format_p, character(1)),
  check.names = FALSE
)

bin_heterogeneity <- data.frame(
  Metric = c(
    "Hospital random-intercept variance", "Latent-scale ICC",
    "Median odds ratio", "Pearson ratio"
  ),
  Value = c(tau2_bin, icc_latent, mor, pearson_ratio),
  check.names = FALSE
)

knitr::kable(
  bin_results,
  digits = 3,
  caption = "Conditional ORs and Wald intervals from the hospital-infection logistic GLMM"
)
Conditional ORs and Wald intervals from the hospital-infection logistic GLMM
Term Conditional OR 95% CI lower 95% CI upper p value
(Intercept) Usual care, elective admission, age 60 0.126 0.094 0.169 <0.001
checklistChecklist Infection checklist vs usual care 0.571 0.380 0.856 0.00676
emergencyEmergency Emergency vs elective admission 1.562 1.197 2.037 0.001
age_c10 Per 10-year increase in age 1.254 1.130 1.393 <0.001 
knitr::kable(
  bin_heterogeneity,
  digits = 3,
  caption = "Summary of hospital heterogeneity and binary-model diagnostics"
)
Summary of hospital heterogeneity and binary-model diagnostics
Metric Value
Hospital random-intercept variance 0.347
Latent-scale ICC 0.095
Median odds ratio 1.754
Pearson ratio 0.905

The model explicitly uses nAGQ = 1, the Laplace approximation used by lme4::glmer(). A formal analysis should report the integration approximation and optimizer. With rare events, very small clusters, or large random effects, increase the number of adaptive Gauss–Hermite quadrature points and examine sensitivity.

The checklist’s conditional OR is 0.571 (95% CI 0.38 to 0.856). Among patients of the same age and emergency status and with the same latent hospital propensity, the odds of infection are lower in the checklist group.

For a simple logistic random-intercept model, the latent ICC is:

ICClatent≈τ00τ00+π2/3. ICC_{latent}\approx\frac{\tau_{00}}{\tau_{00}+\pi^2/3}.

In this example it is 0.095. It belongs to a latent logistic scale and cannot be read directly as a percentage of the observed 0/1 infection variance. The median odds ratio is 1.754: if the same patient were moved from a randomly selected lower-risk hospital to a randomly selected higher-risk hospital, the median multiplicative hospital effect on the odds would be approximately this value.

The Pearson ratio is 0.905 and does not suggest obvious overdispersion, but a single ratio cannot exclude incorrect functional form, zero inflation, or an omitted level.

7.3 Setting the hospital effect to zero is not the same as integrating over hospital heterogeneity

beta <- lme4::fixef(m_bin)
eta_usual <- unname(beta["(Intercept)"])
eta_check <- eta_usual + unname(beta["checklistChecklist"])
p_typical <- stats::plogis(c(usual = eta_usual, checklist = eta_check))

integrand <- function(u, eta, sd) {
  stats::plogis(eta + u) * stats::dnorm(u, 0, sd)
}

p_marginal <- c(
  usual = stats::integrate(
    integrand,
    -Inf,
    Inf,
    eta = eta_usual,
    sd = sqrt(tau2_bin)
  )$value,
  checklist = stats::integrate(
    integrand,
    -Inf,
    Inf,
    eta = eta_check,
    sd = sqrt(tau2_bin)
  )$value
)

probability_table <- data.frame(
  Scale = c(
    "Reference patient, typical hospital: u=0",
    "Reference patient, integrated over hospital random effect"
  ),
  `Usual-care risk` = c(p_typical["usual"], p_marginal["usual"]),
  `Checklist risk` = c(p_typical["checklist"], p_marginal["checklist"]),
  `Risk difference` = c(diff(p_typical), diff(p_marginal)),
  check.names = FALSE
)

knitr::kable(
  probability_table,
  digits = 3,
  caption = "Infection probabilities for the reference patient at a typical hospital and after integration over the hospital random effect"
)
Infection probabilities for the reference patient at a typical hospital and after integration over the hospital random effect
Scale Usual-care risk Checklist risk Risk difference
Reference patient, typical hospital: u=0 0.112 0.067 -0.045
Reference patient, integrated over hospital random effect 0.125 0.077 -0.048

predict(m_bin, re.form = NA, type = "response") sets the hospital random effect to zero and returns the fixed-part probability for a typical hospital; it does not automatically integrate over the distribution of hospitals. The second row of the table integrates only over the hospital random intercept and still refers to a 60-year-old patient with an elective admission. Estimating an overall marginal risk for a target population would additionally require standardization over that population’s age, emergency-status, and other covariate distributions.

prob_matrix <- rbind(
  `Reference patient: typical hospital u=0` = p_typical,
  `Reference patient: hospital effect integrated` = p_marginal
)
barplot(
  t(prob_matrix),
  beside = TRUE,
  col = c(palette_ml["orange"], palette_ml["teal"]),
  ylim = c(0, max(prob_matrix) * 1.25),
  ylab = "Predicted infection probability",
  main = "Fixed-part and integrated probabilities",
  las = 1
)
legend(
  "topright",
  legend = c("Usual care", "Checklist"),
  fill = c(palette_ml["orange"], palette_ml["teal"]),
  bty = "n"
)
A grouped bar chart compares usual care with the checklist. Both approaches show lower risk with the checklist; both bars after integration over the hospital random effect are slightly higher than the corresponding typical-hospital probabilities.

Infection probabilities for a 60-year-old reference patient with an elective admission, at a typical hospital and after integration over the hospital random effect. The nonlinear logit link makes the integrated result differ from the probability obtained by setting the random effect to zero.

Check your understanding: Can a conditional OR of 0.57 be reported as a 43% reduction in infection risk? No. An OR compares odds, not risks, and this is a conditional OR given the hospital random effect. Separately calculate standardized risks, risk ratios, or risk differences for the target population and its covariate distribution, and state whether the calculation integrates over random effects.

8 More complex cluster structures

8.1 Nesting, cross-classification, and multiple membership require different syntax

If schools are strictly nested in communities, (1 | community/school) may be appropriate. If students’ residential communities and schools cross one another, estimate separate random-intercept terms. Forcing a crossed structure into a nested model assigns variance to the wrong social structure.

# Schools are strictly nested in communities; school_id need only be unique
# within a community.
m_nested <- lme4::lmer(
  score ~ student_ses + school_resources +
    (1 | community_id/school_id),
  data = education_data
)

# Schools and residential communities are cross-classified.
m_crossed <- lme4::lmer(
  score ~ student_ses + school_resources + neighborhood_deprivation +
    (1 | school_id) + (1 | neighborhood_id),
  data = education_data
)

# Multiple membership: if several therapists jointly treat one patient, the
# membership weight for each therapist must be specified. A single ordinary
# (1 | therapist) term cannot represent this structure.

Multiple-membership models generally require weights for every upper-level member and software that supports the structure. If a patient changes therapists, retaining only the “last therapist” discards the actual dependence process.

9 ML, REML, and inference

9.1 Match the estimation method to the comparison

fitting_reference <- data.frame(
  Situation = c(
    "Final LMM variance components and fixed effects",
    "Comparing LMMs with different fixed effects",
    "Comparing covariance structures with the same fixed effects",
    "Logistic or Poisson GLMM",
    "Few upper-level clusters"
  ),
  Recommendation = c(
    "Commonly use REML and report both fixed and random parts",
    "Refit with ML before comparison",
    "REML may be used, but boundary tests need care",
    "Use maximum likelihood; state the integration approximation and optimizer",
    "Emphasize effects and intervals; consider bootstrap or small-sample corrections"
  ),
  Reason = c(
    "REML reduces finite-sample bias in variance components",
    "REML likelihoods with different fixed effects are not directly comparable",
    "Zero random variance lies on the boundary of the parameter space",
    "Non-Gaussian models do not have ordinary LMM REML",
    "Asymptotic z or t approximations may be too optimistic"
  ),
  check.names = FALSE
)

knitr::kable(fitting_reference, caption = "Choosing ML, REML, and uncertainty methods")
Choosing ML, REML, and uncertainty methods
Situation Recommendation Reason
Final LMM variance components and fixed effects Commonly use REML and report both fixed and random parts REML reduces finite-sample bias in variance components
Comparing LMMs with different fixed effects Refit with ML before comparison REML likelihoods with different fixed effects are not directly comparable
Comparing covariance structures with the same fixed effects REML may be used, but boundary tests need care Zero random variance lies on the boundary of the parameter space
Logistic or Poisson GLMM Use maximum likelihood; state the integration approximation and optimizer Non-Gaussian models do not have ordinary LMM REML
Few upper-level clusters Emphasize effects and intervals; consider bootstrap or small-sample corrections Asymptotic z or t approximations may be too optimistic

lme4::lmer() does not provide fixed-effect p values by default because finite-sample degrees of freedom are not uniquely defined. This tutorial describes results with estimates, Wald 95% CIs, and t values. A formal study can prespecify Satterthwaite or Kenward–Roger approximations, profile likelihood, or a parametric bootstrap. Do not choose whichever interval method looks most favorable after seeing the results.

The null value for a random-effect variance lies on the boundary of the parameter space, so the usual chi-squared approximation for a likelihood-ratio test may be inaccurate. The random structure should primarily follow the data levels, repeated-measures arrangement, and scientific question—not a stepwise search for the smallest p value.

More lower-level observations cannot substitute for more upper-level clusters Adding patients to each hospital improves estimation of that hospital’s mean, but precision for hospital-level policies, treatments, or resources is driven mainly by the number of hospitals. Whenever reporting the total row count, also report the number of units at every level and the distribution of cluster sizes.

9.2 Conditional predictions, a typical cluster, and new clusters

# Conditional prediction for observed people and clusters: include their
# estimated random effects.
predict(m_psy, newdata = existing_people, re.form = NULL)

# Set all random effects to zero: fixed-part or typical-cluster prediction.
predict(m_psy, newdata = prediction_grid, re.form = NA)

# A new person treated by an existing therapist: retain the therapist effect
# but omit the unknown person effect.
predict(
  m_psy,
  newdata = new_people_existing_therapist,
  re.form = ~(1 | therapist),
  allow.new.levels = TRUE
)

# A new person and a new therapist: both random effects are unknown, so use
# the fixed part.
predict(
  m_psy,
  newdata = new_people_new_therapist,
  re.form = NA,
  allow.new.levels = TRUE
)

re.form = NA removes every random effect in the model; it cannot represent the partially known structure “new person, existing therapist.” Retain only the therapist term as in the example for that scenario. With an identity-link LMM, averaging over zero-mean random effects aligns the mean with the fixed part. With a nonlinear-link GLMM, setting random effects to zero generally differs from integrating over their distribution. Predictions should also distinguish a confidence interval for the mean, a conditional prediction interval for an existing person, and a prediction interval for a new cluster.

10 Diagnostics: a model that runs is not necessarily trustworthy

10.1 Checking all four primary models consistently

diagnostic_flags <- data.frame(
  Model = c(
    "Medical continuous-outcome LMM", "Sociology random-slope LMM",
    "Psychology three-level LMM", "Medical infection GLMM"
  ),
  `Singular fit` = c(
    lme4::isSingular(m_med),
    lme4::isSingular(m_soc),
    lme4::isSingular(m_psy),
    lme4::isSingular(m_bin)
  ),
  `Optimizer reports convergence` = vapply(
    list(m_med, m_soc, m_psy, m_bin),
    function(model) {
      optimizer_code <- model@optinfo$conv$opt
      all(optimizer_code == 0) &&
        is.null(model@optinfo$conv$lme4$messages)
    },
    logical(1)
  ),
  check.names = FALSE
)

knitr::kable(diagnostic_flags, caption = "Singularity and convergence flags for the four teaching models")
Singularity and convergence flags for the four teaching models
Model Singular fit Optimizer reports convergence
Medical continuous-outcome LMM FALSE TRUE
Sociology random-slope LMM FALSE TRUE
Psychology three-level LMM FALSE TRUE
Medical infection GLMM FALSE TRUE

All four teaching models converge successfully and are not singular fits, but that is only a minimum requirement. A complete diagnostic review also asks whether:

  • IDs, repeated records, and cluster sizes are correct at every level;
  • a continuous predictor varies enough within clusters to support its random slope;
  • the fixed part has nonlinearity, omitted interactions, or heteroskedasticity;
  • influential observations or clusters occur at either level;
  • the random-effect distribution has severe departures;
  • a binary model has overdispersion, zero inflation, or an omitted level;
  • longitudinal residuals retain serial correlation; and
  • a few exceptionally large or small clusters drive the conclusion.

10.1.1 Convergence warnings and singular fits

A convergence warning can result from differences in scale, gradients, an overly complex structure, or weak information. First check the data, centering or scaling, within-cluster variation, and model design; then treat an alternate optimizer as a numerical sensitivity check. A warning’s disappearance does not resolve the scientific problem.

A singular fit commonly means that a random variance is near zero or a correlation is near ±1, placing a dimension at the boundary under the current data. It is not an automatic instruction to remove every random slope. Make a transparent decision based on the prespecified question, design support, alternative reasonable structures, and sensitivity results.

11 Boundaries involving missingness, selection, and causal interpretation

11.1 A mixed model does not automatically repair missing data

A maximum-likelihood mixed model can use individuals with different numbers of observed occasions, but unbiased interpretation still depends on the missingness mechanism and model. Standard likelihood inference commonly assumes that missingness is reasonably considered missing at random given the observed history and model covariates.

  • If people with worsening symptoms are more likely to stop completing diaries, an ordinary model may be biased.
  • Withdrawal of an entire hospital and loss to follow-up of one patient are selection processes at different levels.
  • Missing-not-at-random dropout requires pattern-mixture, selection-model, or other sensitivity analyses.
  • Single mean imputation distorts variance and multilevel relationships.
  • If cluster size is associated with potential outcomes, informative cluster size may be present.

11.2 Representing dependence does not identify a causal effect

A multilevel model can represent correlation, but it does not automatically provide:

  • a well-defined intervention strategy and common time zero;
  • conditional exchangeability or absence of unmeasured confounding;
  • positivity;
  • consistency and absence of interference;
  • appropriate treatment of time-varying treatment and time-varying confounding; or
  • a remedy for selection into institutions, loss to follow-up, or measurement error.

A hospital-level intervention particularly requires control of hospital-level confounding and adequate overlap across hospitals. A statistically clear patient-level coefficient should not automatically be translated into the causal effect of changing that exposure for a patient. See Causal Inference in Depth for the full framework.

A random intercept is not an “unmeasured-confounding controller” A random intercept describes a distribution of unobserved cluster deviations. It neither guarantees that those deviations are independent of model covariates nor automatically controls all shared cluster-level causes. When random effects may be associated with covariates, consider a correlated-random-effects or Mundlak decomposition that includes cluster means, and state the additional assumptions.

12 Integrating results across the four applications

12.1 Dynamic results summary

case_summary <- data.frame(
  Application = c(
    "Patients–clinics", "Students–schools",
    "Occasions–people–therapists", "Patients–hospital infection"
  ),
  `Data structure` = c(
    sprintf("%d people / %d clinics", nrow(med), length(unique(med$clinic))),
    sprintf("%d people / %d schools", nrow(soc), length(unique(soc$school))),
    sprintf(
      "%d occasions / %d people / %d therapists",
      nrow(psy), length(unique(psy$person)), length(unique(psy$therapist))
    ),
    sprintf("%d people / %d hospitals", nrow(bin), length(unique(bin$hospital)))
  ),
  `Primary result` = c(
    sprintf("Intervention mean difference %.2f mmHg", med_fixed$Estimate[2]),
    sprintf(
      "Within-school SES %.2f; cross-level interaction %.2f",
      soc_fixed$Estimate[2], soc_fixed$Estimate[5]
    ),
    sprintf("Additional intervention session slope %.2f", psy_fixed$Estimate[6]),
    sprintf(
      "Checklist conditional OR %.2f; reference-patient risk difference after hospital-effect integration %.3f",
      bin_results[["Conditional OR"]][2], diff(p_marginal)
    )
  ),
  `Key heterogeneity` = c(
    sprintf("Empty-model ICC %.3f", med_empty_icc),
    sprintf(
      "Empty-model ICC %.3f; SES slope SD %.2f",
      soc_empty_icc, soc_vc_table$Value[2]
    ),
    sprintf(
      "Therapist SD %.2f; person slope SD %.2f",
      psy_vc_table$Value[4], psy_vc_table$Value[2]
    ),
    sprintf("Latent ICC %.3f; MOR %.2f", icc_latent, mor)
  ),
  check.names = FALSE
)

knitr::kable(case_summary, caption = "Dynamic summary of the four multilevel-model applications")
Dynamic summary of the four multilevel-model applications
Application Data structure Primary result Key heterogeneity
Patients–clinics 1425 people / 48 clinics Intervention mean difference -4.11 mmHg Empty-model ICC 0.174
Students–schools 2098 people / 70 schools Within-school SES 2.77; cross-level interaction -1.01 Empty-model ICC 0.437; SES slope SD 1.60
Occasions–people–therapists 1120 occasions / 160 people / 20 therapists Additional intervention session slope -0.40 Therapist SD 1.64; person slope SD 0.29
Patients–hospital infection 2501 people / 56 hospitals Checklist conditional OR 0.57; reference-patient risk difference after hospital-effect integration -0.048 Latent ICC 0.095; MOR 1.75

These results use different estimands and scales and should not be compared by numerical magnitude. The continuous-outcome models report mean differences and slopes, whereas the logistic GLMM reports a conditional OR. ICC definitions also differ between random-slope and non-Gaussian models.

13 Auditable reporting template

We analyzed [N] observations from [P] people and [J] clusters. The outcome was [definition and scale], and time was centered at [zero point]. The fixed part included [variables, nonlinear terms, and interactions], and the random part included [cluster intercepts, person intercepts, random slopes, and correlations]. The continuous outcome was fit with [REML/ML], and models were compared with [method]; the non-Gaussian outcome used a [link and family]. Cluster, person, and residual standard deviations were [values], and the [ICC/VPC/MOR] was [value with explicit definition]. The primary [mean difference/slope/conditional OR/marginal risk difference] was [estimate and 95% CI]. Diagnostics showed [convergence, singularity, residual behavior, influential clusters, and overdispersion]. Inference primarily generalizes to [target people and clusters] and depends on [random-effect, missingness, selection, and causal assumptions].

At minimum, report all of the following:

  • unit counts at every level, the range of cluster sizes, and imbalance;
  • each predictor’s level, centering convention, and zero point;
  • the complete fixed and random formulas;
  • random-effect standard deviations, correlations, and the residual scale;
  • whether effects are conditional or marginal and the scale on which they are expressed;
  • ML or REML, optimizer, interval method, and degrees-of-freedom method;
  • convergence, singularity, residual, and influential-cluster diagnostics; and
  • limitations involving missingness, serial correlation, unmeasured confounding, and extrapolation.

14 Common errors: quick guide

common_errors <- data.frame(
  Error = c(
    "Treating every row as independent after ignoring clusters",
    "Mechanically adding a random intercept whenever an ID appears",
    "Interpreting fixed and random as whether a variable changes",
    "Ignoring clustering because the ICC is small",
    "Calling the ICC the cluster's causal contribution",
    "Using the person count to assess information for an upper-level effect",
    "Treating an ID unique only within clusters as globally unique",
    "Writing a crossed structure incorrectly as nested",
    "Failing to decompose within- and between-person time-varying effects",
    "Assuming grand-mean centering has separated the two effects",
    "Using only a random intercept when a random slope is needed",
    "Treating a random slope as an AR(1) residual structure",
    "Comparing REML likelihoods across different fixed effects",
    "Selecting the random structure solely by stepwise p values",
    "Ignoring convergence or singular-fit warnings",
    "Treating a GLMM conditional OR as a risk ratio or marginal OR",
    "Using random effects to give institutions definitive ranks",
    "Assuming a mixed model automatically handles all missingness",
    "Treating a random intercept as an unmeasured-confounding controller"
  ),
  `Better practice` = c(
    "Represent dependence according to the design and audit sample size at every level",
    "First explain the shared mechanism and unit of inference",
    "Describe common coefficients and a distribution of cluster deviations",
    "Consider cluster size, cluster count, and the level of the exposure",
    "Describe it as a model variance or correlation structure",
    "Recognize that upper-level precision mainly comes from upper-level units",
    "Create a globally unique ID or an explicit interaction ID",
    "Model school and community random effects separately",
    "Include the within-cluster deviation and cluster mean",
    "Choose a group-mean decomposition that matches the question",
    "Fit and diagnose a random slope when the design supports it",
    "Separately assess remaining serial correlation",
    "Refit with ML before comparing fixed parts",
    "Prespecify candidate structures from design and scientific questions",
    "Investigate scaling, information, and complexity and report transparently",
    "Report the correct scale and estimate the target marginal quantity",
    "Show shrinkage, uncertainty, and case mix",
    "State assumptions such as MAR and conduct MNAR sensitivity analyses",
    "State additional assumptions such as random-effect–covariate independence"
  ),
  check.names = FALSE
)

knitr::kable(common_errors, caption = "Common multilevel-model errors and better practices")
Common multilevel-model errors and better practices
Error Better practice
Treating every row as independent after ignoring clusters Represent dependence according to the design and audit sample size at every level
Mechanically adding a random intercept whenever an ID appears First explain the shared mechanism and unit of inference
Interpreting fixed and random as whether a variable changes Describe common coefficients and a distribution of cluster deviations
Ignoring clustering because the ICC is small Consider cluster size, cluster count, and the level of the exposure
Calling the ICC the cluster’s causal contribution Describe it as a model variance or correlation structure
Using the person count to assess information for an upper-level effect Recognize that upper-level precision mainly comes from upper-level units
Treating an ID unique only within clusters as globally unique Create a globally unique ID or an explicit interaction ID
Writing a crossed structure incorrectly as nested Model school and community random effects separately
Failing to decompose within- and between-person time-varying effects Include the within-cluster deviation and cluster mean
Assuming grand-mean centering has separated the two effects Choose a group-mean decomposition that matches the question
Using only a random intercept when a random slope is needed Fit and diagnose a random slope when the design supports it
Treating a random slope as an AR(1) residual structure Separately assess remaining serial correlation
Comparing REML likelihoods across different fixed effects Refit with ML before comparing fixed parts
Selecting the random structure solely by stepwise p values Prespecify candidate structures from design and scientific questions
Ignoring convergence or singular-fit warnings Investigate scaling, information, and complexity and report transparently
Treating a GLMM conditional OR as a risk ratio or marginal OR Report the correct scale and estimate the target marginal quantity
Using random effects to give institutions definitive ranks Show shrinkage, uncertainty, and case mix
Assuming a mixed model automatically handles all missingness State assumptions such as MAR and conduct MNAR sensitivity analyses
Treating a random intercept as an unmeasured-confounding controller State additional assumptions such as random-effect–covariate independence

15 Exercises and answers

  1. An empty hospital model has ICC=0.08. Can you write that “hospitals caused 8% of infections”?
  2. Which design better supports an effect of school-level resources: 50 schools with 20 students each or five schools with 200 students each?
  3. Which two questions are mixed when diary stress enters a model without decomposition?
  4. Does a random-intercept model allow the time slope to vary across people?
  5. If the intercept–slope correlation is negative, why must the time zero point be stated first?
  6. In a logistic GLMM, does predict(..., re.form = NA) return the overall marginal risk?
  7. If a model has a singular fit, should every random slope be removed automatically?
  8. If a mixed model uses every available follow-up record, has attrition bias been eliminated?
  9. If students belong to both schools and residential communities and neither is nested in the other, what structure is appropriate?
  10. When can an adjusted mixed-model treatment coefficient receive a causal interpretation?
Show exercise answers
  1. No. An ICC describes model-scale variance or correlation, not the causal contribution of hospitals.
  2. Usually the 50-school design because it provides more independent school-level information.
  3. The within-person relationship comparing a person’s stress today with their own usual level, and the between-person relationship comparing chronically higher- and lower-stress people.
  4. No. A person-specific random time slope is required.
  5. The intercept represents time=0; changing that zero point changes the intercept and the interpretation of its correlation with the slope.
  6. No. It sets random effects to zero. A marginal risk requires integration over the random-effect distribution and a clearly defined covariate distribution.
  7. No. First assess within-cluster variation, scaling, sample support, the correlation structure, and the prespecified scientific question.
  8. No. Inference commonly still relies on a missing-at-random assumption given model information; missing-not-at-random mechanisms require sensitivity analysis.
  9. Use crossed random effects for school and community rather than forcing a nested structure.
  10. Only when the intervention, time zero, exchangeability, positivity, consistency, absence of interference, and relevant selection mechanisms support identification.

16 Quick reference

quick_reference <- data.frame(
  Goal = c(
    "Two-level random-intercept LMM", "Random slope", "Three-level longitudinal model",
    "Binary-outcome GLMM", "Extract fixed effects", "Variance components",
    "Conditional modes of random effects", "Singularity", "Typical-cluster prediction",
    "Existing-cluster conditional prediction"
  ),
  `R entry point` = c(
    "lme4::lmer(y ~ x + (1 | cluster), data=d)",
    "lme4::lmer(y ~ x + (1 + x | cluster), data=d)",
    "lme4::lmer(y ~ time*A + (1+time|person) + (1|site), data=d)",
    "lme4::glmer(y ~ x + (1|cluster), family=binomial, data=d)",
    "lme4::fixef(model)",
    "lme4::VarCorr(model)",
    "lme4::ranef(model, condVar=TRUE)",
    "lme4::isSingular(model)",
    "predict(model, re.form=NA)",
    "predict(model, re.form=NULL)"
  ),
  Reminder = c(
    "Report sample sizes at every level and the ICC",
    "Requires within-cluster variation in x and enough clusters",
    "IDs must be correct; assess residual serial correlation",
    "The effect is usually a conditional OR",
    "A fixed effect is not synonymous with an unconditional marginal effect",
    "Report SDs, correlations, and the residual scale",
    "These are not error-free cluster truths",
    "A boundary diagnostic, not an automatic term-deletion command",
    "In a GLMM this is not the integrated marginal mean",
    "Use only for clusters with estimated random effects"
  ),
  check.names = FALSE
)

knitr::kable(quick_reference, caption = "Core R entry points and interpretation reminders for multilevel models")
Core R entry points and interpretation reminders for multilevel models
Goal R entry point Reminder
Two-level random-intercept LMM lme4::lmer(y ~ x + (1 | cluster), data=d) Report sample sizes at every level and the ICC
Random slope lme4::lmer(y ~ x + (1 + x | cluster), data=d) Requires within-cluster variation in x and enough clusters
Three-level longitudinal model lme4::lmer(y ~ time*A + (1+time|person) + (1|site), data=d) IDs must be correct; assess residual serial correlation
Binary-outcome GLMM lme4::glmer(y ~ x + (1|cluster), family=binomial, data=d) The effect is usually a conditional OR
Extract fixed effects lme4::fixef(model) A fixed effect is not synonymous with an unconditional marginal effect
Variance components lme4::VarCorr(model) Report SDs, correlations, and the residual scale
Conditional modes of random effects lme4::ranef(model, condVar=TRUE) These are not error-free cluster truths
Singularity lme4::isSingular(model) A boundary diagnostic, not an automatic term-deletion command
Typical-cluster prediction predict(model, re.form=NA) In a GLMM this is not the integrated marginal mean
Existing-cluster conditional prediction predict(model, re.form=NULL) Use only for clusters with estimated random effects

Minimal analysis workflow

# 1. Audit units, IDs, cluster sizes, and within-cluster predictor variation
# at every level.
table(d$cluster_id)
tapply(d$x, d$cluster_id, stats::var)

# 2. Fit an empty model first and decompose the variance.
m0 <- lme4::lmer(y ~ 1 + (1 | cluster_id), data = d, REML = TRUE)
lme4::VarCorr(m0)

# 3. Decompose within and between components according to the question, and
# specify a scientifically supported random structure.
m1 <- lme4::lmer(
  y ~ x_within + x_cluster_mean + z_cluster +
    x_within:z_cluster + (1 + x_within | cluster_id),
  data = d,
  REML = TRUE,
  control = lme4::lmerControl(optimizer = "bobyqa")
)

# 4. Check fixed effects, variance components, singularity, residuals, and
# influential clusters.
lme4::fixef(m1)
lme4::VarCorr(m1)
lme4::isSingular(m1)

# 5. Distinguish conditional predictions for existing clusters, predictions
# for a typical cluster, and predictions for a new cluster.
predict(m1, re.form = NULL)
predict(m1, re.form = NA)

16.1 Final checklist

  • Are the meaning of each row, each level, and the target of inference explicit?
  • Are IDs globally unique, and are nested, crossed, or multiple-membership structures represented correctly?
  • Are sample sizes at every level, the range of cluster sizes, and imbalance reported?
  • Are predictors decomposed into within and between components according to the scientific question?
  • Do the zero points and centering conventions make intercepts and interactions interpretable?
  • Does the design provide within-cluster variation and enough clusters for each random slope?
  • Were ML or REML, the optimizer, interval method, and degrees-of-freedom method specified in advance?
  • Were convergence, singularity, residuals, serial correlation, and influential clusters checked?
  • For a GLMM, are conditional, typical-cluster, and marginal results distinguished?
  • Are random effects interpreted with shrinkage and uncertainty rather than used as direct ranks?
  • Were missingness, informative cluster size, measurement error, and selection assessed separately?
  • Is causal language aligned with the study design and identification assumptions?

Where to go next

Next topics include nonlinear time trajectories, spline random slopes, heteroskedasticity and AR(1) residuals, cross-classified and multiple-membership models, Bayesian multilevel models, joint longitudinal–survival models, multilevel mediation, cluster-randomized trials, small-sample degrees-of-freedom corrections, parametric bootstrap, survey-weighted multilevel models, external validation, and dynamic prediction.

sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Tahoe 26.5.1
##
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
##
## locale:
## [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
##
## time zone: America/Edmonton
## tzcode source: internal
##
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base
##
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-169     cli_3.6.6        knitr_1.51       rlang_1.3.0
##  [5] xfun_0.60        reformulas_0.4.4 jsonlite_2.0.0   minqa_1.2.8
##  [9] htmltools_0.5.9  lme4_2.0-6       sass_0.4.10      rmarkdown_2.31
## [13] grid_4.6.1       evaluate_1.0.5   jquerylib_0.1.4  MASS_7.3-65
## [17] fastmap_1.2.0    yaml_2.3.12      lifecycle_1.0.5  compiler_4.6.1
## [21] Rcpp_1.1.2       lattice_0.22-9   digest_0.6.39    nloptr_2.2.1
## [25] R6_2.6.1         Rdpack_2.6.6     splines_4.6.1    rbibutils_2.4.1
## [29] bslib_0.12.0     Matrix_1.7-5     tools_4.6.1      boot_1.3-32
## [33] cachem_1.1.0

For educational use: This module supports learning about the modeling, diagnosis, interpretation, and reporting of multilevel data. It does not replace a study-specific statistical analysis plan, data audit, or professional statistical consultation. All data and numerical results are simulated.