Estimating Intraclass Correlation Reliably

Why the Estimation Method Matters

The intraclass correlation (ICC) describes how similar observations from the same donor are, relative to the total variation among observations. In single-cell data, a large ICC means that two cells from one donor tend to be more alike than two cells drawn from different donors, under the model being used. The estimate is not just a property of a calculation: it depends on the model, the number of donors, and how many cells were observed from each donor.

This tutorial focuses on estimating the ICC and judging whether the estimate is stable enough to interpret. We use a one-way random-intercept model, the same basic framework used in the earlier material on separating within- and between-donor variance. The formulas for method-of-moments estimation below assume a balanced design: the same number of cells is observed for every donor.

Model and assumptions: For donor $i=1,\ldots,D$ and cell $j=1,\ldots,m$, suppose $$ Y_{ij}=\mu+b_i+\epsilon_{ij}, $$ where $\mu$ is the overall mean, donor effects $b_i$ are independent with mean zero and variance $\tau^2$, and cell-level errors $\epsilon_{ij}$ are independent with mean zero and variance $\sigma^2$. Donor effects and errors are mutually independent. For the Gaussian random-intercept model, additionally assume $b_i\sim N(0,\tau^2)$ and $\epsilon_{ij}\sim N(0,\sigma^2)$. The variance components satisfy $\tau^2\geq 0$ and $\sigma^2\geq 0$.

The model treats donors as independent sampling units and cells as observations nested within donors. It assumes a common residual variance and a common donor variance for the observations being modeled. These assumptions are modeling choices, not facts guaranteed by collecting single-cell data. For example, if cell type or condition changes the variance structure, a single common ICC may not summarize the dependence adequately.

Definition: When $\tau^2+\sigma^2>0$, the ICC is $$ \rho=\frac{\tau^2}{\tau^2+\sigma^2}. $$ It is the correlation between two distinct cells from the same donor under the random-intercept model. With nonnegative variance components and $\tau^2+\sigma^2>0$, $0\leq\rho\leq1$. If both variance components are zero, the ratio is undefined; there is no variation from which to define this ICC.

Method of Moments in a Balanced Design

Method-of-moments estimation matches observed mean squares to their expected values under the model. In a balanced design, define the donor mean $\bar{Y}_{i\cdot}=m^{-1}\sum_{j=1}^{m}Y_{ij}$ and the grand mean $\bar{Y}_{\cdot\cdot}=(Dm)^{-1}\sum_{i=1}^{D}\sum_{j=1}^{m}Y_{ij}$. The between-donor and within-donor mean squares are

MS_B=\frac{m\sum_{i=1}^{D}(\bar{Y}_{i\cdot}-\bar{Y}_{\cdot\cdot})^2}{D-1}, \qquad MS_W=\frac{\sum_{i=1}^{D}\sum_{j=1}^{m}(Y_{ij}-\bar{Y}_{i\cdot})^2}{D(m-1)}.

The numerator for $MS_B$ measures variation among donor means, scaled by the number of cells per donor. The numerator for $MS_W$ measures variation among cells around their own donor means. Under the model,

E(MS_W)=\sigma^2,\qquad E(MS_B)=\sigma^2+m\tau^2.

The second expectation has two contributions: within-donor variation and donor-to-donor variation, with the latter multiplied by $m$. Matching these expectations gives $\hat{\sigma}^2=MS_W$ and the unconstrained estimate $\hat{\tau}^2=(MS_B-MS_W)/m$. Substituting those estimates into the ICC formula gives the method-of-moments ICC.

\hat{\rho}_{MM} =\frac{MS_B-MS_W}{MS_B+(m-1)MS_W}, \qquad\text{when }MS_B\geq MS_W.

If $MS_B<MS_W$, the unconstrained method-of-moments estimate of $\tau^2$ is negative. A negative variance is not allowed by the model. A common constrained method-of-moments report sets $\hat{\tau}^2=0$ and reports $\hat{\rho}_{MM}=0$. This boundary estimate should be understood as evidence that the data do not support a positive between-donor variance with this estimator—not as proof that the population ICC is exactly zero.

Worked Example: A Positive Method-of-Moments Estimate

Suppose a hypothetical balanced data set has $D=4$ donors and $m=3$ cells per donor. For simplicity, the measured values are: donor A, 8, 9, 10; donor B, 10, 11, 12; donor C, 12, 13, 14; and donor D, 14, 15, 16. These are hypothetical measurements on a common scale.

