The independence assumption: clustered data

Most statistical tests assume that observations are independent: once the model has done its work, one measurement’s error says nothing about another’s. Repeated measurements on the same unit, and units grouped in clusters, break that assumption, and the main damage is not to the estimates but to the standard errors, so to every interval and p-value built on them.

This post explains what independence means, where it fails, and how to judge it. It then works one real example through a deterministic engine: tree growth data analysed twice, once as if every measurement were independent and once with a random intercept for each tree. The point estimates match; the intervals do not.

What independence means

Independence is a statement about errors, not about raw values. Two measurements can both be large because both trees received the same treatment; the model explains that. Independence fails when two measurements share something the model does not include, so their errors move together.

The usual sources are few and easy to name:

• Repeated measures: the same animal, plant or participant measured several times.
• Clusters: pupils within schools, cells within animals, samples within batches, plots within fields.
• Time and space: measurements close in time or in location that share conditions.

The assumption applies within each group and between groups. In a two-sample t-test, the people in one arm must not influence each other, and nobody may sit in both arms unless you use a paired test. One-way ANOVA asks the same of every group. In regression it is stated about the residuals: the errors around the fitted line should not move together.

Hurlbert called the analysis that ignores this structure pseudoreplication: testing for a treatment effect with replicates that are not statistically independent. Lazic lists the same three sources. Using one issue of Nature Neuroscience as a case study, he found pseudoreplication in 12% of the papers and suspected it in a further 36%, where too little information was reported to be sure.

Why it matters: fewer units than rows

A cluster of five measurements on one tree does not carry five trees’ worth of information. If trees differ a lot from each other and little from one measurement to the next, the five rows mostly repeat what the first one told you. The effective sample size for a question about trees is closer to the number of trees than to the number of rows.

An analysis that counts rows treats that repetition as new evidence. Its standard errors shrink as if the sample were larger than it is. Aarts and colleagues showed that ignoring nesting can push the false-positive rate for a t-test far above the nominal 5%, and that improving such a design nearly always means collecting more independent units, not more observations per unit.

Worked example: Sitka spruce, with and without the tree

R’s MASS package includes the Sitka data: log size of 79 Sitka spruce trees, 54 grown in ozone-enriched chambers and 25 controls, each measured five times in 1988 at roughly monthly intervals. That is 395 rows from 79 trees.

The question is whether ozone exposure changes log size, allowing for growth over the season. Two predictors play different roles:

• Ozone treatment varies between trees. Each tree is either ozone or control for the whole season.
• Time varies within trees. Every tree is measured on every day.

The engine fitted a linear mixed model with a random intercept per tree: size ~ treat + Time + (1 | tree). For comparison, a second engine run fitted the same fixed effects with ordinary least squares, ignoring the tree. Both runs use the same data and seed.

	Ordinary regression (ignores tree)	Mixed model (random intercept per tree)
Ozone effect on log size	-0.211	-0.211
Its standard error	0.0686	0.149
95% interval	-0.346 to -0.0763	-0.507 to 0.085
p-value	0.00223	0.159
Time: standard error of the daily slope	0.00085	0.00027

The estimated ozone effect is the same in both: on average, ozone trees were about 0.211 smaller on the log scale. What changes is the uncertainty. Ignoring the tree, the standard error is under half the mixed model’s, the interval excludes zero and the p-value is small. Modelling the tree, the standard error is 0.149, the interval runs from -0.507 to 0.085 and the Kenward-Roger F test gives F = 2.02 on 1 and 77 degrees of freedom, p = 0.159.

The ordinary regression is not a second opinion here; it is the wrong model for these data, and its interval is too narrow. The mixed model’s result is inconclusive: the data are compatible with no ozone effect and with a reduction of about half a unit of log size. It does not show that ozone has no effect.

The degrees of freedom tell the same story. The test for ozone has 77 residual degrees of freedom: the 79 trees minus the two group means, not the 395 rows. The engine uses the Kenward-Roger method, which adjusts the standard errors and degrees of freedom of fixed effects in mixed models for small samples.

How strong is the clustering?

