Advanced Statistics and Machine Learning for Health Research

Mark Khurana; Neil Scheidwasser

8 Mixed-Effects Models for Clustered and Longitudinal Data

Almost every clinical dataset hides a quiet problem: the rows are not independent of one another. You measure the same patient’s blood pressure at five visits; you recruit patients from twelve different ICUs; a multi-site trial pools records from six hospitals. Standard regression assumes every observation is a fresh, unrelated draw — as if each blood-pressure reading came from a different, unrelated person. When that assumption is false, the analysis quietly fools itself into thinking it has more independent information than it really does, and it hands you standard errors that are too small and p-values that are too impressive. This chapter introduces mixed-effects models, the standard tool for getting clustered and repeated data right, and shows how to read their output in plain clinical terms.

8.1 Introduction

Two situations recur constantly in health research, and both break the “independent observations” assumption that ordinary regression relies on:

Repeated measurements on the same patient. You record HbA1c every three months, or blood pressure at every clinic visit. Two readings from the same patient resemble each other far more than two readings from different patients — a patient who runs high tends to run high every time.
Patients grouped (nested) within a cluster. Patients treated in the same ICU, registered with the same GP practice, or recruited at the same trial site share an environment, staff, and case mix. Patients within a cluster resemble each other more than patients from different clusters.

In both cases we say the data are clustered (or hierarchical, or multilevel): observations come in groups, and observations within a group are correlated.

Why ignoring clustering is dangerous

If you fit an ordinary regression that pretends every row is independent, you overstate how much information you have. The result is standard errors that are too small, confidence intervals that are too narrow, and p-values that are too small — so you declare effects “significant” that are not. A study of 50 patients each measured 10 times has 500 rows, but it does not contain 500 patients’ worth of independent information. Mixed models correct for exactly this.

flowchart TD
    H1["Hospital A"] --> P1["Patient 1"]
    H1 --> P2["Patient 2"]
    H2["Hospital B"] --> P3["Patient 3"]
    P1 --> V1["Visit 1, 2, 3 ..."]
    P2 --> V2["Visit 1, 2, 3 ..."]
    P3 --> V3["Visit 1, 2, 3 ..."]

Figure 8.1: Clustered data: patients are nested within hospitals, and repeated visits are nested within patients. Observations inside the same box resemble each other.

8.2 Fixed effects and random effects

The phrase “mixed-effects” means the model mixes two kinds of effect. Understanding the difference is the heart of this chapter.

A fixed effect is an effect you care about estimating directly and want to generalise to everyone: the average change in blood pressure per year of age, the average difference between a drug and placebo. These are the familiar regression coefficients. There is one number for the whole population, and that number is the answer to your clinical question.

A random effect is different. Instead of estimating a separate, free coefficient for every cluster, we assume the clusters’ values are drawn from a distribution — typically a normal distribution centred on the overall average. We do not care about the individual value for “Hospital B” as such; we treat the hospitals as a random sample from the larger population of hospitals, and we estimate the spread between them. The random effects are how the model acknowledges that clusters differ without spending a coefficient on each one.

A useful way to hear it: fixed effects answer “what is the effect?”; random effects answer “how much do the clusters vary around that effect?”

8.2.1 Random intercepts: each cluster its own baseline

The simplest mixed model gives each cluster its own intercept. Imagine systolic blood pressure measured across several GP practices. Practice A’s patients might average 138 mmHg, Practice B’s 129, Practice C’s 144 — different baselines, perhaps because of differing populations or measurement habits. A random intercept lets every practice slide its own regression line up or down to match its baseline, while the slope (the effect of, say, age or treatment) stays shared across practices.

Picture a set of parallel lines, one per cluster: same tilt, different heights.

8.2.2 Random slopes: each cluster its own trend

Sometimes clusters differ not just in baseline but in how strongly a predictor acts. Perhaps a drug lowers blood pressure more in some hospitals than others, or blood pressure climbs faster with age in one practice than another. A random slope lets the slope of a predictor vary by cluster too. Now the lines differ in both height and tilt.

You add random slopes when you believe the effect of a variable genuinely differs between clusters — and you want your uncertainty to reflect that.

flowchart LR
    A["Random intercept<br/>= each cluster its<br/>own baseline<br/>(parallel lines)"]
    B["Random slope<br/>= each cluster its<br/>own trend<br/>(lines fan out)"]

Figure 8.2: Random intercepts give each cluster its own baseline (parallel lines); random slopes also let each cluster’s trend differ (lines fan out).

8.2.3 Partial pooling: borrowing strength

Here is the clever part. You could analyse each cluster completely separately (no pooling) — but a practice with only three patients would give a wildly unreliable estimate. Or you could lump everyone together and ignore clusters (complete pooling) — the mistake we are trying to avoid. Mixed models do something in between, called partial pooling: each cluster’s estimate is pulled (or shrunk) towards the overall average, and the pull is stronger for small or noisy clusters and weaker for large, well-measured ones.

The clinical intuition: a practice that contributes only a handful of patients should not be trusted to define its own baseline from scratch; we sensibly assume it is probably not too far from the average until its own data prove otherwise. Partial pooling does this automatically, which is why mixed models give stabler, more honest cluster estimates than fitting each cluster alone.

The same structure turns up again in Chapter 24, where the clusters are entire published studies rather than clinics, and the “random-effects model” of meta-analysis is this model with the within-cluster variances already known.

8.2.4 The intraclass correlation (ICC)

Once we let clusters differ, a natural question is: how much of the variation is between clusters versus within them? The intraclass correlation coefficient (ICC) answers exactly this. It is the proportion of the total variation that is between clusters:

\[ \text{ICC} = \frac{\sigma^2_{\text{between}}}{\sigma^2_{\text{between}} + \sigma^2_{\text{within}}} \]

where \(\sigma^2_{\text{between}}\) is the variance of the cluster-level baselines and \(\sigma^2_{\text{within}}\) is the leftover (residual) variation within a cluster. Read it as a number from 0 to 1:

ICC near 0: clusters are essentially interchangeable; clustering hardly matters.
ICC near 1: almost all the variation is between clusters; which cluster you are in matters enormously.

An ICC of 0.20 for blood pressure across GP practices means 20% of the variation in readings is explained by which practice a patient attends. The ICC is also a tidy summary of “how correlated are two observations from the same cluster” — which is the very thing ordinary regression ignores. For a conceptual, clinician-friendly walkthrough of how the ICC and multilevel models quantify this “contextual” clustering, the brief tutorial series by Merlo et al. (2005) works entirely through a worked blood-pressure example.

8.2.5 Continuous, binary, or count outcomes: LMM vs GLMM

When the outcome is continuous (blood pressure, HbA1c, weight) we fit a linear mixed model (LMM) — a linear regression with random effects bolted on. When the outcome is binary (died/survived, readmitted/not), a count (number of admissions), or otherwise non-normal, we fit a generalised linear mixed model (GLMM) — the mixed-model version of logistic or Poisson regression. The random-effects idea is identical; only the link between the linear predictor and the outcome changes, exactly as in ordinary generalised linear models.