The donor means are 9, 11, 13, and 15, and the grand mean is 12. The within-donor squared deviations sum to $2$ for each donor, so the within-donor sum of squares is $8$. Thus $MS_W=8/[4(3-1)]=8/8=1$.

For the between-donor mean square, the squared deviations of the donor means from 12 sum to $(9-12)^2+(11-12)^2+(13-12)^2+(15-12)^2=20$. Multiplying by $m=3$ gives a between-donor sum of squares of 60, so $MS_B=60/(4-1)=20$.

The variance-component estimates are $\hat{\sigma}^2=1$ and $\hat{\tau}^2=(20-1)/3=6.333$. Therefore $\hat{\rho}_{MM}=6.333/(6.333+1)=0.8636$. Equivalently, the direct formula gives $(20-1)/[20+(3-1)(1)]=19/22=0.8636$. The estimate is high because the observed between-donor variation is large relative to the within-donor variation. Four donors are still a small basis for judging how precisely this population quantity has been estimated.

Model-Based Estimation and Boundary Behavior

A model-based approach estimates variance components by optimizing a likelihood under the random-intercept model. Maximum likelihood (ML) estimates parameters by maximizing the likelihood of the observed data. Restricted maximum likelihood (REML) instead constructs a likelihood from contrasts that remove the fixed overall mean; this accounts for the estimation of that mean when estimating variances. Both approaches enforce the model’s nonnegative variance-component constraints.

In a balanced Gaussian design, the unconstrained REML estimates equal the method-of-moments estimates above. With the nonnegative variance-component constraint, they coincide when the unconstrained estimate of $\tau^2$ is nonnegative. If that estimate is negative, the constrained REML fit places $\tau^2$ at the boundary zero, rather than reporting a negative variance; the residual variance is estimated as part of that constrained fit. Thus it is important to distinguish the unconstrained moment calculation from the constrained model fit.

For unbalanced cell counts, the balanced-design mean-square formula should not be applied as if every donor had the same $m$. A model-based fit can use the observed donor-specific cell counts under its stated assumptions. The estimates can nevertheless be sensitive to imbalance, especially if some donors contribute very few cells or if cell count is related to donor characteristics or the outcome.

Worked Example: A Negative Unconstrained Estimate

Consider a hypothetical balanced data set with four donors and three cells per donor. The donor means are 10, 10.1, 9.9, and 10.0. Within each donor, the three values are the donor mean minus 1, the donor mean, and the donor mean plus 1. These data have $MS_W=1$, since each donor contributes within-donor squared deviations totaling $2$, for a total of $8$, and $8/[4(3-1)]=1$.

The grand mean is 10. The squared deviations of the donor means sum to $0.01+0.01=0.02$, so the between-donor sum of squares is $3(0.02)=0.06$, and $MS_B=0.06/(4-1)=0.02$.

The unconstrained estimate is $\hat{\tau}^2=(0.02-1)/3=-0.3267$. This is not an admissible variance-component estimate. The constrained method-of-moments estimate is $\hat{\tau}^2=0$, giving $\hat{\rho}_{MM}=0$. A constrained REML fit likewise cannot assign a negative variance to donors; it fits on the zero-variance boundary if that is the likelihood optimum. The result indicates that these observations do not resolve a positive donor component with much support. It does not establish that the true ICC is zero.

Worked Example: Interpreting a Model-Based Fit

Suppose a hypothetical Gaussian random-intercept REML fit to expression measurements reports $\hat{\tau}^2=4.00$ and $\hat{\sigma}^2=2.25$, both in squared expression units. The implied ICC is

\hat{\rho}_{REML} =\frac{4.00}{4.00+2.25} =\frac{4.00}{6.25} =0.6400.

Under the fitted model, the estimated correlation between two distinct cells from the same donor is 0.6400. This is a dimensionless correlation, not a statement that 64% of each cell’s expression is caused by the donor. It also does not describe causal influence. The variance components depend on the modeled outcome scale, so changing that scale may change their numerical values even when the ICC remains comparable under a common linear rescaling.

Why Few Donors and Few Cells Make Estimates Unstable

The donor count and the cell count per donor provide different kinds of information. With $D$ donors, the between-donor mean square has $D-1$ degrees of freedom. When $D$ is small, that estimate of between-donor variation can change substantially if the observed donors happen to be more or less dispersed than the target population. Adding many cells to those same donors can improve each donor’s mean, but it does not create more independent donors.