The engine reports the intraclass correlation (ICC): the share of the variation left after ozone and time that lies between trees rather than between measurements on the same tree. Here it is 0.904 (95% profile interval 0.869 to 0.932). The tree-level standard deviation is 0.608, against 0.198 for the measurement-level residual. Trees differ from each other far more than one tree’s measurements differ from the fitted line, which is exactly the situation where counting rows overstates the information.

Why the two intervals move in opposite directions

The ozone interval widened, but the Time standard error fell from 0.00085 to 0.00027. Both changes are correct, and they come from the same mechanism.

A between-cluster predictor (ozone) is compared across trees. Every tree’s own level adds noise to that comparison, and the ordinary regression pretends the noise averages out over 395 rows when it only averages out over 79 trees. Moulton’s paper on regressions with group-level predictors makes the same point for economic data: ordinary least squares can give standard errors that are seriously biased downward.

A within-cluster predictor (Time) is compared inside each tree. Once each tree has its own intercept, tree-to-tree differences drop out of that comparison, much as a paired t-test removes the differences between subjects. The ordinary regression leaves them in the residual, which inflates the Time standard error.

So the rule of thumb “ignoring clustering makes standard errors too small” holds for predictors that differ between clusters. For predictors that vary within clusters, the error can go the other way.

A caveat about the model itself

A random intercept lets each tree start at its own level but forces every tree to grow at the same rate. If trees grow at different rates, the Time slope varies between trees, and the random-intercept model’s Time standard error is likely too small. The next check is a model with a random slope for Time, (1 + Time | tree), before relying on that interval.

There is also a possible higher level. The data say the trees were grown in ozone-enriched chambers, but the file does not record which chamber each tree was in. If several trees shared a chamber, the chamber is a further cluster, and the ozone comparison has fewer independent units than 79 trees.

The other checks

The engine’s assumption checks on the mixed model: Shapiro-Wilk on the residuals was failed (W 0.970, p < 0.001), and Shapiro-Wilk on the tree intercepts did not detect non-normality (p 0.173), which is not evidence that they are normal. With 395 residuals the test detects even small departures, so the useful next step is a residual plot by tree rather than a verdict from the p-value.

How to check independence

Independence is decided by the design, before any output exists.

1. Name the unit of analysis. What was randomised, or sampled? In Sitka the treatment was assigned to trees (or, if trees shared chambers, to chambers), not to measurement days.
2. Count rows per unit. If any unit contributes more than one row, the rows are not independent. A column of IDs that repeat is the clearest sign.
3. Look for shared conditions. Same batch, same day, same observer, same site. Each is a candidate cluster.
4. Check the predictors against the clusters. Note which predictors vary between clusters and which within, because their standard errors are affected differently.

A test can support this judgement but not replace it. Durbin and Watson built their test for serial correlation: correlation between successive residuals in one ordering, which matters for data collected over time. By construction it cannot see clustering in a different grouping. A residual plot coloured by cluster, or the ICC from a mixed model, shows how much clustering there is once you know where to look.

Random sampling or assignment of independent units is what makes the assumption plausible. Convenience samples, snowball recruitment and families enrolled together invite dependence.

What to do about it

• Model the clustering. A mixed model with a random effect for each cluster, as above, is the general tool. The repeated-measures ANOVA vs mixed model guide explains when each fits.
• Analyse one number per unit. Average each tree’s measurements, or fit each tree’s growth slope, then compare those numbers across trees with an ordinary test. It is less efficient but hard to get wrong.
• Use cluster-robust standard errors. They correct the standard errors of an ordinary regression without modelling the structure. Cameron and Miller’s practitioner’s guide recommends them when the number of clusters is large and treats few clusters as a separate complication.
• Use population-averaged models. Generalised estimating equations (GEE) estimate population-averaged effects for continuous, binary or count outcomes (Zeger, Liang and Albert 1988), with standard errors that stay consistent under mild assumptions about the within-cluster correlation (Liang and Zeger 1986). Like cluster-robust errors, they need many clusters.
• For correlation in time or space, where neighbouring measurements are alike, use a model with a correlated error structure instead of ordinary regression.
• For two measurements per unit, a paired test is the simplest version of the same idea. The paired t-test vs Wilcoxon guide covers it.

Reporting independence in a Methods section