8.2.6 A note on GEE: marginal vs conditional models

There is a sibling approach worth knowing about: generalised estimating equations (GEE). Both GEE and mixed models handle clustered data, but they answer subtly different questions.

Mixed models are conditional (subject-specific). A coefficient describes the effect for a given patient or cluster — “for this patient, treatment lowers the odds of the event by this much.”
GEE is marginal (population-average). A coefficient describes the effect averaged over the whole population — “across the population, the proportion with the event shifts by this much.”

For continuous outcomes with an identity link the two coincide, but for logistic models they genuinely differ, and the population-average effect is usually closer to the null. Use a mixed model when you care about individual or cluster-level trajectories and the variance structure itself; use GEE when you want a population-average effect and mainly regard the clustering as a nuisance to be corrected for. GEE is fitted with geepack in R. For a clinician-oriented tutorial that sets the summary-statistic, repeated-measures ANOVA, GEE and mixed-model approaches side by side — with worked examples from published studies — see Schober and Vetter (2018).

8.3 Packages

For R, lme4 (with its lmer function for LMMs and glmer for GLMMs) is the long-standing standard, paired with lmerTest to obtain p-values for the fixed effects (base lme4 deliberately omits them). glmmTMB is a newer, fast alternative that handles trickier cases (zero-inflation, complex variance structures) with the same formula syntax, and geepack provides GEE. For Python, statsmodels offers MixedLM for linear mixed models, taking data in long format (one row per observation, not per patient).

8.4 Worked example: blood pressure over repeated visits across clinics

We simulate a longitudinal study: patients seen at several clinics, each with systolic blood pressure (SBP) measured at multiple visits over time. We expect (a) clinics to differ in baseline SBP, (b) SBP to drift with visit number, and (c) patients to differ from one another. We will fit a random-intercept model, then a random-slope model, and interpret the variance components.

Code

library(tidyverse) # data wrangling and simulation
library(lme4) # lmer() for linear mixed models
library(lmerTest) # adds p-values to lmer() output

set.seed(42)

# --- Simulate a longitudinal, clustered dataset ---
n_clinic <- 8 # 8 clinics
n_patient <- 25 # 25 patients per clinic
n_visit <- 5 # 5 visits per patient

# Each clinic has its own baseline SBP (random intercept, SD = 6)
clinic_effect <- rnorm(n_clinic, 0, 6)
# Each patient has their own baseline on top of the clinic (SD = 8)
patient_effect <- rnorm(n_clinic * n_patient, 0, 8)

bp <- expand_grid(
  clinic = 1:n_clinic,
  patient = 1:n_patient,
  visit = 1:n_visit
) |>
  mutate(
    patient_id = interaction(clinic, patient), # unique patient label
    pid_index = as.integer(patient_id),
    # SBP = grand mean + clinic baseline + patient baseline
    #       + a small downward drift per visit + noise
    sbp = 135 +
      clinic_effect[clinic] +
      patient_effect[pid_index] +
      -1.5 * (visit - 1) +
      rnorm(n(), 0, 5)
  )

# --- Random-intercept model ---
# (1 | clinic/patient_id) = random intercept for clinic,
#   and for patient nested within clinic
m_intercept <- lmer(
  sbp ~ visit + (1 | clinic / patient_id),
  data = bp
)
summary(m_intercept)

What the code shows. We build a dataset of 8 clinics, 25 patients each, 5 visits each (1000 rows). Each clinic gets its own baseline shift, each patient another shift on top, and SBP drifts down about 1.5 mmHg per visit. The model sbp ~ visit + (1 | clinic/patient_id) has one fixed effect (visit, the shared trend we care about) and random intercepts at two levels: clinics, and patients nested within clinics. In the summary() output, look first at the Random effects table: it reports a variance (and standard deviation) for clinic, for patient_id within clinic, and a Residual. These tell you how much SBP varies between clinics, between patients, and visit-to-visit within a patient. Then the Fixed effects table gives the visit coefficient — expect roughly \(-1.5\), meaning SBP falls about 1.5 mmHg per visit on average, with a standard error and (thanks to lmerTest) a p-value. The clinical reading: after properly accounting for the fact that visits cluster within patients and patients within clinics, blood pressure declines modestly over follow-up.

Code

library(tidyverse)
library(lme4)

set.seed(42)
# (Same simulation as above, repeated so this chunk stands alone)
n_clinic <- 8
n_patient <- 25
n_visit <- 5
clinic_effect <- rnorm(n_clinic, 0, 6)
patient_effect <- rnorm(n_clinic * n_patient, 0, 8)
bp <- expand_grid(
  clinic = 1:n_clinic,
  patient = 1:n_patient,
  visit = 1:n_visit
) |>
  mutate(
    patient_id = interaction(clinic, patient),
    pid_index = as.integer(patient_id),
    sbp = 135 +
      clinic_effect[clinic] +
      patient_effect[pid_index] +
      -1.5 * (visit - 1) +
      rnorm(n(), 0, 5)
  )

m_intercept <- lmer(sbp ~ visit + (1 | clinic / patient_id), data = bp)

# --- Extract variance components and compute the ICC by hand ---
vc <- as.data.frame(VarCorr(m_intercept))
var_clinic <- vc$vcov[vc$grp == "clinic"]
var_patient <- vc$vcov[vc$grp == "patient_id:clinic"]
var_resid <- vc$vcov[vc$grp == "Residual"]

icc_clinic <- var_clinic / (var_clinic + var_patient + var_resid)
cat("ICC (clinic level):", round(icc_clinic, 3), "\n")

What the code shows. VarCorr() pulls out the three estimated variances. We then compute the clinic-level ICC — the clinic variance divided by the total variance — by hand, so the formula is transparent. The printed value (around 0.25–0.30 given our simulated variances: a clinic SD of 6 against a patient SD of 8 and residual SD of 5, so roughly \(36/(36+64+25) \approx 0.29\)) is read as: “this fraction of the total variation in blood pressure is attributable to differences between clinics.” Even a modest ICC like this is enough to badly distort an ordinary regression’s standard errors, which is precisely why we used a mixed model.

Code

library(tidyverse)
library(lme4)
library(lmerTest)

set.seed(42)
n_clinic <- 8
n_patient <- 25
n_visit <- 5
clinic_effect <- rnorm(n_clinic, 0, 6)
patient_effect <- rnorm(n_clinic * n_patient, 0, 8)
bp <- expand_grid(
  clinic = 1:n_clinic,
  patient = 1:n_patient,
  visit = 1:n_visit
) |>
  mutate(
    patient_id = interaction(clinic, patient),
    pid_index = as.integer(patient_id),
    sbp = 135 +
      clinic_effect[clinic] +
      patient_effect[pid_index] +
      -1.5 * (visit - 1) +
      rnorm(n(), 0, 5)
  )

# --- Random-intercept model (refit here so this chunk stands alone) ---
# (1 | clinic/patient_id) expands to (1|clinic) + (1|clinic:patient_id)
m_intercept <- lmer(sbp ~ visit + (1 | clinic / patient_id), data = bp)

