Separating Algorithmic From Sampling Variability

Why Separate the Sources of Instability?

A clustering result can change for at least two distinct reasons. On the same subjects and with the same input data, a randomized optimizer may reach different solutions from different initializations. Even if optimization were completely reproducible, a different sample of subjects could produce a different partition. These sources call for different responses: more starts may address optimization instability, whereas more representative data or a more robust analysis may be needed to address sampling instability.

A comparison across one run per resample mixes the two sources: a change in the result might reflect a new sample, a new initialization, or both. The design here runs several randomized starts within each subject sample, so the within-sample variation informs algorithmic instability and the variation among sample means informs sampling instability. The analysis is performed on a defined numerical summary of each fitted result; it does not decompose an entire partition into two variance components.

Definition: Let $b=1,\ldots,B$ index independently generated subject samples, and let $r=1,\ldots,R$ index randomized initializations within sample $b$. Let $Y_{br}$ be a scalar summary of the clustering result from that sample and initialization. Under the balanced nested random-effects model, $$ Y_{br}=\mu+S_b+A_{br}, $$ where $\mu$ is the overall mean, $S_b$ is the sampling effect, and $A_{br}$ is the within-sample algorithmic effect. The model assumes independent effects with mean zero, variances $\sigma_S^2$ and $\sigma_A^2$, and a common variance structure across samples. The two effects are independent of one another.

Define a Comparable Output Before Repeating the Fits

The $Y_{br}$ values must measure the same thing on the same scale. One option is a label-invariant disagreement score between each fitted partition and a prespecified reference partition, evaluated on a common set of subjects. For example, define the score as the fraction of subject pairs whose co-clustering status differs from the reference. It lies between 0 and 1, with larger values indicating more disagreement. The reference and evaluation set should be fixed before comparing runs.

For a clustering method that does not naturally assign new subjects to a fitted partition, specify a consistent evaluation rule. This could use a prespecified out-of-sample assignment rule, or restrict comparisons to a common set of subjects that can be evaluated in every fit. If the evaluation set or the definition of the score changes across resamples, variation in $Y_{br}$ may reflect that change rather than either source of instability. A single score is only one view of a partition: it is often useful to repeat the decomposition for several well-defined summaries or pairwise co-clustering indicators.

Use the same clustering pipeline within each sample. If preprocessing is part of the intended analysis, refit it within each resample as specified by the analysis design. To attribute within-sample differences specifically to random initialization, hold other sources of randomness fixed or control them separately. Otherwise, the within-sample component represents all uncontrolled within-sample randomness, not initialization alone.

Worked Example: A Pairwise Co-Clustering Output

Suppose four subject samples are generated, and three randomized starts are run within each. For a fixed pair of evaluation subjects, record $Y_{br}=1$ if the pair is placed in the same cluster and $Y_{br}=0$ otherwise. The hypothetical results are: sample 1, $1,1,0$; sample 2, $1,1,1$; sample 3, $0,0,1$; and sample 4, $0,0,0$. This pairwise indicator is a simple, comparable output, though a full analysis would usually examine many pairs or an aggregate score.

The sample means are $2/3,1,1/3,0$, and the grand mean is $\bar{Y}_{..}=6/12=0.5$. Within sample 1, the sum of squared deviations is $2(1-2/3)^2+(0-2/3)^2=2/3$; sample 3 has the same value, and samples 2 and 4 have zero. Thus $SS_A=4/3$, with $B(R-1)=4(2)=8$ degrees of freedom, so $MS_A=(4/3)/8=1/6$.

For the sample means, $SS_S=R\sum_b(\bar{Y}_{b.}-\bar{Y}_{..})^2=3(5/9)=5/3$, with $B-1=3$ degrees of freedom. Therefore $MS_S=(5/3)/3=5/9$. The estimated components are $\hat{\sigma}_A^2=1/6$ and $\hat{\sigma}_S^2=(5/9-1/6)/3=7/54\approx0.1296$. These are estimates in squared indicator units. They describe this pair under this sampling and fitting design, not the stability of every feature of the partitions.

Estimating the Variance Components

For the balanced design, calculate each sample mean $\bar{Y}_{b.}$, the grand mean $\bar{Y}_{..}$, and the two mean squares below. $SS_A$ measures dispersion among starts within samples; $SS_S$ measures dispersion among sample means, scaled by the number of starts.

SS_A=\sum_{b=1}^{B}\sum_{r=1}^{R}(Y_{br}-\bar{Y}_{b.})^2,\qquad MS_A=\frac{SS_A}{B(R-1)}

SS_S=R\sum_{b=1}^{B}(\bar{Y}_{b.}-\bar{Y}_{..})^2,\qquad MS_S=\frac{SS_S}{B-1}

The method-of-moments estimates are

\hat{\sigma}_A^2=MS_A,\qquad \hat{\sigma}_S^2=\frac{MS_S-MS_A}{R}.

Derivation: Within a sample, subtracting its mean removes $S_b$ and leaves $A_{br}-\bar{A}_{b.}$. Under independent, equal-variance algorithmic effects, the expected within-sample sum of squares is $(R-1)\sigma_A^2$. Summing over $B$ samples gives $E(MS_A)=\sigma_A^2$. For the sample means, $\bar{Y}_{b.}=\mu+S_b+\bar{A}_{b.}$, so their variance is $\sigma_S^2+\sigma_A^2/R$. The between-sample sum of squares includes a factor $R$, giving $E(MS_S)=\sigma_A^2+R\sigma_S^2$. Subtracting the first expectation from the second and dividing by $R$ yields the sampling-component estimator.

The estimated total variance of a result under this model is $\hat{\sigma}_A^2+\hat{\sigma}_S^2$. The corresponding square roots express typical variation in the units of $Y$, but they do not replace inspection of the actual scores or partitions. The estimates are conditional on the chosen sampling scheme, pipeline, evaluation rule, and output summary.