Say what the unit of analysis was, how many units there were against how many rows, and how clustering was handled. For example: “Each of 79 trees was measured five times; tree was modelled as a random intercept. Chamber was not recorded, so chamber-level clustering could not be modelled.” State an unresolved cluster like that as a limitation rather than leaving it out.

For the assumptions that sit alongside this one, normality and equal variances, see the guide to parametric test assumptions. The t-tests explained post shows how independence enters the simplest tests, and the guide to choosing a statistical test starts from the design questions above. The method reference pages describe the mixed model and linear regression outputs in detail.

Rerunning the example

Below the article is the R script exactly as the engine generated it, with the seed and the bundled data path, and the output it printed. It refits the mixed model with lme4 (Bates and colleagues), computes the Kenward-Roger F tests and the coefficient table with intervals, and computes the ICC with its profile interval. It does not recompute the assumption checks or the variance-component intervals.

The script does not include the ordinary regression. Those figures come from a separate engine run on the same data; you can reproduce them with lm(size ~ treat + Time, data = data) and confint() on the same file. The guide to reproducible statistical analysis explains what a reviewer needs in order to rerun a result like this one.

# Results

- Test: Linear mixed-effects model (random intercept)
- n: 395
- statistic: 2.019
- p-value: 0.159
- effect size: intraclass correlation = 0.9041 (95% CI 0.8692 to 0.9319)
- 95% CI: [-0.5071, 0.08477]
- seed: 20261022

## Assumption checks
- normality (Shapiro-Wilk): passed=FALSE (evidence of non-normality)
- random-intercept normality (Shapiro-Wilk on BLUPs): passed=TRUE (no evidence against normality)

## Term tests
| Term | Test | Statistic | df | Residual df | p-value |
|---|---|---|---|---|---|
| treat | F | 2.019 | 1 | 77 | 0.159 |
| Time | F | 2284 | 1 | 315 | < 0.001 |

## Random effects
- degrees of freedom: Kenward-Roger
- cluster sizes: 5 to 5 rows (median 5); 0 cluster(s) with a single row
- predictors: treat (between cluster), Time (within cluster)
| Group | Term | With | Variance | SD or correlation | 95% CI |
|---|---|---|---|---|---|
| tree | (Intercept) | - | 0.3696 | 0.6079 | 0.5154 to 0.7095 |
| Residual | - | - | 0.03921 | 0.198 | 0.1832 to 0.2142 |
- variance component intervals: profile likelihood

## Power
- Power is not calculated for this test.

# Celeus - reproduce.R  (run from this package folder; bundled data in data/)
# Test: Linear mixed-effects model (random intercept)  |  engine 0.30.0  |  catalog c39-d3cd
# Research context: not provided for this run.
# Reproduces the recorded result deterministically.