# --- Random-slope model: let the visit trend vary by clinic ---
# Same nesting, but clinic now also gets its own visit slope.
m_slope <- lmer(
  sbp ~ visit + (visit | clinic) + (1 | clinic:patient_id),
  data = bp
)
summary(m_slope)

# Compare the two models with a likelihood-ratio test
anova(m_intercept, m_slope)

What the code shows. The term (visit | clinic) now gives each clinic both its own intercept and its own visit slope, so the blood-pressure trend is allowed to differ from clinic to clinic. In the summary(), the random-effects table gains a variance for the visit slope and a correlation between intercept and slope. The anova() call performs a likelihood-ratio test comparing the simpler random-intercept model with this richer random-slope model: a small p-value would say the trends genuinely differ between clinics and the slopes are worth keeping. Because our data were simulated with a single shared slope, expect the test to be non-significant here — a reassuring reminder that you should add random slopes only when the data (and clinical sense) support them, since needless complexity can make models fail to converge.

Code

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

rng = np.random.default_rng(42)

# --- Simulate the same longitudinal, clustered dataset (long format) ---
n_clinic, n_patient, n_visit = 8, 25, 5
clinic_effect = rng.normal(0, 6, n_clinic)

rows = []
pid = 0
for c in range(n_clinic):
    for p in range(n_patient):
        patient_effect = rng.normal(0, 8)   # each patient's own baseline
        for v in range(n_visit):
            sbp = (135 + clinic_effect[c] + patient_effect
                   - 1.5 * v + rng.normal(0, 5))
            rows.append({"clinic": c, "patient_id": pid, "visit": v, "sbp": sbp})
        pid += 1

bp = pd.DataFrame(rows)

# --- Random-intercept model: random intercept for patient ---
# statsmodels MixedLM groups by ONE variable; we group by patient.
m_intercept = smf.mixedlm("sbp ~ visit", data=bp, groups=bp["patient_id"])
res_intercept = m_intercept.fit()
print(res_intercept.summary())

# ICC at the patient level = group variance / (group variance + residual)
var_group = float(res_intercept.cov_re.iloc[0, 0])
var_resid = float(res_intercept.scale)
icc = var_group / (var_group + var_resid)
print(f"\nPatient-level ICC: {icc:.3f}")

What the code shows. We build the same data in long format (one row per visit) because that is what statsmodels expects. mixedlm("sbp ~ visit", groups=bp["patient_id"]) fits a linear mixed model with visit as the fixed effect and a random intercept per patient. In the printed summary, the visit coefficient should sit near \(-1.5\) (SBP falling per visit), and the Group Var line reports the between-patient variance. We then compute the patient-level ICC — between-patient variance over total variance — which expresses how much of the variation is due to patients being persistently different from one another. Note one practical limitation: statsmodels MixedLM groups by a single level at a time, so to model both clinics and patients you would typically choose the level you care most about or use the vc_formula argument for additional variance components; R’s lme4 handles nested levels more naturally with (1 | clinic/patient_id).

Code

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

rng = np.random.default_rng(42)
n_clinic, n_patient, n_visit = 8, 25, 5
clinic_effect = rng.normal(0, 6, n_clinic)
rows, pid = [], 0
for c in range(n_clinic):
    for p in range(n_patient):
        patient_effect = rng.normal(0, 8)
        for v in range(n_visit):
            sbp = 135 + clinic_effect[c] + patient_effect - 1.5 * v + rng.normal(0, 5)
            rows.append({"clinic": c, "patient_id": pid, "visit": v, "sbp": sbp})
        pid += 1
bp = pd.DataFrame(rows)

# --- Random-slope model: each patient gets their own visit slope ---
m_slope = smf.mixedlm("sbp ~ visit", data=bp,
                      groups=bp["patient_id"], re_formula="~visit")
res_slope = m_slope.fit()
print(res_slope.summary())

What the code shows. Adding re_formula="~visit" tells statsmodels to give each patient their own slope on visit as well as their own intercept — the Python equivalent of a random-slope model. The summary now reports an extra variance for the visit random effect (how much patients’ trajectories differ in steepness) and a covariance between each patient’s intercept and slope. Interpret it clinically as: not only do patients start at different blood-pressure levels, their pressures may also change at different rates over follow-up. As in R, only keep random slopes if they are justified; superfluous random slopes often cause the optimiser to warn that the model did not converge.

Putting the variance components into words

The hardest part of a mixed model for newcomers is reading the random-effects table. Translate each variance into a clinical sentence: the clinic variance is “how much clinics differ in baseline,” the patient variance is “how much patients differ from each other,” and the residual is “how much a single patient’s readings bounce around visit to visit.” If the between-cluster variance is large relative to the residual, clustering matters a lot — and analysing the data as if it were independent would have been seriously misleading.

“Model failed to converge”

Mixed models, especially with random slopes or many variance components, sometimes report a convergence warning. This usually means the model is too complex for the data to support — there is not enough information to estimate every variance. The remedy is usually to simplify: drop a random slope, or reduce the random-effects structure, rather than to ignore the warning.

8.5 Exercises

Exercise 1: Spot the clustering (conceptual)

For each study below, state (a) what the clusters are, (b) whether the data are also longitudinal, and (c) what would go wrong if you analysed it with ordinary regression ignoring the clustering.

A trial of a new inhaler recruiting asthma patients from 15 GP practices.
A study measuring fasting glucose monthly for a year in 200 diabetics.
A multi-site ICU study of a sepsis bundle across 6 hospitals, one outcome per patient.

Solution

Code

# =============================================================================
# Chapter 6b, Exercise 1: Spot the clustering (conceptual)
# Identify clusters, longitudinal structure, and the risk of ignoring both.
# =============================================================================

# This is a conceptual exercise. The answer is written as structured comments.
# For each study, state (a) the clusters, (b) whether the data are also
# longitudinal, and (c) what goes wrong if you use ordinary regression that
# ignores the clustering.

# -----------------------------------------------------------------------------
# (a) A trial of a new inhaler recruiting asthma patients from 15 GP practices.
# -----------------------------------------------------------------------------
#   Clusters:     The 15 GP practices. Patients are nested within practices;
#                 patients from the same practice share staff, prescribing
#                 habits, catchment population and case mix, so their outcomes
#                 resemble one another more than outcomes from other practices.
#   Longitudinal? Not as described. Each patient contributes (presumably) a
#                 single outcome, so the data are CLUSTERED but not longitudinal.
#                 (It would become longitudinal if each patient were measured
#                 repeatedly over time.)
#   Ignoring it:  Ordinary regression treats all patients as independent and so
#                 overstates the amount of information. Because patients within
#                 a practice are positively correlated (ICC > 0), the standard
#                 error on the inhaler effect is too small, the confidence
#                 interval too narrow and the p-value too small -- you risk
#                 declaring the inhaler effective when the evidence does not
#                 support it (inflated type-I error). The fix is a random
#                 intercept for practice: sbp/outcome ~ treatment + (1 | practice).

