Clustered and repeated-measures data

The assumption the output cannot show you

Regression diagnostics lists four assumptions, and three of them have a plot: residuals against fitted values, the QQ plot, influence. The fourth — observations are independent of one another — is where that page stops, with the remark that this one cannot be seen. Poisson regression twice says “use GEE or a mixed-effects model” and the site has never had that page. This page settles all three debts.

Independence cannot be seen because it is not a property of the data, it is a property of how the data came about. In the same CSV, a row that is one patient and a row that is the same patient’s third visit look identical. You have to know it from the study design; no model will tell you.

Violations are more common in clinical research than people expect:

• Repeated measures: the same patient scored for pain at 30 minutes, 90 minutes, 4 hours and the next morning
• Cluster randomisation: clinics or wards are randomised, and patients within a clinic share the same physicians, the same processes, the same case mix
• Multicentre studies: patients from one hospital resemble each other in quality of care more than they resemble patients from another
• Matched designs: left and right eye, twins, matched cases and controls
• Multistage sampling: sample schools, then sample students within schools

The common thread is that the data has two levels: the patient (or clinic), and the measurements inside the patient. Every complication follows from one fact — observations within a cluster are more alike than observations from different clusters.

The example on this page

medicaldata::licorice_gargle is a randomised trial of a licorice gargle given before intubation to prevent post-operative sore throat. 235 patients were randomised, 117 to control and 118 to licorice, and sore throat was scored at four time points.

Reshaping wide to long gives 940 rows, of which 8 scores are missing.

There are two things to read off this plot. First, a patient resembles himself — that is what clustering means. But this is not something to judge by eye, and this dataset happens to show why. Most of what fills the frame is the crowd sitting near 0: of 233 patients, 122 scored exactly the same value all four times, and 120 of those scored 0 every time, so of course their lines are flat. The high-scoring lines do swing violently, but there are only 8 patients whose within-patient range reaches 4 points. The impression that “the lines are mostly flat” is carried by the patients scoring zero; it is not that the data lacks movement. What actually supports the claim that a patient resembles himself is the intraclass correlation (ICC) of 0.540 computed later on this page, not the visual. A trajectory plot shows you the shape of the variation; quantifying it takes variance components.

Second, the individual lines are too crowded to reveal the group difference, which is why the mean trajectories are drawn on top — but the means are computed from all of the data, not from the lines you can pick out. The individual lines tell you how large the variation is; the mean trajectories tell you where the effect is. A plot with only the means invites the reader to assume everyone looks like the average.

Running it yourself

library(medicaldata)
library(nlme)

lg <- licorice_gargle
lg$id <- seq_len(nrow(lg))

tp <- c("pacu30min_throatPain", "pacu90min_throatPain",
        "postOp4hour_throatPain", "pod1am_throatPain")

# Wide to long: one row per patient-by-time-point
long <- do.call(rbind, lapply(seq_along(tp), function(i)
  data.frame(id = lg$id, treat = lg$treat, t = i, pain = lg[[tp[i]]])))
L <- long[!is.na(long$pain), ]
L$idf <- factor(L$id)

# Look at the missingness pattern BEFORE analysing
table(rowSums(is.na(lg[, tp])))

# Naive model: 932 rows treated as 932 independent observations
summary(lm(pain ~ treat, data = L))

# Mixed model: one random intercept per patient
m <- lme(pain ~ treat, random = ~ 1 | idf, data = L, method = "REML")
summary(m)
VarCorr(m)                      # variance components; the ICC comes from here

# A third estimator: average each patient, then do a two-sample comparison
pm <- aggregate(pain ~ id + treat, data = L, FUN = mean)
summary(lm(pain ~ treat, data = pm))

# Within-subject contrasts: time
summary(lm(pain ~ factor(t) * treat, data = L))                    # naive
summary(lme(pain ~ factor(t) * treat, random = ~ 1 | idf, data = L))

# Cluster-robust (sandwich) standard errors, written out:
#   (X'X)^-1 [ sum_g X_g' u_g u_g' X_g ] (X'X)^-1
# one cluster per patient, times the CR1 small-sample factor.
cluster_robust_se <- function(fit, cluster) {
  X  <- model.matrix(fit)
  u  <- as.numeric(residuals(fit))
  cl <- factor(cluster)
  G <- nlevels(cl); N <- nrow(X); K <- ncol(X)
  bread <- solve(crossprod(X))
  meat  <- matrix(0, K, K)
  for (g in levels(cl)) {
    idx <- which(cl == g)
    sg  <- crossprod(X[idx, , drop = FALSE], u[idx])   # this cluster's score vector
    meat <- meat + tcrossprod(sg)
  }
  adj <- (G / (G - 1)) * ((N - 1) / (N - K))           # CR1
  setNames(sqrt(diag(bread %*% meat %*% bread * adj)), colnames(X))
}