set.seed(20261022)
.data_path <- "data/independence-assumption.csv"
if (!file.exists(.data_path)) stop("Run reproduce.R from the package folder (the bundled data is in data/).")
data <- readr::read_csv(.data_path, show_col_types = FALSE)  # needs the readr package
clean_cells <- function (df) 
{
    ctrl_collapse <- "[\r\n\t\v\f]+"
    zero_width_chars <- intToUtf8(c(8203L, 65279L), multiple = TRUE)
    uni_space_chars <- intToUtf8(c(160L, 8239L, 8199L), multiple = TRUE)
    clean_text <- function(x) {
        out <- x
        for (ch in zero_width_chars) out <- gsub(ch, "", out, 
            fixed = TRUE)
        out <- gsub(ctrl_collapse, " ", out)
        for (ch in uni_space_chars) out <- gsub(ch, " ", out, 
            fixed = TRUE)
        trimws(gsub(" +", " ", out))
    }
    changed_cols <- character(0)
    cells_cleaned <- 0L
    retyped <- character(0)
    for (j in seq_along(df)) {
        col <- df[[j]]
        if (!is.character(col)) 
            next
        cleaned <- clean_text(col)
        n_changed <- sum(!is.na(col) & !is.na(cleaned) & col != 
            cleaned, na.rm = TRUE)
        if (n_changed > 0L) {
            changed_cols <- c(changed_cols, names(df)[j])
            cells_cleaned <- cells_cleaned + n_changed
        }
        df[[j]] <- cleaned
        non_blank <- cleaned[!is.na(cleaned) & nzchar(cleaned)]
        if (length(non_blank) && all(grepl("^[+-]?(0|[1-9][0-9]*)([.][0-9]+)?$", 
            non_blank))) {
            parsed <- suppressWarnings(as.numeric(cleaned))
            was_blank <- is.na(cleaned) | !nzchar(cleaned)
            lost <- (is.na(parsed) | is.infinite(parsed)) & !was_blank
            if (!any(lost)) {
                finite <- parsed[!is.na(parsed)]
                all_whole <- all(finite == floor(finite))
                in_range <- all(abs(finite) <= .Machine$integer.max)
                df[[j]] <- if (all_whole && in_range) 
                  as.integer(parsed)
                else parsed
                retyped <- c(retyped, names(df)[j])
            }
        }
    }
    names(df) <- clean_text(names(df))
    if (length(changed_cols)) {
        attr(df, "ingest_clean_columns") <- changed_cols
        attr(df, "ingest_clean_cells_cleaned") <- cells_cleaned
    }
    if (length(retyped)) 
        attr(df, "ingest_clean_retyped") <- retyped
    df
}
data <- clean_cells(data)
# Complete-case analysis: drop rows missing any analysis value.
.acols <- c("size", "treat", "Time", "tree")
data <- data[stats::complete.cases(data[, .acols, drop = FALSE]), , drop = FALSE]
result <- { .full <- lme4::lmer(`size` ~ `treat` + `Time` + (1 | `tree`), data = data, REML = TRUE); .tt <- { .X <- lme4::getME(.full, "X"); .asg <- attr(.X, "assign"); .I <- diag(ncol(.X)); t(vapply(sort(unique(.asg[.asg > 0L])), function(j) unlist(pbkrtest::KRmodcomp(.full, .I[.asg == j, , drop = FALSE])$test["Ftest", c("stat", "ndf", "ddf", "p.value")]), numeric(4L))) }; .p <- .tt[1L, "p.value"]; .ct <- { .V <- pbkrtest::vcovAdj(.full); .V0 <- as.matrix(stats::vcov(.full)); .L <- diag(ncol(.V0)); .b <- lme4::fixef(.full); .se <- sqrt(diag(as.matrix(.V))); .df <- vapply(seq_len(ncol(.V0)), function(i) pbkrtest::Lb_ddf(L = .L[i, , drop = FALSE], V0 = .V0, Vadj = .V), numeric(1L)); .tv <- .b / .se; .q <- stats::qt(1 - (1 - 0.95) / 2, .df); cbind(Estimate = .b, `Std. Error` = .se, df = .df, `t value` = .tv, `Pr(>|t|)` = 2 * stats::pt(abs(.tv), .df, lower.tail = FALSE), conf.low = .b - .q * .se, conf.high = .b + .q * .se) }; list(model = summary(.full), p_value = unname(.p), estimate = unname(.ct[2L, "Estimate"]), conf_int = unname(.ct[2L, c("conf.low", "conf.high")]), icc = { .vc <- as.data.frame(lme4::VarCorr(.full)); .vc$vcov[1L] / (.vc$vcov[1L] + .vc$vcov[2L]) }, icc_ci = tryCatch({ .dev <- lme4::lmer(`size` ~ `treat` + `Time` + (1 | `tree`), data = data, REML = TRUE, devFunOnly = TRUE); .th <- lme4::getME(.full, "theta")[[1L]]; .g <- function(t) .dev(t) - .dev(.th) - stats::qchisq(0.95, 1); .lo <- if (.g(0) <= 0) 0 else stats::uniroot(.g, c(0, .th), tol = 1e-10)$root; .up <- max(.th, 1); .k <- 0L; while (.g(.up) <= 0 && .k < 60L) { .up <- 2 * .up; .k <- .k + 1L }; .hi <- if (.g(.up) > 0) stats::uniroot(.g, c(.th, .up), tol = 1e-10)$root else NA_real_; c(.lo, .hi)^2 / (1 + c(.lo, .hi)^2) }, error = function(e) c(NA_real_, NA_real_)), coefficients = .ct, term_tests = .tt) }
print(result)