# -----------------------------------------------------------------------------
# (b) A study measuring fasting glucose monthly for a year in 200 diabetics.
# -----------------------------------------------------------------------------
#   Clusters:     The 200 patients. The repeated monthly readings are nested
#                 within each patient; a patient who runs high tends to run high
#                 at every visit.
#   Longitudinal? YES. Each patient has 12 monthly measurements ordered in
#                 time (200 x 12 = 2400 rows), so this is a classic repeated-
#                 measures / longitudinal design.
#   Ignoring it:  You have 2400 rows but NOT 2400 independent observations --
#                 you effectively have ~200 patients' worth of independent
#                 information. Ordinary regression would give standard errors
#                 that are far too small and p-values that are far too
#                 impressive, and it could not separate within-patient change
#                 over time from between-patient differences. A mixed model with
#                 a random intercept (and possibly a random slope for time) per
#                 patient handles this correctly: glucose ~ month + (month | id).

# -----------------------------------------------------------------------------
# (c) A multi-site ICU study of a sepsis bundle across 6 hospitals,
#     one outcome per patient.
# -----------------------------------------------------------------------------
#   Clusters:     The 6 hospitals. Patients are nested within hospital and share
#                 protocols, staffing and case mix.
#   Longitudinal? NO. There is a single outcome per patient, so the data are
#                 cross-sectional but clustered.
#   Ignoring it:  Patients within the same hospital are correlated, so ordinary
#                 regression again gives standard errors that are too small and
#                 overstates the significance of the sepsis-bundle effect.
#                 A caution specific to this study: with only 6 clusters there
#                 are very few hospitals from which to estimate the between-
#                 hospital variance, so a random-effects estimate of that
#                 variance will itself be imprecise (some analysts prefer to
#                 treat so few sites as fixed effects, or to use a cluster-robust
#                 approach). The clustering must still be accounted for either
#                 way -- ignoring it is not an option.

cat("Exercise 1 is conceptual -- see the commented answers above.\n")

Code

# =============================================================================
# Chapter 6b, Exercise 1: Spot the clustering (conceptual)
# Identify clusters, longitudinal structure, and the risk of ignoring both.
# =============================================================================

# This is a conceptual exercise. The answer is written as structured comments.
# For each study, state (a) the clusters, (b) whether the data are also
# longitudinal, and (c) what goes wrong if you use ordinary regression that
# ignores the clustering.

# -----------------------------------------------------------------------------
# (a) A trial of a new inhaler recruiting asthma patients from 15 GP practices.
# -----------------------------------------------------------------------------
#   Clusters:     The 15 GP practices. Patients are nested within practices;
#                 patients from the same practice share staff, prescribing
#                 habits, catchment population and case mix, so their outcomes
#                 resemble one another more than outcomes from other practices.
#   Longitudinal? Not as described. Each patient contributes (presumably) a
#                 single outcome, so the data are CLUSTERED but not longitudinal.
#                 (It would become longitudinal if each patient were measured
#                 repeatedly over time.)
#   Ignoring it:  Ordinary regression treats all patients as independent and so
#                 overstates the amount of information. Because patients within
#                 a practice are positively correlated (ICC > 0), the standard
#                 error on the inhaler effect is too small, the confidence
#                 interval too narrow and the p-value too small -- you risk
#                 declaring the inhaler effective when the evidence does not
#                 support it (inflated type-I error). The fix is a random
#                 intercept for practice: outcome ~ treatment + (1 | practice).

# -----------------------------------------------------------------------------
# (b) A study measuring fasting glucose monthly for a year in 200 diabetics.
# -----------------------------------------------------------------------------
#   Clusters:     The 200 patients. The repeated monthly readings are nested
#                 within each patient; a patient who runs high tends to run high
#                 at every visit.
#   Longitudinal? YES. Each patient has 12 monthly measurements ordered in
#                 time (200 x 12 = 2400 rows), so this is a classic repeated-
#                 measures / longitudinal design.
#   Ignoring it:  You have 2400 rows but NOT 2400 independent observations --
#                 you effectively have ~200 patients' worth of independent
#                 information. Ordinary regression would give standard errors
#                 that are far too small and p-values that are far too
#                 impressive, and it could not separate within-patient change
#                 over time from between-patient differences. A mixed model with
#                 a random intercept (and possibly a random slope for time) per
#                 patient handles this correctly.

# -----------------------------------------------------------------------------
# (c) A multi-site ICU study of a sepsis bundle across 6 hospitals,
#     one outcome per patient.
# -----------------------------------------------------------------------------
#   Clusters:     The 6 hospitals. Patients are nested within hospital and share
#                 protocols, staffing and case mix.
#   Longitudinal? NO. There is a single outcome per patient, so the data are
#                 cross-sectional but clustered.
#   Ignoring it:  Patients within the same hospital are correlated, so ordinary
#                 regression again gives standard errors that are too small and
#                 overstates the significance of the sepsis-bundle effect.
#                 A caution specific to this study: with only 6 clusters there
#                 are very few hospitals from which to estimate the between-
#                 hospital variance, so a random-effects estimate of that
#                 variance will itself be imprecise (some analysts prefer to
#                 treat so few sites as fixed effects, or to use a cluster-robust
#                 approach). The clustering must still be accounted for either
#                 way -- ignoring it is not an option.

print("Exercise 1 is conceptual -- see the commented answers above.")

Exercise 2: Fit and interpret a random-intercept model

Using the simulated blood-pressure dataset from this chapter, fit a random-intercept model with visit and a simulated binary treatment variable as fixed effects.

Report the treatment fixed-effect coefficient and interpret it clinically.
Compute the ICC and put it into one plain-English sentence.
Refit the same data with an ordinary linear regression (lm) that ignores clustering. Compare the standard error on treatment. Which is larger, and why does that matter?

Solution

Code

# =============================================================================
# Chapter 6b, Exercise 2: Fit and interpret a random-intercept model
# Recreate the BP data, add a binary treatment, fit a random-intercept LMM,
# compute the ICC, and compare the treatment SE against ordinary lm().
# =============================================================================

library(tidyverse)   # data wrangling and simulation
library(lme4)        # lmer() for linear mixed models
# lmerTest adds p-values to lmer(); use it if available, otherwise fall back
# to estimates + profile confidence intervals from base lme4.
has_lmerTest <- requireNamespace("lmerTest", quietly = TRUE)
if (has_lmerTest) suppressMessages(library(lmerTest))

set.seed(42)

# --- Recreate the chapter's longitudinal, clustered BP dataset --------------
n_clinic  <- 8     # 8 clinics
n_patient <- 25    # 25 patients per clinic
n_visit   <- 5     # 5 visits per patient

clinic_effect  <- rnorm(n_clinic, 0, 6)              # clinic baselines (SD 6)
patient_effect <- rnorm(n_clinic * n_patient, 0, 8)  # patient baselines (SD 8)

