Answer

How do you interpret the ICC when random effects are nested?

Inspired by a question on Cross Validated · October 16, 2026

multilevel modelsintraclass correlation

The short answer

Suppose you take repeated measurements at several stations inside each of many sites (lakes, schools, hospitals) and fit an intercept-only mixed model with a random effect for site and for station within site. The model splits the total variance into three parts: between sites, between stations within a site, and between readings at the same station. The intraclass correlation (ICC) is a ratio of these parts, but which parts go in the numerator depends on the question you are asking.

• "Are readings from the same site more alike, whatever the station?" Use site variance / total variance. This is the correlation between two readings taken at the same site but at different stations.
• "Are readings from the same station alike?" Use (site variance + station variance) / total variance. Two readings from one station share both random effects, so they are more strongly correlated.
• "Within a site, how much do stations differ?" Use station variance / (station variance + residual variance), which ignores between-site differences altogether.

Many summaries (for example the "adjusted ICC" from the R package performance) add up every random-effect variance, so with nested effects they report the same-station correlation. That number is always the larger one, and it overstates how alike two readings from the same site are when they come from different stations.

Why the components add up this way

The model writes each reading as an overall mean plus a site effect, plus a station effect, plus a residual, with all three independent and with variances σ²_site, σ²_station and σ²_resid. The total variance of a single reading is their sum. The covariance between two readings equals the variance of whatever random effects they share:

• Same site, different stations: they share only the site effect, so the covariance is σ²_site.
• Same station: they share the site and the station effect, so the covariance is σ²_site + σ²_station.
• Different sites: they share nothing, so the covariance is zero.

Dividing each covariance by the total variance turns it into a correlation, which is all an ICC is. So with nested random effects you get one ICC per level of sharing, and it is worth reporting the variance components themselves so readers can form whichever ratio they need. Comparing the size of the site and station variances (as people often do by eye) tells you where the variation sits, but the ICC is what translates that into how alike two readings are.

See it in R and Python

The code simulates a balanced design with 30 sites, 6 stations per site and 20 readings per station. The true variances are 0.49 (site), 0.36 (station) and 0.81 (residual), so the true ICCs are about .30 for the same site, .51 for the same station, and .31 for stations within a site. It estimates the components from the nested ANOVA mean squares, which can be computed by hand in any language, and then checks them against a REML fit from the lme4 package. For balanced data like these the two methods give the same estimates as long as none comes out negative.

set.seed(677417)

# Balanced three-level design: 30 sites, 6 stations per site, 20 readings per station
a <- 30; b <- 6; n <- 20
site    <- rep(1:a, each = b * n)
station <- rep(1:(a * b), each = n)
# True SDs: site 0.7, station 0.6, residual 0.9 (variances 0.49, 0.36, 0.81)
y <- 1 + rnorm(a, 0, 0.7)[site] + rnorm(a * b, 0, 0.6)[station] + rnorm(a * b * n, 0, 0.9)

# Variance components from the nested ANOVA mean squares
m_st <- tapply(y, station, mean); m_si <- tapply(y, site, mean); gm <- mean(y)
ms_site <- b * n * sum((m_si - gm)^2) / (a - 1)
ms_stat <- n * sum((m_st - rep(m_si, each = b))^2) / (a * (b - 1))
ms_res  <- sum((y - m_st[station])^2) / (a * b * (n - 1))
v <- c(site = (ms_site - ms_stat) / (b * n),
       station = (ms_stat - ms_res) / n,
       residual = ms_res)
round(v, 3)

# Three different "ICCs" from the same model
tot <- sum(v)
round(c(same_site_diff_station = v[["site"]] / tot,
        same_station           = (v[["site"]] + v[["station"]]) / tot,
        station_within_site    = v[["station"]] / (v[["station"]] + v[["residual"]])), 3)

# The same model fitted by REML with lme4
library(lme4)
fit <- lmer(y ~ 1 + (1 | site / station),
            data = data.frame(y, site = factor(site), station = factor(station)))
vc <- as.data.frame(VarCorr(fit))
setNames(round(vc$vcov, 3), vc$grp)

import numpy as np

rng = np.random.default_rng(677417)

