Tutorials › Multiple Imputation by Chained Equations (MICE) › Few Clusters: Stabilizing Multilevel Imputation

Nested and Multilevel Data · Tutorial 349 of 497

Few Clusters: Stabilizing Multilevel Imputation

See why the number of clusters—not just the total sample size—limits variance estimation, and compare practical ways to stabilize multilevel imputation models.

Advanced 11 min read

What You'll Learn

  • Explain why between-cluster variance is hard to estimate when the number of clusters is small.
  • Derive a simple random-intercept variance estimate and identify its assumptions.
  • Interpret boundary estimates without treating them as proof that cluster variation is absent.
  • Compare parsimonious model choices with regularizing prior information.
  • Diagnose instability using model behavior, chain traces, and sensitivity checks.
  • Distinguish stabilizing an estimate from adding information to the data.

Why the Number of Clusters Matters

A data set can contain many individuals and still provide little information about how much outcomes vary between clusters. For a random-intercept model, information about the variance of the cluster effects comes mainly from the number of clusters, not simply from the total number of individuals. With few clusters, estimated variance components may change substantially when one cluster changes, and a fitted variance may land at the boundary of zero.

This matters in multiple imputation because uncertainty in the imputation-model parameters should be reflected in the generated values. If a random-intercept variance is poorly determined, the implied amount of between-cluster variation—and the degree of shrinkage for cluster-specific effects—can be uncertain too. A model that is too elaborate for the available clusters can therefore produce unstable imputations, even when the total number of records is large.

Definition: In this tutorial, “few clusters” means that the number of independent clusters is small enough, relative to the information needed to estimate the model, that variance components or other cluster-level parameters are weakly determined. There is no universal numerical cutoff: the problem depends on the model, cluster sizes, outcome variation, and the amount and pattern of missingness.

A Simple Model Shows the Source of Instability

Consider a balanced, continuous-outcome setting with \(J\) independent clusters and \(n\) individuals per cluster. Suppose \(Y_{ij}=\beta_0+u_j+e_{ij}\), where \(u_j\sim N(0,\tau^2)\), \(e_{ij}\sim N(0,\sigma^2)\), and all random effects and errors are mutually independent. Here \(\beta_0\) is the population intercept, \(\tau^2\) is the between-cluster variance, and \(\sigma^2\) is the within-cluster variance. Assume for the calculation below that the model is correctly specified and that \(\sigma^2\) is known.

The cluster mean contains both the cluster effect and the average individual error. Keeping the intercept in the expression is essential:

$$ \bar{Y}_j=\beta_0+u_j+\bar{e}_j,\qquad \operatorname{Var}(\bar{Y}_j)=\tau^2+\frac{\sigma^2}{n}. $$

The second equality follows because \(u_j\) and \(\bar e_j\) are independent, with variances \(\tau^2\) and \(\sigma^2/n\). If \(\bar Y_j\) is observed in every cluster, its sample variance \(S^2_{\bar Y}\) estimates the variance of the cluster means. A simple method-of-moments estimate subtracts the sampling variance of each cluster mean:

$$ \hat{\tau}^2_{\mathrm{mom}}= \max\left(0,\ S^2_{\bar Y}-\frac{\sigma^2}{n}\right). $$

The truncation at zero reflects the constraint that a variance cannot be negative; it does not establish that the true between-cluster variance is zero. Under the stated normal model, the cluster means are independent and normally distributed. After estimating their common mean, their sample variance has \(J-1\) degrees of freedom:

$$ \frac{(J-1)S^2_{\bar Y}}{\tau^2+\sigma^2/n}\sim\chi^2_{J-1}. $$

This result makes the limitation explicit: the uncertainty in the cluster-mean variance is governed by \(J-1\), not by \(Jn-1\). Increasing the number of individuals within each cluster reduces \(\sigma^2/n\), but it does not create more independent clusters from which to learn about \(\tau^2\).

Worked Example: Estimating Between-Cluster Variation

Suppose a hypothetical study has \(J=5\) clusters with \(n=4\) individuals per cluster. Their observed outcome means are 8, 9, 10, 11, and 12 units. Assume the within-cluster variance is known to be \(\sigma^2=4\) square units and the random-intercept model above is appropriate.