# Add a simulated binary treatment, assigned at the PATIENT level (a patient is
# either treated or not, and keeps that assignment across all their visits).
# True treatment effect built into the simulation: -5 mmHg.
treatment_by_patient <- rbinom(n_clinic * n_patient, 1, 0.5)

bp <- expand_grid(
  clinic  = 1:n_clinic,
  patient = 1:n_patient,
  visit   = 1:n_visit
) |>
  mutate(
    patient_id = interaction(clinic, patient),   # unique patient label
    pid_index  = as.integer(patient_id),
    treatment  = treatment_by_patient[pid_index],
    sbp = 135 +
          clinic_effect[clinic] +
          patient_effect[pid_index] +
          -1.5 * (visit - 1) +      # downward drift per visit
          -5   * treatment +        # true treatment effect: -5 mmHg
          rnorm(n(), 0, 5)          # residual noise
  )

cat("Dataset:", nrow(bp), "rows |",
    n_clinic, "clinics x", n_patient, "patients x", n_visit, "visits\n\n")

# --- Fit the random-intercept model -----------------------------------------
# Fixed effects: visit + treatment. Random intercepts for clinic and for
# patient nested within clinic.
m_ri <- lmer(sbp ~ visit + treatment + (1 | clinic/patient_id), data = bp)

cat("=== Random-intercept model summary ===\n")
print(summary(m_ri))

# --- (a) Treatment fixed effect and its interpretation ----------------------
coefs   <- summary(m_ri)$coefficients
beta_trt <- coefs["treatment", "Estimate"]
se_trt   <- coefs["treatment", "Std. Error"]

cat("\n=== (a) Treatment fixed effect ===\n")
cat(sprintf("Estimate: %.3f mmHg   SE: %.3f\n", beta_trt, se_trt))
if (has_lmerTest && "Pr(>|t|)" %in% colnames(coefs)) {
  cat(sprintf("p-value (lmerTest): %.4g\n", coefs["treatment", "Pr(>|t|)"]))
} else {
  ci <- confint(m_ri, parm = "treatment", method = "Wald")
  cat(sprintf("95%% CI (Wald): %.3f to %.3f\n", ci[1], ci[2]))
}
cat("Interpretation: after accounting for visit and the clustering of visits\n")
cat("  within patients and patients within clinics, being on treatment is\n")
cat(sprintf(
  "  associated with a %.2f mmHg change in systolic BP (a reduction),\n",
  beta_trt))
cat("  holding visit number constant.\n")

# --- (b) Intraclass correlation ---------------------------------------------
vc <- as.data.frame(VarCorr(m_ri))
var_clinic  <- vc$vcov[vc$grp == "clinic"]
var_patient <- vc$vcov[vc$grp == "patient_id:clinic"]
var_resid   <- vc$vcov[vc$grp == "Residual"]
var_total   <- var_clinic + var_patient + var_resid

icc_clinic  <- var_clinic / var_total
icc_patient <- var_patient / var_total
# Proportion of variance at or above the patient level (clinic + patient):
icc_cluster <- (var_clinic + var_patient) / var_total

cat("\n=== (b) Variance components and ICC ===\n")
cat(sprintf("Between-clinic variance : %.2f\n", var_clinic))
cat(sprintf("Between-patient variance: %.2f\n", var_patient))
cat(sprintf("Residual variance       : %.2f\n", var_resid))
cat(sprintf("Clinic-level ICC : %.3f\n", icc_clinic))
cat(sprintf("Patient-level ICC: %.3f\n", icc_patient))
cat(sprintf(
  "Plain English: about %.0f%% of the total variation in blood pressure lies\n",
  100 * icc_cluster))
cat(sprintf(
  "  between clinics and patients rather than within a patient's own visits\n"))
cat(sprintf(
  "  (clinic alone accounts for roughly %.0f%%).\n", 100 * icc_clinic))

# --- (c) Compare with ordinary lm() that ignores clustering -----------------
m_lm <- lm(sbp ~ visit + treatment, data = bp)
lm_coefs <- summary(m_lm)$coefficients
se_trt_lm <- lm_coefs["treatment", "Std. Error"]

cat("\n=== (c) Ordinary lm() ignoring clustering ===\n")
cat(sprintf("Treatment estimate (lm)  : %.3f\n", lm_coefs["treatment", "Estimate"]))
cat(sprintf("Treatment SE (lm)        : %.3f\n", se_trt_lm))
cat(sprintf("Treatment SE (mixed)     : %.3f\n", se_trt))
cat(sprintf("Ratio mixed/lm           : %.2f\n", se_trt / se_trt_lm))
cat("\nWhich is larger, and why it matters:\n")
cat("  The MIXED-model SE is larger. Treatment is a BETWEEN-PATIENT variable\n")
cat("  (constant across a patient's 5 visits), so the real amount of\n")
cat("  independent information about treatment is roughly the number of\n")
cat("  patients, not the 1000 rows. Ordinary lm() pretends all 1000 rows are\n")
cat("  independent, so it reports an SE that is too small -- an over-optimistic\n")
cat("  confidence interval and p-value. Ignoring the clustering would make the\n")
cat("  treatment effect look more certain than the data can support.\n")

Code

# =============================================================================
# Chapter 6b, Exercise 2: Fit and interpret a random-intercept model
# Recreate the BP data, add a binary treatment, fit a random-intercept MixedLM,
# compute the ICC, and compare the treatment SE against ordinary OLS.
# =============================================================================

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

np.random.seed(42)

# --- Recreate the chapter's longitudinal, clustered BP dataset (long format) -
n_clinic, n_patient, n_visit = 8, 25, 5
clinic_effect = np.random.normal(0, 6, n_clinic)   # clinic baselines (SD 6)

rows = []
pid = 0
for c in range(n_clinic):
    for p in range(n_patient):
        patient_effect = np.random.normal(0, 8)        # patient baseline (SD 8)
        # Treatment assigned at the PATIENT level (constant across the patient's
        # visits). True built-in effect: -5 mmHg.
        treatment = np.random.binomial(1, 0.5)
        for v in range(n_visit):
            sbp = (135 + clinic_effect[c] + patient_effect
                   - 1.5 * v            # downward drift per visit
                   - 5.0 * treatment    # true treatment effect: -5 mmHg
                   + np.random.normal(0, 5))
            rows.append({"clinic": c, "patient_id": pid, "visit": v,
                         "treatment": treatment, "sbp": sbp})
        pid += 1

bp = pd.DataFrame(rows)
print(f"Dataset: {len(bp)} rows | {n_clinic} clinics x {n_patient} patients "
      f"x {n_visit} visits\n")

# --- Fit the random-intercept model -----------------------------------------
# statsmodels MixedLM groups by ONE level; we use a random intercept per
# patient (the level at which treatment varies), which absorbs the repeated-
# visit correlation. Fixed effects: visit + treatment.
m_ri = smf.mixedlm("sbp ~ visit + treatment", data=bp, groups=bp["patient_id"])
res_ri = m_ri.fit()

print("=== Random-intercept model summary ===")
print(res_ri.summary())