Worked Example: Decomposing Variation in a Disagreement Score

Consider hypothetical normalized disagreement scores, measured on a common evaluation set, from three subject samples with two starts per sample: sample 1 has $0.1,0.3$; sample 2 has $0.5,0.7$; and sample 3 has $0.2,0.4$. The sample means are $0.2,0.6,0.3$, and the grand mean is $0.3667$.

Within each sample, the two values differ from their mean by $0.1$, so each within-sample sum of squares is $0.1^2+0.1^2=0.02$. Hence $SS_A=0.06$, $df_A=3(2-1)=3$, and $MS_A=0.06/3=0.0200$ squared score units.

For the sample means, $SS_S=2[(0.2-0.3667)^2+(0.6-0.3667)^2+(0.3-0.3667)^2]\approx0.1733$. Thus $df_S=3-1=2$ and $MS_S\approx0.08667$. The component estimates are $\hat{\sigma}_A^2=0.0200$ and $\hat{\sigma}_S^2=(0.08667-0.0200)/2\approx0.03333$. The estimated total variance is $0.0200+0.03333=0.05333$ squared score units. The estimated standard deviations are $\sqrt{0.0200}\approx0.1414$ score units for within-sample algorithmic variation and $\sqrt{0.03333}\approx0.1826$ score units for sampling variation. In this design, the estimated sampling component is larger, but three samples provide limited information about that component.

Assumptions, Diagnostics, and Comparisons

The formulas above require a balanced design with the same number of starts per sample and a scalar response defined consistently across all fits. The random-effects interpretation also assumes independent subject-sample effects, independent algorithmic effects, and comparable variances across samples. These conditions are design assumptions, not facts guaranteed by the calculations. If samples are stratified by cohort or batch, generate them according to the intended stratified scheme; the sampling component then pertains to that scheme.

Inspect the individual $Y_{br}$ values and sample means, not only the two estimates. Look for outlying samples, a few extreme starts, different within-sample spreads, and patterns by cohort or batch. If within-sample spreads vary greatly, a single $\sigma_A^2$ may be misleading. If scores are highly discrete, such as a binary pairwise indicator, the normal random-effects model may be a rough approximation; the mean-square decomposition remains a descriptive method-of-moments summary but should not be treated as an exact probability model. A small number of subject samples also makes the sampling component imprecise.

A negative method-of-moments estimate of $\sigma_S^2$ can occur when $MS_S<MS_A$. Since a true variance cannot be negative, interpret it as no clear evidence for a positive sampling component under this design; sampling noise may explain the estimate. It is not negative variability. Report the raw estimate transparently. A constrained variance-components fit may place the estimate at zero, but that boundary estimate is a different estimation procedure and does not erase the uncertainty in the data.

The nested design is more informative than running one start per sample, which confounds the two sources. Running many starts on just one fixed data set can reveal optimization variation but says little about subject-sampling variation. Agreement measures such as the adjusted Rand index, Jaccard stability, or variation of information compare partitions directly; they are useful summaries but do not by themselves separate the sources. The variance-components approach instead quantifies how a chosen scalar output varies under a specified nested design.

Worked Example: Interpreting a Negative Sampling Estimate

Suppose three samples each receive two starts, giving scores $0.1,0.5$; $0.2,0.6$; and $0.3,0.7$. The sample means are $0.3,0.4,0.5$, with grand mean $0.4$. Each within-sample sum of squares is $0.2^2+0.2^2=0.08$, so $SS_A=0.24$, $df_A=3$, and $MS_A=0.08$.

The between-sample sum of squares is $SS_S=2[(0.3-0.4)^2+(0.4-0.4)^2+(0.5-0.4)^2]=0.04$, so $MS_S=0.04/2=0.02$. The raw sampling estimate is $(0.02-0.08)/2=-0.03$ squared score units. This is not a negative variance. It indicates no clear evidence for a positive sampling component from this small design: within-sample variation is larger than the observed variation among sample means. More samples and inspection of the raw results are appropriate before drawing a substantive conclusion.

Limitations and Practical Use

The decomposition is specific to the chosen score and sampling design. A stable average disagreement score can conceal an unstable rare subgroup, while one pairwise indicator cannot represent all aspects of a partition. Examine subgroup-level and feature-level summaries as appropriate, and avoid interpreting a variance component as biological validity or evidence that a discovered subgroup is meaningful.

The method-of-moments formulas are simplest for a balanced nested design. With unequal numbers of starts, missing runs, or additional sources of randomness, these formulas generally do not apply unchanged; a suitable mixed-effects model or a redesigned balanced experiment is needed. The sampling effects may also be dependent if samples are deliberately linked or reused in ways not represented by the model. In all cases, state which randomness is repeated, which is held fixed, and what population or resampling scheme the sampling component represents.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

Why are multiple starts needed within each subject sample?
In a balanced design with $R=4$, $MS_A=0.03$, and $MS_S=0.11$, calculate both variance-component estimates.
What does a raw estimate $\hat{\sigma}_S^2<0$ mean?
Why must the output summary be comparable across resamples and starts?

Solutions. (1) The repeated starts estimate within-sample variation attributable to optimization or other uncontrolled within-sample randomness; without them, that variation is mixed with sampling variation. (2) $\hat{\sigma}_A^2=0.03$, and $\hat{\sigma}_S^2=(0.11-0.03)/4=0.02$, in squared output units. (3) The mean square among sample means was smaller than the within-sample mean square; this is no clear evidence for a positive sampling component under the design, not negative variability. (4) If the score, evaluation subjects, or assignment rule changes, variation in the output can reflect those changes rather than differences caused by sampling or optimization.