# Balanced three-level design: 30 sites, 6 stations per site, 20 readings per station
a, b, n = 30, 6, 20
site = np.repeat(np.arange(a), b * n)
station = np.repeat(np.arange(a * b), n)
# True SDs: site 0.7, station 0.6, residual 0.9 (variances 0.49, 0.36, 0.81)
y = (1 + rng.normal(0, 0.7, a)[site] + rng.normal(0, 0.6, a * b)[station]
     + rng.normal(0, 0.9, a * b * n))

# Variance components from the nested ANOVA mean squares
m_st = y.reshape(a * b, n).mean(axis=1)
m_si = y.reshape(a, b * n).mean(axis=1)
gm = y.mean()
ms_site = b * n * np.sum((m_si - gm) ** 2) / (a - 1)
ms_stat = n * np.sum((m_st - np.repeat(m_si, b)) ** 2) / (a * (b - 1))
ms_res = np.sum((y - m_st[station]) ** 2) / (a * b * (n - 1))
v_site = (ms_site - ms_stat) / (b * n)
v_station = (ms_stat - ms_res) / n
v_res = ms_res
print("site", round(v_site, 3), "station", round(v_station, 3), "residual", round(v_res, 3))

# Three different "ICCs" from the same model
tot = v_site + v_station + v_res
print("same_site_diff_station", round(v_site / tot, 3))
print("same_station", round((v_site + v_station) / tot, 3))
print("station_within_site", round(v_station / (v_station + v_res), 3))

    site  station residual 
   0.448    0.371    0.818 
same_site_diff_station           same_station    station_within_site 
                 0.274                  0.500                  0.312 
station:site         site     Residual 
       0.371        0.448        0.818 

The output is from one seeded run of the R code. statsmodels is not used, so the Python version computes only the ANOVA estimates, not a REML fit. Its random numbers differ from R's, so its figures differ, and the site variance differs most (0.80 instead of 0.45), because it is estimated from only 30 sites. The pattern is the same: the same-station ICC (.59 in Python) is well above the same-site ICC (.41), and the within-site station ICC is close to .30 in both.

• The two methods agree. The hand-computed ANOVA components (0.448, 0.371, 0.818) match the lme4 REML estimates to three decimals.
• One model, three ICCs. Two readings from the same site but different stations correlate at .274; two readings from the same station correlate at .500; and within a site, stations account for .312 of the variance.
• A software summary that adds up all random effects would print .500 here, almost twice the answer to the question "are readings from the same site alike regardless of station?".

Is the clustering big enough to matter?

• There is no universal cut-off. An ICC of .27 means a site explains about a quarter of the variation in single readings, which is a lot in most fields. Even an ICC of .05 matters for inference when each cluster holds many readings, because the effective sample size shrinks well below the number of readings.
• For a formal test that the site variance is above zero, compare the model with and without the site random effect using a likelihood ratio test. Because a variance cannot be negative, the usual chi-square p-value is conservative; halving it is the common correction for one variance component.
• The ICC belongs to the scale you analysed. If you transformed the outcome, for example with log(x + 1) because of zeros and extreme values, the ICC describes the transformed values, and it can change with the constant you add. See Log-transforming data: when it helps, and what it changes.
• Unbalanced designs are fine for a REML fit; only the hand-computed ANOVA formulas above assume equal numbers per group.

For a refresher on what random intercepts are, see Fixed vs random effects in mixed models. Not sure whether you need a mixed model at all? Try the test chooser, or browse the DASS blog for plain-language guides.

How to report the ICC in APA style (7th edition)

Name the levels, give the variance components, and say which ICC you mean. Using the R example:

"An intercept-only model with random intercepts for sites and for stations nested within sites, fitted by REML, showed that 27% of the variance lay between sites (ICC = .27) and a further 23% between stations within sites. Two readings from the same station therefore correlated at ICC = .50."

Put the variance components (here 0.448 for sites, 0.371 for stations and 0.818 for the residual) in a table or a note, with the number of readings, stations and sites. ICCs cannot exceed 1, so they have no leading zero.

Related tools and guides

• Multilevel power and MDES calculator
• Fixed vs random effects in mixed models
• Log-transforming data: when it helps, and what it changes
• Which statistical test should I use?
• Intraclass correlation (Wikipedia)

More answered questions

• Fixed vs random effects in mixed models: what is the difference?

Written by AskStats with AI assistance. This is general information, not advice for your specific data or study. When the results matter, check your approach with a qualified statistician.