The model also gives $\operatorname{Var}(\bar{Y}_{i\cdot})=\tau^2+\sigma^2/m$. Increasing $m$ reduces the contribution of within-donor noise to the donor mean, but does not remove the between-donor term $\tau^2$. With few cells per donor, the within-donor variance is itself estimated imprecisely. With one cell per donor, there is no within-donor replication, so this design alone cannot separate $\tau^2$ from $\sigma^2$.

An ICC point estimate should therefore be examined alongside its uncertainty. A model-based interval, for example one obtained by profiling the likelihood for $\rho$, can show whether a high point estimate is compatible with a broad range of values. Near the zero-variance boundary, the sampling distribution can be asymmetric, so a routine symmetric interval may be misleading. A donor-level resampling analysis is another possible sensitivity check when its assumptions and small donor count are considered carefully. Neither an interval nor a resampling procedure manufactures information absent from the design.

Useful diagnostics include checking the number of donors, the distribution of cell counts by donor, and whether the variance estimates are at or near a boundary. Compare method-of-moments and model-based estimates when both are appropriate: a large discrepancy can flag sensitivity to assumptions, imbalance, or the fitting method. Also inspect whether a few donors dominate the estimated between-donor spread. An estimate that changes sharply under reasonable modeling choices should be reported as unstable rather than presented with false precision.

State the target and model.
Specify the outcome scale, donor and cell indexing, and the assumptions that define the variance components and ICC.

Check the design.
Count donors and cells per donor. Use balanced-design moment formulas only when the cell count is common across donors.

Estimate and check constraints.
Calculate mean squares for the balanced method-of-moments estimate, or fit the stated model. Do not report a negative variance component as valid.

Assess stability before interpretation.
Examine uncertainty, boundary behavior, donor count, cell-count imbalance, and sensitivity to the estimation method.

Comparison and Common Pitfalls

Approach	Strength	Important limitation
Method of moments	Transparent calculation from between- and within-donor mean squares in a balanced design.	Can yield a negative unconstrained variance estimate; the simple formula does not directly handle unequal cell counts.
REML	Model-based estimation with nonnegative variance constraints; can use unequal cell counts under the specified model.	Still depends on the model and available donors; a boundary fit or precise-looking output does not guarantee a precise ICC.
ML	Likelihood-based estimation of model parameters.	Variance estimates can differ from REML, particularly in small samples; fitting by ML does not remove limited donor information.

A common mistake is to treat the number of cells as the effective sample size for estimating donor heterogeneity. For example, “thousands of cells means the ICC is precise” ignores that those cells are clustered within donors. The correction is to report both the donor count and the cells-per-donor pattern, then assess uncertainty at the level relevant to the model. Another mistake is interpreting a zero boundary estimate as proof of no donor variation. It means the fitted constrained procedure selected the boundary for these data; few donors may leave substantial uncertainty about the population ICC.

The ICC is useful for describing dependence, but it is not a universal quality score and it is not enough to diagnose every feature of single-cell data. It does not by itself establish whether a downstream analysis has valid uncertainty estimates, and it cannot compensate for a poorly specified outcome model. Interpret it in the context of the outcome, the sampling design, the model assumptions, and the uncertainty in its estimate.

Check Your Understanding

Try each question, then compare your reasoning with the solutions provided.

In the random-intercept model, what does the ICC represent, and under what condition is it defined?
Solution: It is the correlation between two distinct cells from the same donor under the model. It is defined when $\tau^2+\sigma^2>0$; if both components are zero, the ratio is undefined.
A balanced design has $MS_B=8$, $MS_W=2$, and $m=3$. Calculate the method-of-moments ICC.
Solution: $\hat{\rho}_{MM}=(8-2)/[8+(3-1)2]=6/12=0.5000$.
Why is a negative unconstrained estimate of $\tau^2$ not a valid population variance estimate?
Solution: Variances in the model are constrained to be nonnegative. A negative moment estimate signals that the observed between-donor mean square is below the within-donor mean square; a constrained procedure places the estimate at the zero boundary rather than interpreting the negative value as a variance.
Why can observing more cells from the same donors fail to resolve uncertainty caused by having few donors?
Solution: More cells can reduce within-donor sampling noise, but the independent donor count remains unchanged. The between-donor component is still informed by the observed set of donors.
In a balanced Gaussian design, when do unconstrained REML and method-of-moments estimates coincide, and what changes if the unconstrained donor variance estimate is negative?
Solution: The unconstrained estimates coincide. With nonnegative variance constraints, they coincide when the unconstrained donor variance estimate is nonnegative. If it is negative, constrained REML fits the donor variance at the zero boundary and estimates the remaining parameters under that constraint.