# --- (a) Treatment fixed effect and its interpretation ----------------------
beta_trt = res_ri.params["treatment"]
se_trt = res_ri.bse["treatment"]
p_trt = res_ri.pvalues["treatment"]

print("\n=== (a) Treatment fixed effect ===")
print(f"Estimate: {beta_trt:.3f} mmHg   SE: {se_trt:.3f}   p-value: {p_trt:.4g}")
print("Interpretation: after accounting for visit and the fact that repeated")
print("  visits are clustered within patients, being on treatment is")
print(f"  associated with a {beta_trt:.2f} mmHg change in systolic BP (a")
print("  reduction), holding visit number constant.")

# --- (b) Intraclass correlation ---------------------------------------------
# ICC at the patient level = between-patient variance / total variance.
var_group = float(res_ri.cov_re.iloc[0, 0])   # between-patient variance
var_resid = float(res_ri.scale)               # residual (within-patient) var
icc = var_group / (var_group + var_resid)

print("\n=== (b) Variance components and ICC ===")
print(f"Between-patient variance: {var_group:.2f}")
print(f"Residual variance       : {var_resid:.2f}")
print(f"Patient-level ICC       : {icc:.3f}")
print(f"Plain English: about {100*icc:.0f}% of the total variation in blood")
print("  pressure is due to persistent differences between patients rather than")
print("  visit-to-visit fluctuation within a patient.")

# --- (c) Compare with ordinary OLS that ignores clustering ------------------
m_ols = smf.ols("sbp ~ visit + treatment", data=bp).fit()
se_trt_ols = m_ols.bse["treatment"]

print("\n=== (c) Ordinary OLS ignoring clustering ===")
print(f"Treatment estimate (OLS): {m_ols.params['treatment']:.3f}")
print(f"Treatment SE (OLS)      : {se_trt_ols:.3f}")
print(f"Treatment SE (mixed)    : {se_trt:.3f}")
print(f"Ratio mixed/OLS         : {se_trt / se_trt_ols:.2f}")
print("\nWhich is larger, and why it matters:")
print("  The MIXED-model SE is larger. Treatment is a BETWEEN-PATIENT variable")
print("  (constant across a patient's 5 visits), so the real amount of")
print("  independent information about treatment is roughly the number of")
print("  patients, not the 1000 rows. Ordinary OLS pretends all 1000 rows are")
print("  independent, so it reports an SE that is too small -- an over-")
print("  optimistic confidence interval and p-value. Ignoring the clustering")
print("  would make the treatment effect look more certain than the data")
print("  can support.")

Exercise 3: Random intercept vs random slope

Extend the model in Exercise 2 to allow the visit trend to vary by clinic (a random slope).

Use a likelihood-ratio test (anova in R) to compare the random-intercept and random-slope models.
Does allowing clinic-specific trends improve the model? What does your answer imply clinically?

Solution

Code

# =============================================================================
# Chapter 6b, Exercise 3: Random intercept vs random slope
# Extend the Exercise-2 model with a clinic-specific visit slope and compare
# the random-intercept and random-slope models with a likelihood-ratio test.
# =============================================================================

library(tidyverse)   # data wrangling and simulation
library(lme4)        # lmer() and anova() likelihood-ratio test
has_lmerTest <- requireNamespace("lmerTest", quietly = TRUE)
if (has_lmerTest) suppressMessages(library(lmerTest))

set.seed(42)

# --- Recreate the chapter's BP dataset with a binary treatment --------------
n_clinic  <- 8
n_patient <- 25
n_visit   <- 5

clinic_effect  <- rnorm(n_clinic, 0, 6)
patient_effect <- rnorm(n_clinic * n_patient, 0, 8)
treatment_by_patient <- rbinom(n_clinic * n_patient, 1, 0.5)

bp <- expand_grid(
  clinic  = 1:n_clinic,
  patient = 1:n_patient,
  visit   = 1:n_visit
) |>
  mutate(
    patient_id = interaction(clinic, patient),
    pid_index  = as.integer(patient_id),
    treatment  = treatment_by_patient[pid_index],
    # NOTE: the data are simulated with a SINGLE shared visit slope (-1.5),
    # i.e. no genuine clinic-to-clinic variation in the trend.
    sbp = 135 + clinic_effect[clinic] + patient_effect[pid_index] +
          -1.5 * (visit - 1) + -5 * treatment + rnorm(n(), 0, 5)
  )

# --- Random-intercept model (from Exercise 2) -------------------------------
m_ri <- lmer(sbp ~ visit + treatment + (1 | clinic/patient_id), data = bp)

# --- Random-slope model: let the visit trend vary by clinic -----------------
# (visit | clinic) gives each clinic its own intercept AND its own visit slope.
m_rs <- lmer(sbp ~ visit + treatment + (visit | clinic) + (1 | clinic:patient_id),
             data = bp)

cat("=== Random-slope model summary ===\n")
print(summary(m_rs))

# --- (a) Likelihood-ratio test comparing the two models ---------------------
# anova() on lmer objects refits both models with ML and performs the LRT.
cat("\n=== (a) Likelihood-ratio test (random intercept vs random slope) ===\n")
lrt <- anova(m_ri, m_rs)
print(lrt)

p_lrt <- lrt$`Pr(>Chisq)`[2]

# --- (b) Interpretation ------------------------------------------------------
cat("\n=== (b) Does a clinic-specific trend improve the model? ===\n")
cat(sprintf("LRT p-value: %.4g\n", p_lrt))
if (p_lrt < 0.05) {
  cat("The random slope significantly improves fit: the visit trend genuinely\n")
  cat("differs between clinics.\n")
} else {
  cat("The random slope does NOT significantly improve fit (p > 0.05).\n")
  cat("Clinically: there is no evidence that blood pressure changes at\n")
  cat("different rates across clinics -- the single shared downward trend is\n")
  cat("adequate. This is expected, because the data were simulated with one\n")
  cat("common visit slope. The lesson: add random slopes only when the data\n")
  cat("(and clinical sense) support them; needless complexity can also trigger\n")
  cat("convergence warnings.\n")
}

Code

# =============================================================================
# Chapter 6b, Exercise 3: Random intercept vs random slope
# Extend the Exercise-2 model with a clinic-specific visit slope and compare
# the random-intercept and random-slope models with a likelihood-ratio test.
# =============================================================================

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from scipy import stats

np.random.seed(42)

# --- Recreate the chapter's BP dataset with a binary treatment --------------
n_clinic, n_patient, n_visit = 8, 25, 5
clinic_effect = np.random.normal(0, 6, n_clinic)

rows = []
pid = 0
for c in range(n_clinic):
    for p in range(n_patient):
        patient_effect = np.random.normal(0, 8)
        treatment = np.random.binomial(1, 0.5)
        for v in range(n_visit):
            # NOTE: simulated with a SINGLE shared visit slope (-1.5), i.e. no
            # genuine clinic-to-clinic variation in the trend.
            sbp = (135 + clinic_effect[c] + patient_effect
                   - 1.5 * v - 5.0 * treatment + np.random.normal(0, 5))
            rows.append({"clinic": c, "patient_id": pid, "visit": v,
                         "treatment": treatment, "sbp": sbp})
        pid += 1