The mean of the five cluster means is \(\bar{\bar Y}=10\) units. Their sample variance is \(S^2_{\bar Y}=[(8-10)^2+(9-10)^2+(10-10)^2+(11-10)^2+(12-10)^2]/(5-1)=10/4=2.5\) square units.

The sampling variance of a cluster mean is \(\sigma^2/n=4/4=1\) square unit. Therefore, \(\hat{\tau}^2_{\mathrm{mom}}=\max(0,2.5-1)=1.5\) square units. The estimated between-cluster standard deviation is \(\sqrt{1.5}=1.225\) units.

This is a point estimate from only five clusters, not a precise measurement of the variance component. In this simplified model, the chi-square result has \(J-1=4\) degrees of freedom, so the sampling variability of \(S^2_{\bar Y}\) is substantial. The calculation also depends on the assumptions that clusters are independent, the model is normal with a common within-cluster variance, and the supplied \(\sigma^2\) is correct.

Worked Example: A Boundary Estimate Is Not Proof of No Clustering

Consider another hypothetical balanced data set with five cluster means: 9.6, 9.8, 10.0, 10.2, and 10.4 units. Suppose \(n=4\) and \(\sigma^2=2\) square units. Assume again that the clusters are independent and the random-intercept model is suitable.

The mean of these cluster means is 10.0 units. Their sample variance is \(S^2_{\bar Y}=[0.16+0.04+0+0.04+0.16]/4=0.10\) square units. The estimated sampling variance of each cluster mean is \(\sigma^2/n=2/4=0.50\) square units.

The untruncated subtraction is \(0.10-0.50=-0.40\) square units, so the constrained method-of-moments estimate is \(\hat{\tau}^2_{\mathrm{mom}}=0\). This boundary result means that the observed spread in these five means is not large enough, under this calculation, to yield a positive estimate after subtracting within-cluster sampling variation. It does not prove that the population variance is exactly zero. With so few clusters, sampling variability can readily obscure real between-cluster differences.

Practical Ways to Stabilize the Imputation Model

Stabilization begins by matching model complexity to what the clusters can support. The aim is not to force a preferred answer, but to avoid estimating more weakly informed parameters than are needed for the imputation target. The earlier material on matching imputation to the multilevel estimand is relevant: simplification must preserve the relationships needed for the intended analysis.

1
Retain essential hierarchy.
When clustering is part of the data structure and the estimand requires it, retain a defensible representation of within- and between-cluster variation. Do not remove the random intercept merely because its estimate is inconvenient or imprecise.
2
Reduce unsupported complexity.
Use a parsimonious conditional model: avoid unnecessary predictors, interactions, or variance parameters, especially terms whose information is largely at the cluster level. Keep variables needed to preserve the analysis relationships and missingness assumptions.
3
Consider regularizing prior information.
In a Bayesian multilevel imputation model, a calibrated prior on a standard deviation, such as a half-normal or half-\(t\) distribution, can discourage implausibly large estimates while allowing uncertainty. Choose its scale using context and assess sensitivity to reasonable alternatives.
4
Check whether the stabilization is consequential.
Compare plausible model specifications and prior scales, inspect parameter and imputation-chain behavior, and check whether conclusions for the intended analysis materially change.

A prior is not a substitute for clusters. With limited cluster-level information, the posterior for a variance component may depend appreciably on the prior. Reporting a single regularized estimate without describing its sensitivity can conceal that dependence. Conversely, using an extremely diffuse prior is not automatically safer: it may allow poorly identified variance values and slow or unstable computation.

Worked Example: How a Prior Can Matter With Few Clusters

Suppose four hypothetical cluster-mean residuals, measured relative to a known intercept, are \(1,-1,0,0\) units. Assume \(n=4\) and \(\sigma^2=4\) square units, so the sampling variance of each cluster mean is \(v=\sigma^2/n=1\) square unit. Treat the residuals as independent normal observations with variance \(v+\tau^2\), and put a half-normal prior with scale 1 unit on \(\tau\): its density is proportional to \(\exp(-\tau^2/2)\) for \(\tau\geq 0\). These are deliberately simplifying assumptions for illustrating prior influence, not a general recipe for choosing a prior.