cluster_robust_se(lm(pain ~ treat, data = L), L$id)
cluster_robust_se(lm(pain ~ factor(t) * treat, data = L), L$id)
# In practice: sandwich::vcovCL(fit, cluster = L$id, type = "HC1").
# This site installs no extra packages, hence the hand-written version.

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

lg = sm.datasets.get_rdataset("licorice_gargle", "medicaldata").data
lg["id"] = np.arange(len(lg))

tp = ["pacu30min_throatPain", "pacu90min_throatPain",
      "postOp4hour_throatPain", "pod1am_throatPain"]

L = lg.melt(id_vars=["id", "treat"], value_vars=tp,
            var_name="tp", value_name="pain").dropna(subset=["pain"])
L["t"] = L["tp"].map({v: i + 1 for i, v in enumerate(tp)})

print(smf.ols("pain ~ treat", data=L).fit().summary())        # naive

md = smf.mixedlm("pain ~ treat", data=L, groups=L["id"])      # random intercept
print(md.fit(reml=True).summary())

pm = L.groupby(["id", "treat"], as_index=False)["pain"].mean()
print(smf.ols("pain ~ treat", data=pm).fit().summary())

# Cluster-robust: OLS with a sandwich variance, clustering on patient
print(smf.ols("pain ~ treat", data=L)
        .fit(cov_type="cluster", cov_kwds={"groups": L["id"]}).summary())

Direction one: the BETWEEN-subject standard error comes out too small

Start with the simplest question. How much lower is the sore-throat score in the licorice arm?

The naive approach throws all 932 rows into lm(). The point estimate is -0.582 points (lower under licorice), with a standard error of 0.0686, t = -8.48, p < 0.001, on 930 residual degrees of freedom.

The mixed model gives exactly the same point estimate, but a standard error of 0.1113, t = -5.23, p < 0.001, and the degrees of freedom for treat drop to 231.

The naive standard error is 38.3% smaller than the correct one.

Direction two: the WITHIN-subject standard errors come out too large

Most textbooks stop after the previous section, which leaves a rule that is easy to misremember: “ignoring clustering makes standard errors too small.” That sentence is only half true, and the half it omits is the more useful one.

Refit as pain ~ factor(time) * treat and look at the contrasts between time points:

The last column says what happened. The ratio is 1 for the intercept and for treat — the mixed model and the naive model return identical standard errors there, because in this interaction model treat is the between-subject contrast at the PACU 30 min time point alone, where each patient contributes exactly one row and clustering has nothing to act on. For every other row the ratio is 0.673 — the mixed-model standard error is 32.7% smaller than the naive one.

What the mixed model is doing

lme(pain ~ treat, random = ~ 1 | patient) adds one term: each patient gets his own intercept $b_i$, and those $b_i$ are assumed to be drawn from a normal distribution with mean zero. The model becomes

$b_i$ is a random effect: it is not a parameter to be estimated but a distribution to be estimated. That is what separates it from treating the patient id as a fixed-effect factor with 233 levels — the latter would burn 232 degrees of freedom, and would swallow the between-subject effect entirely, since treat is constant within a patient and therefore perfectly collinear with the patient id.

That standard error is just “average each patient, then compare”

Collapse each patient’s four scores to a single mean, which leaves 233 rows, and run the most ordinary two-sample comparison there is. The standard error is 0.1113 — equal to the mixed model’s 0.1113 to nine decimal places. The two differ by 2.4e-10, which is the REML optimiser’s convergence tolerance, not a difference between the estimators.

ICC: how much of the variation is “people differ”

The mixed model splits the variation in two. The between-patient variance is 0.5947 and the within-patient variance is 0.5068. The intraclass correlation coefficient (ICC) is the first as a share of the total:

which comes to 0.540.