bp = pd.DataFrame(rows)

# statsmodels MixedLM groups by ONE level, so to let the visit trend vary "by
# clinic" we group by clinic here. (R's lme4 can nest clinic/patient in one
# call; in statsmodels we choose the level of interest for this comparison.)

# --- Random-intercept model: random intercept per clinic --------------------
m_ri = smf.mixedlm("sbp ~ visit + treatment", data=bp,
                   groups=bp["clinic"]).fit(reml=True)

# --- Random-slope model: each clinic gets its own visit slope ----------------
m_rs = smf.mixedlm("sbp ~ visit + treatment", data=bp,
                   groups=bp["clinic"], re_formula="~visit").fit(reml=True)

print("=== Random-slope model summary ===")
print(m_rs.summary())

# --- (a) Likelihood-ratio test comparing the two models ---------------------
# The random-slope model adds a slope variance and an intercept-slope
# covariance -> 2 extra parameters. Both models share the same fixed effects,
# so a REML-based LRT on the random structure is valid.
# (Clamp at 0: when the extra random-effect variance collapses to the boundary
# the richer model can land on a marginally lower REML log-likelihood, giving a
# tiny negative statistic that simply means "no improvement".)
lr_stat = max(0.0, 2 * (m_rs.llf - m_ri.llf))
df_diff = 2
p_lrt = stats.chi2.sf(lr_stat, df_diff)

print("\n=== (a) Likelihood-ratio test (random intercept vs random slope) ===")
print(f"logLik (random intercept): {m_ri.llf:.2f}")
print(f"logLik (random slope)    : {m_rs.llf:.2f}")
print(f"LR statistic: {lr_stat:.3f} on {df_diff} df")
print(f"LRT p-value : {p_lrt:.4g}")

# --- (b) Interpretation ------------------------------------------------------
print("\n=== (b) Does a clinic-specific trend improve the model? ===")
if p_lrt < 0.05:
    print("The random slope significantly improves fit: the visit trend")
    print("genuinely differs between clinics.")
else:
    print("The random slope does NOT significantly improve fit (p > 0.05).")
    print("Clinically: there is no evidence that blood pressure changes at")
    print("different rates across clinics -- the single shared downward trend")
    print("is adequate. This is expected, because the data were simulated with")
    print("one common visit slope. The lesson: add random slopes only when the")
    print("data (and clinical sense) support them; needless complexity can also")
    print("trigger convergence warnings.")

Exercise 4: GLMM and the marginal/conditional distinction (conceptual)

You have a binary outcome (30-day readmission) for patients clustered within hospitals.

Would you fit an LMM or a GLMM? Why?
A colleague fits the same question with GEE and gets a smaller coefficient. Explain, in plain language, why the population-average (GEE) and subject-specific (mixed model) estimates differ for a logistic model.
Which would you report if a hospital manager asks “what happens to this hospital’s readmission rate if we adopt the intervention”?

Solution

Code

# =============================================================================
# Chapter 6b, Exercise 4: GLMM and the marginal/conditional distinction
# Conceptual: LMM vs GLMM, why GEE and mixed-model estimates differ, and which
# to report for a hospital-specific question.
# =============================================================================

# This is a conceptual exercise. The answer is written as structured comments.
# Setting: a binary outcome (30-day readmission) for patients clustered within
# hospitals.

# -----------------------------------------------------------------------------
# (a) LMM or GLMM? Why?
# -----------------------------------------------------------------------------
#   Fit a GLMM (a generalised linear mixed model -- here a logistic mixed
#   model). The outcome is BINARY (readmitted / not), not continuous and
#   normally distributed, so a linear mixed model (LMM) is inappropriate: an
#   LMM assumes a continuous, roughly normal outcome with constant variance and
#   can predict probabilities outside 0-1. A GLMM keeps the random-effects idea
#   (e.g. a random intercept for hospital) but uses the logit link and a
#   binomial error, exactly as ordinary logistic regression extends linear
#   regression. In R: glmer(readmit ~ intervention + (1 | hospital),
#   family = binomial).

# -----------------------------------------------------------------------------
# (b) Why do the population-average (GEE) and subject-specific (mixed) estimates
#     differ for a logistic model, and why is the GEE coefficient smaller?
# -----------------------------------------------------------------------------
#   The two answer different questions:
#     * The mixed model is CONDITIONAL / SUBJECT-SPECIFIC. Its coefficient is
#       the effect of the intervention FOR A GIVEN HOSPITAL (holding that
#       hospital's random intercept fixed).
#     * GEE is MARGINAL / POPULATION-AVERAGE. Its coefficient is the effect
#       AVERAGED OVER the whole population of hospitals.
#   For a linear (identity-link) model the two coincide. For a logistic model
#   they genuinely differ because the logit link is NON-LINEAR: the average of
#   the individual (hospital-specific) log-odds effects is not equal to the
#   effect on the population-averaged probability. Averaging a curved (S-shaped)
#   relationship over the spread of hospital random effects flattens it, pulling
#   the population-average effect TOWARD THE NULL. So the marginal (GEE)
#   coefficient is smaller in magnitude than the conditional (mixed-model) one.
#   (This is the well-known non-collapsibility of the odds ratio.) The larger
#   the between-hospital variance, the bigger the gap.
#   Plain language: "the average of the individual effects is not the same as
#   the effect on the average patient", and the logit link's curvature shrinks
#   the averaged effect toward zero.

# -----------------------------------------------------------------------------
# (c) Which would you report if a hospital manager asks: "what happens to THIS
#     hospital's readmission rate if we adopt the intervention?"
# -----------------------------------------------------------------------------
#   Report the SUBJECT-SPECIFIC / CONDITIONAL estimate from the MIXED MODEL
#   (GLMM). The manager is asking a cluster-specific question -- the change for
#   one particular hospital -- which is exactly what the conditional coefficient
#   describes. The marginal (GEE) estimate answers a different question: the
#   average shift in readmission across the whole population of hospitals, which
#   is what you would report to a health-system regulator planning a
#   population-wide rollout. General good practice when reporting a mixed model:
#   state the fixed-effect estimates with confidence intervals (on the odds-
#   ratio scale for a logistic GLMM), report the random-effects variance /
#   between-hospital SD (or the ICC) to convey how much hospitals differ, and be
#   explicit that the estimates are conditional (subject-specific).

cat("Exercise 4 is conceptual -- see the commented answers above.\n")

Code

# =============================================================================
# Chapter 6b, Exercise 4: GLMM and the marginal/conditional distinction
# Conceptual: LMM vs GLMM, why GEE and mixed-model estimates differ, and which
# to report for a hospital-specific question.
# =============================================================================

# This is a conceptual exercise. The answer is written as structured comments.
# Setting: a binary outcome (30-day readmission) for patients clustered within
# hospitals.