The residual sum of squares is \(Q=1^2+(-1)^2+0^2+0^2=2\) square units. Apart from factors common to all values of \(\tau\), the likelihood is \(L(\tau)\propto(1+\tau^2)^{-4/2}\exp[-2/\{2(1+\tau^2)\}]\).

Compare \(\tau=0.5\) units with \(\tau=1\) unit. The likelihood ratio is \(L(0.5)/L(1)=[1.25^{-2}\exp(-2/2.5)]/[2^{-2}\exp(-2/4)]\). Numerically, the numerator is approximately \(0.64(0.4493)=0.2876\), and the denominator is approximately \(0.25(0.6065)=0.1516\), giving a likelihood ratio of about \(1.897\).

The prior density ratio is \(\pi(0.5)/\pi(1)=\exp[-0.5^2/2]/\exp[-1^2/2]=\exp(0.375)=1.455\). Thus the posterior density ratio at these two points is approximately \(1.897(1.455)=2.760\). This calculation shows that both the data and the prior affect the relative support for variance values. It does not by itself give a posterior mean, interval, or probability; those require integrating and normalizing the full posterior density.

Diagnostics, Alternatives, and Limitations

No single diagnostic proves that a multilevel imputation model is stable. Use several checks, keeping the small number of independent clusters in view:

  • Inspect variance-component behavior. Repeated estimates at zero, near-boundary estimates, very large changes after modest model edits, or convergence warnings are signs that the data may weakly identify a component.
  • Examine chain traces and mixing. In chained equations, inspect whether relevant parameter traces move across their plausible ranges and whether chains show persistent drift or long periods at a boundary. A stable-looking trace does not establish model adequacy.
  • Check sensitivity. Compare defensible parsimonious specifications and, where applicable, reasonable prior scales. Large changes in imputed distributions or substantive estimates indicate that conclusions rely on modeling choices.
  • Review the cluster-level information. Count independent clusters, examine whether a small number dominate the information, and check whether the cluster-level predictors are supported across the observed clusters. A large individual-level sample does not resolve a lack of cluster-level variation.

The main alternatives have different costs. A parsimonious random-intercept model retains a model for clustering but may omit important structure if simplified too far. Regularizing priors retain a probabilistic model and can stabilize weakly identified parameters, but make conclusions more dependent on prior choices. Treating cluster effects as fixed effects avoids estimating a random-effect variance in that particular formulation, but changes what can be estimated and can be unsuitable for cluster-level predictors or for imputing new clusters. A single-level model is simpler but may fail to preserve dependence and the multilevel relationships needed for the analysis. None is a universal remedy; the choice depends on the estimand and the data-generating structure.

A common mistake is to regard a zero variance estimate as evidence that clustering can be ignored, or to add a strong prior until the computation looks stable. The first confuses weak evidence with evidence of absence; the second can conceal how much the result is being driven by assumptions. Stabilization should be justified by the model and target, and its influence should be made visible through sensitivity checks.

Key takeaway: With few clusters, between-cluster variance and related random-effect quantities are often weakly estimated. Parsimony and calibrated regularization can improve numerical and statistical stability, but neither creates information. Preserve the multilevel structure needed for the estimand, state assumptions, and report sensitivity to modeling choices.

Check Your Understanding

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

  1. Why can 1,000 individuals in a small number of clusters still provide weak information about the between-cluster variance?
  2. In the balanced random-intercept model, what two variance components make up \(\operatorname{Var}(\bar Y_j)\)?
  3. A method-of-moments estimate of a variance component is truncated to zero. What does that establish, and what does it not establish?
  4. Name one benefit and one limitation of using a regularizing prior for the random-intercept standard deviation.
  5. Why is a single-level imputation model not automatically an appropriate solution to instability in a multilevel model?

Solutions. (1) Information about between-cluster variation comes from independent clusters; adding individuals within existing clusters mainly improves the precision of each cluster mean. (2) The variance is \(\tau^2+\sigma^2/n\): between-cluster variation plus within-cluster sampling variation in the mean. (3) It establishes that the constrained calculation returned its boundary value because the subtraction was nonpositive; it does not prove that the population variance is zero. (4) A calibrated prior can regularize weakly identified values; with few clusters, however, the posterior may depend substantially on the prior. (5) A single-level model may fail to represent dependence or preserve relationships needed by the analysis, so it can trade an estimation problem for model misspecification.

Clinical Biostats