Reviewer: “Why Didn’t You Use a Multilevel Model?”
Choose an analysis that respects the unit of assignment and dependence, then show what changes in the inference.
Symptoms
- Participants share schools, clinics, therapists, or households, but your analysis treats every row as independent.
- A reviewer asks for a multilevel model, and you are unsure whether a small ICC excuses the original analysis.
What this usually means
The concern is dependence, not the name of the model. In a cluster-randomized study, treatment is assigned to clusters; additional participants within those clusters do not create additional independent treatment assignments. A mixed model is one response, but an appropriate cluster-level analysis or suitably corrected marginal analysis may also answer the question.
Common causes
- The data export has one participant per row, obscuring the assignment unit.
- The proposed effect is between clusters while the analysis uses participant-level degrees of freedom.
- Repeated measurements and nested or crossed groups have been treated as interchangeable.
Run these checks
- Map assignment, measurement, and sampling units before fitting anything.
- Count independent clusters in each arm and describe their sizes; distinguish cluster-level from within-cluster predictors.
- Check whether cluster size is related to outcomes and whether the target averages participants or clusters.
- Compare effect estimates and uncertainty under an analysis that respects the design.
What not to do
Do not decide from an arbitrary ICC cutoff, or add a random intercept solely to satisfy a reviewer. A nonsignificant variance estimate does not undo cluster randomization. Avoid relying on uncorrected asymptotic inference with very few clusters.
Treatment options
For balanced parallel cluster randomization with a continuous outcome, compare cluster means. For participant covariates or more complex dependence, consider mixed models, GEE, or cluster-robust inference with suitable small-sample corrections. Unequal cluster sizes, informative size, binary outcomes, and observational exposures need additional decisions about weighting and the estimand.
Worked example
This synthetic study has 20 clusters, 25 participants per cluster, and equal allocation. The script compares an independent-participant t test with a t test on the cluster means. Both estimate the intervention-minus-control difference as 0.264.
| Analysis | SE | 95% CI | p-value |
|---|---|---|---|
| Independent participants | 0.120 | 0.029 to 0.499 | 0.0281 |
| Cluster means; 18 df | 0.436 | −0.652 to 1.179 | 0.5525 |
The point estimate agrees because the clusters are equally sized. The uncertainty changes sharply. These simulated results demonstrate the assignment-unit issue; they do not establish that cluster means are optimal for every clustered dataset.
What to tell the reviewers
Adapt the following only after running the corresponding checks:
Treatment was assigned to clusters. We therefore reanalyzed the continuous outcome using cluster means with equal weight in this balanced design. The estimated difference was unchanged, but its confidence interval widened. We revised the inference and reported the cluster counts, sizes, and analysis unit.
See it in R and Python
Python uses the same synthetic participant data as R and reproduces both the independent-participant and cluster-mean t tests.
Python dependencies: NumPy and SciPy. Install with python -m pip install numpy scipy.
Download the shared synthetic CSV beside the Python script before running it. Identical seeds in R and Python do not generally generate identical samples.
# Synthetic balanced cluster-randomized example; base R only.
# Treatment is allocated to clusters. Cluster-mean t test respects that unit.
set.seed(4004)
k <- 20; m <- 25
cluster <- rep(seq_len(k), each=m)
arm <- rep(rep(0:1, each=k/2), each=m)
u <- rnorm(k, sd=1)
y <- 0.5*arm + u[cluster] + rnorm(k*m)
d <- data.frame(cluster, arm, y)
naive <- t.test(y ~ arm, data=d, var.equal=TRUE)
means <- aggregate(y ~ cluster + arm, d, mean)
valid <- t.test(y ~ arm, data=means, var.equal=TRUE)
cat(sprintf('Participants=%d; clusters=%d; cluster size=%d\n', nrow(d),k,m))
cat(sprintf('Difference=%.3f; naive SE=%.3f; cluster SE=%.3f\n', diff(tapply(y,arm,mean)), naive$stderr, valid$stderr))
cat(sprintf('Naive CI=[%.3f, %.3f]; cluster CI=[%.3f, %.3f]\n', -naive$conf.int[2],-naive$conf.int[1],-valid$conf.int[2],-valid$conf.int[1]))
cat(sprintf('Naive p=%.4f; cluster p=%.4f; cluster df=%.0f\n',naive$p.value,valid$p.value,valid$parameter))
stopifnot(valid$stderr > naive$stderr, nrow(means)==k)
# Maintainer-only export; run from repository root. Synthetic data only.
if (Sys.getenv("CLINIC_EXPORT_DATA") == "1") {
write.csv(d,'analysis-clinic/examples/case-004-shared-data.csv',row.names=FALSE,quote=FALSE)
}
# DASS Analysis Clinic Case 004; Python dependencies: numpy, scipy.
from pathlib import Path
import numpy as np
from scipy import stats, optimize
HERE = Path(__file__).resolve().parent
def ols(X, y):
b = np.linalg.lstsq(X, y, rcond=None)[0]
residual = y - X @ b
df = len(y) - X.shape[1]
cov = (residual @ residual / df) * np.linalg.inv(X.T @ X)
se = np.sqrt(np.diag(cov))
ci = np.column_stack((b-stats.t.ppf(.975, df)*se,b+stats.t.ppf(.975,df)*se))
return b, cov, ci
def power(n, d):
# Matches R power.t.test(strict=FALSE): rejection tail in effect direction.
return stats.nct.sf(stats.t.ppf(.975,2*n-2),2*n-2,d*np.sqrt(n/2))
d=np.genfromtxt(HERE/'case-004-shared-data.csv',delimiter=',',names=True)
def compare(a,b):
effect=b.mean()-a.mean(); df=len(a)+len(b)-2
pooled=((len(a)-1)*a.var(ddof=1)+(len(b)-1)*b.var(ddof=1))/df
se=np.sqrt(pooled*(1/len(a)+1/len(b)))
return effect,se,effect+np.array([-1,1])*stats.t.ppf(.975,df)*se,2*stats.t.sf(abs(effect/se),df)
naive=compare(d['y'][d['arm']==0],d['y'][d['arm']==1])
clusters=np.unique(d['cluster']); means=np.array([d['y'][d['cluster']==k].mean() for k in clusters]); arms=np.array([d['arm'][d['cluster']==k][0] for k in clusters])
valid=compare(means[arms==0],means[arms==1])
print('Independent participants:',naive); print('Cluster means:',valid)
assert np.isclose(naive[0],valid[0]) and valid[1]>naive[1] and np.isclose(valid[1],.436,atol=.0005)
Download Python script · Download shared synthetic CSV
Every number above comes from one base-R script, with no packages to install.
Download case-004-multilevel-model.R →