# -----------------------------------------------------------------------------
# (a) LMM or GLMM? Why?
# -----------------------------------------------------------------------------
#   Fit a GLMM (a generalised linear mixed model -- here a logistic mixed
#   model). The outcome is BINARY (readmitted / not), not continuous and
#   normally distributed, so a linear mixed model (LMM) is inappropriate: an
#   LMM assumes a continuous, roughly normal outcome with constant variance and
#   can predict probabilities outside 0-1. A GLMM keeps the random-effects idea
#   (e.g. a random intercept for hospital) but uses the logit link and a
#   binomial error, exactly as ordinary logistic regression extends linear
#   regression. In statsmodels this is BinomialBayesMixedGLM or the R route
#   glmer(readmit ~ intervention + (1 | hospital), family = binomial).

# -----------------------------------------------------------------------------
# (b) Why do the population-average (GEE) and subject-specific (mixed) estimates
#     differ for a logistic model, and why is the GEE coefficient smaller?
# -----------------------------------------------------------------------------
#   The two answer different questions:
#     * The mixed model is CONDITIONAL / SUBJECT-SPECIFIC. Its coefficient is
#       the effect of the intervention FOR A GIVEN HOSPITAL (holding that
#       hospital's random intercept fixed).
#     * GEE is MARGINAL / POPULATION-AVERAGE. Its coefficient is the effect
#       AVERAGED OVER the whole population of hospitals.
#   For a linear (identity-link) model the two coincide. For a logistic model
#   they genuinely differ because the logit link is NON-LINEAR: the average of
#   the individual (hospital-specific) log-odds effects is not equal to the
#   effect on the population-averaged probability. Averaging a curved (S-shaped)
#   relationship over the spread of hospital random effects flattens it, pulling
#   the population-average effect TOWARD THE NULL. So the marginal (GEE)
#   coefficient is smaller in magnitude than the conditional (mixed-model) one.
#   (This is the well-known non-collapsibility of the odds ratio.) The larger
#   the between-hospital variance, the bigger the gap.
#   Plain language: "the average of the individual effects is not the same as
#   the effect on the average patient", and the logit link's curvature shrinks
#   the averaged effect toward zero.

# -----------------------------------------------------------------------------
# (c) Which would you report if a hospital manager asks: "what happens to THIS
#     hospital's readmission rate if we adopt the intervention?"
# -----------------------------------------------------------------------------
#   Report the SUBJECT-SPECIFIC / CONDITIONAL estimate from the MIXED MODEL
#   (GLMM). The manager is asking a cluster-specific question -- the change for
#   one particular hospital -- which is exactly what the conditional coefficient
#   describes. The marginal (GEE) estimate answers a different question: the
#   average shift in readmission across the whole population of hospitals, which
#   is what you would report to a health-system regulator planning a
#   population-wide rollout. General good practice when reporting a mixed model:
#   state the fixed-effect estimates with confidence intervals (on the odds-
#   ratio scale for a logistic GLMM), report the random-effects variance /
#   between-hospital SD (or the ICC) to convey how much hospitals differ, and be
#   explicit that the estimates are conditional (subject-specific).

print("Exercise 4 is conceptual -- see the commented answers above.")

8.6 Summary

Clustered and longitudinal data — repeated measurements on patients, patients nested in hospitals or practices, multi-site studies — violate the independence assumption of ordinary regression. Ignoring the clustering produces standard errors and p-values that are too optimistic, leading to false claims of significance. Mixed-effects models solve this by combining fixed effects (the population-level effects you want to estimate) with random effects (which capture how clusters vary). Random intercepts give each cluster its own baseline; random slopes also let each cluster’s trend differ. Through partial pooling, mixed models give stable, honest cluster estimates, and the intraclass correlation summarises how much variation lives between clusters. Linear mixed models handle continuous outcomes and generalised linear mixed models handle binary or count outcomes. GEE offers a population-average alternative when the clustering is a nuisance rather than the focus.

Key Takeaways

If observations cluster (repeated measures, or patients within sites), do not use ordinary regression — it gives standard errors that are too small.
Fixed effects answer “what is the effect?”; random effects answer “how much do clusters vary around it?”
Random intercepts = each cluster its own baseline; random slopes = each cluster its own trend.
The ICC is the share of total variation that is between clusters — a quick gauge of how much clustering matters.
Use a GLMM (not an LMM) for binary or count outcomes.
Mixed models are subject-specific (conditional); GEE is population-average (marginal). For logistic models the two differ, so know which question you are answering.
Keep the random-effects structure as simple as the data support — over-complex models fail to converge.

8.7 References and Further Reading

For mixed-effects models and longitudinal data analysis, see Twisk (2013), Gelman and Hill (2007), Schober and Vetter (2018), and Merlo et al. (2005).
For GEE methods, see Hubbard et al. (2010).
For the R packages used in this chapter, see Bates et al. (2015) (lme4), Kuznetsova et al. (2017) (lmerTest), and Brooks et al. (2017) (glmmTMB).

Bates, Douglas, Martin Mächler, Ben Bolker, and Steve Walker. 2015. “Fitting Linear Mixed-Effects Models Using lme4.” Journal of Statistical Software 67 (1): 1–48. https://doi.org/10.18637/jss.v067.i01.

Brooks, Mollie E, Kasper Kristensen, Koen J van Benthem, et al. 2017. “glmmTMB Balances Speed and Flexibility Among Packages for Zero-Inflated Generalized Linear Mixed Modeling.” The R Journal 9 (2): 378–400.

Gelman, Andrew, and Jennifer Hill. 2007. Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press.

Hubbard, Alan E, Jennifer Ahern, Nancy L Fleischer, et al. 2010. “To GEE or Not to GEE: Comparing Population Average and Mixed Models for Estimating the Associations Between Neighborhood Risk Factors and Health.” Epidemiology 21 (4): 467–74. https://doi.org/10.1097/EDE.0b013e3181caeb90.

Kuznetsova, Alexandra, Per B Brockhoff, and Rune H B Christensen. 2017. “lmerTest Package: Tests in Linear Mixed Effects Models.” Journal of Statistical Software 82 (13): 1–26. https://doi.org/10.18637/jss.v082.i13.

Merlo, Juan, Basile Chaix, Min Yang, John Lynch, and Lennart Råstam. 2005. “A Brief Conceptual Tutorial of Multilevel Analysis in Social Epidemiology: Linking the Statistical Concept of Clustering to the Idea of Contextual Phenomenon.” Journal of Epidemiology and Community Health 59 (6): 443–49. https://doi.org/10.1136/jech.2004.023473.

Schober, Patrick, and Thomas R Vetter. 2018. “Repeated Measures Designs and Analysis of Longitudinal Data: If at First You Do Not Succeed—Try, Try Again.” Anesthesia & Analgesia 127 (2): 569–75. https://doi.org/10.1213/ANE.0000000000003511.

Twisk, Jos W R. 2013. Applied Longitudinal Data Analysis for Epidemiology: A Practical Guide. 2nd ed. Cambridge University Press.