There is a qualifier here that is easy to drop. Those two variances come from lme(pain ~ treat, random = ~1 | patient), a model that already contains treat, so what they partition is the variation that remains once the group difference is taken out, not the total variation in the raw scores. The correct reading is: within a single arm, the differences between patients are slightly larger than one patient’s fluctuation over time (0.540 against 0.460). The qualifier carries weight here, because 0.540 is only just above one half and a different model specification could put it on the other side. Always report which model, and which variance decomposition, an ICC came from.

The same formula has a second use, and it is the one most papers mean by “ICC = 0.85”. The ICC on this page asks how alike two people in the same cluster are — it measures structure. The reliability ICC asks how alike two measurements of the same subject are — it measures the instrument or the rater. The variance decomposition is identical; the question is not, and the reliability side splits further into six forms (single vs average measurement, consistency vs absolute agreement, raters treated as random vs fixed), none of which can be interpreted without being named. See measurement reliability and the ICC.

An ICC can equally be read as the correlation between any two measurements on the same patient; under a random-intercept model those are the same statement.

Design effect: turning an ICC into “how much sample was wasted”

Sample size and power notes that cluster-randomised trials must be inflated by a design effect $1 + (m-1)\times \mathrm{ICC}$, where $m$ is the number of individuals per cluster. On this page the “cluster” is the patient, and $m$ is the number of measurements per patient, 4.

That gives a design effect of 2.620, so the 932 rows are worth roughly 356 independent observations for the between-subject comparison.

Read that number in both directions, or it looks like pure loss. Downward, 932 rows shrink to 356, which is the price of treating repeated measurements as independent. Upward, 356 is more than the 233 patients, and that surplus is exactly the point made in the previous section: repeated measurement did buy something, it just has a ceiling.

Mixed models, GEE and cluster-robust: three different questions

There are three mainstream routes for handling clustering — a mixed model, GEE (generalised estimating equations) and cluster-robust standard errors. They are not three algorithms for the same thing; they answer three different questions. This is the section most often skipped and most often misused.

The conditional-versus-marginal distinction, in one question each:

• The mixed model’s coefficient answers “if this same person had been in the other arm, how different would his score be?” — it is defined holding that patient’s random intercept $b_i$ fixed.
• The GEE and cluster-robust coefficients answer “if the whole population switched to licorice, how much would the average score change?”

In the linear model on this page the two have the same point estimate (-0.582), so the distinction is invisible. That is a special property of linear models. Once the effect measure is non-collapsible — the odds ratio above all — conditional and marginal effects differ numerically, and the conditional OR is always further from 1; this is the non-collapsibility discussed on logistic regression. Not every non-linear link behaves this way: the log link is collapsible, so a Poisson mixed model and a Poisson GEE estimate the same rate ratio and differ only in the intercept. On one dataset a mixed logistic model and a GEE will report different-sized ORs; both are correct, and they answer different questions.

Choosing between them:

• To say “what would happen to this patient on treatment” (clinical decisions, individual prediction) → mixed model
• To say “what would happen if this policy went population-wide” (public health, policy evaluation) → GEE or cluster-robust
• When the clustering is a nuisance rather than the object of interest (centre effects in a multicentre trial, say) → cluster-robust is the least work, because it does not require you to guess the correlation structure correctly

Cluster-robust is not “multiply by a fixed factor”

In the table above, the cluster-robust column is worth a look on its own. For the between-subject effect the cluster-robust standard error 0.1115 and the mixed model’s 0.1113 very nearly coincide. But the three time-point contrasts have cluster-robust standard errors of 0.064, 0.133 and 0.127 — very different from one another, and the first (PACU 90 min versus PACU 30 min) is even smaller than the mixed model’s 0.092.

The reason is that a sandwich estimator imposes no correlation structure at all. The random-intercept model assumes any two measurements on a patient are equally correlated (exchangeable), so every within-subject contrast shares the one correction factor 0.673. The sandwich estimates each contrast’s actual variance from the residuals, so contrasts between adjacent, strongly correlated time points come out more precise and distant ones less so.

Common misuses

Related pages

• Ordinal logistic regression — the same licorice_gargle data with a different kind of outcome (graded cough)
• Regression diagnostics — how to check the three assumptions other than independence
• Sample size and power — how the design effect enters a sample-size calculation
• Logistic regression — non-collapsibility, that is, why conditional and marginal effects part company under a non-collapsible link such as the logit

Reproducing every number on this page

/opt/homebrew/bin/Rscript figures/scripts/B2-07-clustered.R

Report a problem on this page