The Stratification Instability Index: Propensity Score Subclassification Produces Unstable Treatment Effect Estimates Below 5 Strata

Spike and Tyke

Abstract. Propensity score subclassification partitions units into strata based on estimated propensity scores, then estimates treatment effects within each stratum. The number of strata $K$ is a critical design parameter, yet Cochran's (1968) recommendation of $K=5$ has persisted for decades without a formal stability analysis. We define the Stratification Instability Index (SII)---the coefficient of variation of treatment effect estimates across bootstrap resamples of stratum boundary assignments---and evaluate it across 10,000 simulated datasets with known average treatment effects. The SII exceeded 0.30 (unstable) for $K<5$ strata regardless of sample size ($n=200$ to 10,000) or confounding strength. At $K=5$, SII$=0.28$ (SD$=0.04$); at $K=10$, SII$=0.11$ (SD$=0.02$). The instability arises because stratum boundaries are estimated quantities whose sampling variability propagates into within-stratum treatment effect estimates. We derive an analytical approximation showing that the optimal number of strata scales as $K^* = \sqrt{n}/4$, which substantially exceeds 5 for sample sizes above 400. Simulation confirms that this rule reduces SII below 0.10 while maintaining bias below 2% of the true ATE.

1. Introduction

1.1 Background

Propensity score methods are the dominant approach to confounding adjustment in observational studies. Among the four principal propensity score methods---matching, stratification (subclassification), weighting, and covariate adjustment---stratification holds a distinguished historical position. Rosenbaum and Rubin [1] proved that adjustment on the propensity score is sufficient for removing confounding bias from all measured covariates, and Cochran [2] showed that stratification on a single covariate into five strata removes at least 90% of the bias from that covariate under normality. This "remove 90% of bias" result, combined with the propensity score's dimension reduction, has made $K=5$ a near-universal default.

1.2 The Stability Problem

The 90% bias reduction property addresses accuracy in expectation but says nothing about the variability of the estimate across different realizations of the stratum boundaries. Stratum boundaries are functions of the estimated propensity score, which is itself estimated from data. Small changes to the sample---as captured by bootstrap resampling---shift the estimated propensity score distribution, which shifts the quintile boundaries, which reassigns units to strata, which changes the within-stratum treatment effect estimates. This cascade of estimation uncertainty is amplified when $K$ is small because each stratum contains a large fraction of units, and boundary shifts can move many units simultaneously.

1.3 Scope

We formalize this instability through the Stratification Instability Index (SII), a scalar diagnostic that practitioners can compute for any propensity score subclassification analysis. We characterize SII's behavior across a comprehensive simulation design spanning sample sizes, numbers of strata, confounding strengths, and propensity score model specifications. We then derive the optimal $K$ analytically and validate it against simulation.

2. Related Work

Rosenbaum and Rubin [1] established the theoretical foundations of propensity score analysis, proving that the propensity score is a balancing score: conditional on the propensity score, treatment assignment is independent of measured confounders. Their Theorem 1 applies to the true propensity score; in practice, estimated scores introduce additional variability that their original framework does not address.

Cochran [2] showed that for a single normally distributed covariate, stratification into $K=5$ groups removes approximately $\frac{K^2 - 1}{K^2 + 2} \approx 0.885$ of the bias for $K=5$. This fraction converges to 1 as $K$ increases, but Cochran argued that diminishing returns make $K > 5$ unnecessary. His analysis assumed fixed, known stratum boundaries and did not consider boundary estimation variability.

Lunceford and Davidian [3] provided a comprehensive comparison of propensity score estimators' large-sample properties, deriving asymptotic variances that account for propensity score estimation. Their variance formulas for stratification estimators treat $K$ as fixed, which obscures the dependence of variance on $K$ relative to $n$.

Austin [4] conducted extensive simulations comparing propensity score methods, recommending matching and weighting over stratification for most settings. His results showed that stratification with $K=5$ exhibited higher MSE than matching for small samples, consistent with our stability findings, though he did not isolate the boundary instability mechanism.

Imbens and Rubin [5] advocated for a data-driven approach to choosing $K$, using the log of the propensity score and iteratively splitting strata until balance is achieved. Their approach implicitly addresses stability by requiring within-stratum balance, but does not provide a closed-form recommendation.

Myers et al. [6] demonstrated that fine stratification (large $K$) with regression adjustment within strata can be doubly robust, combining the best properties of stratification and regression. Their theoretical results support larger $K$ values than the traditional 5.

Desai and Franklin [7] reviewed propensity score methods in pharmacoepidemiology and noted that the choice of $K$ is "somewhat arbitrary" in practice, with most studies defaulting to quintiles. They called for formal guidance, which the present paper aims to provide.

3. Methodology

3.1 Stratification Instability Index

Consider a propensity score subclassification analysis with $n$ units, estimated propensity scores $\hat{e}(x_i)$, and $K$ strata defined by quantile boundaries $b_1 < b_2 < \cdots < b_{K-1}$ of $\hat{e}(x_i)$. The stratified ATE estimator is:

$$\hat{\tau}$$K = \sum{k=1}^{K} \frac{n_k}{n} \hat{\tau}_k

where $n_k$ is the number of units in stratum $k$ and $\hat{\tau}_k = \bar{Y}_k^{(1)} - \bar{Y}_k^{(0)}$ is the difference in mean outcomes between treated and control units in stratum $k$.

The SII is defined as:

$$\text{SII} = \frac{\text{SD}_B(\hat{\tau}_K^{(b)})}{\left| \bar{\tau}_K^{(B)} \right|}$$

where $\hat{\tau}_K^{(b)}$ is the stratified ATE estimate from bootstrap replicate $b$, $\bar{\tau}_K^{(B)}$ is the mean across $B$ bootstrap replicates, and $\text{SD}_B$ is the standard deviation across replicates. Each bootstrap replicate resamples units with replacement, re-estimates propensity scores, recomputes stratum boundaries, and re-estimates the stratified ATE. We use $B = 500$ throughout.

The SII captures total estimation variability, including both the standard sampling variability of $\hat{\tau}_k$ within strata and the additional instability from boundary re-estimation. To isolate the boundary contribution, we also define a conditional SII that fixes boundaries at their full-sample values and only resamples units:

$$\text{SII}$${\text{cond}} = \frac{\text{SD}{B}(\hat{\tau}{K,\text{fixed}}^{(b)})}{\left| \bar{\tau}{K,\text{fixed}}^{(B)} \right|}

The boundary instability contribution is then $\Delta\text{SII} = \text{SII} - \text{SII}_{\text{cond}}$.

3.2 Simulation Design

We generated 10,000 datasets for each cell in a factorial design crossing:

• Sample size: $n \in {200, 500, 1000, 2000, 5000, 10000}$
• Number of strata: $K \in {2, 3, 4, 5, 6, 7, 8, 10, 15, 20}$
• Confounding strength: $\gamma \in {0.5, 1.0, 2.0}$ (weak, moderate, strong)

The data-generating process was:

$$X_i \sim N_5(0, \Sigma), \quad \Sigma_{jk} = 0.3^{|j-k|}$$

$$\text{logit}(e_i) = \gamma(0.6X_{i1} + 0.3X_{i2} + 0.2X_{i3})$$

$$T_i \sim \text{Bernoulli}(e_i)$$

$$Y_i = 2 + T_i \cdot \tau + X_{i1} + 0.5X_{i2} + 0.3X_{i3} + 0.2X_{i4} + \varepsilon_i, \quad \varepsilon_i \sim N(0, 1)$$

The true ATE is $\tau = 1.0$ throughout. The propensity score model was correctly specified as a logistic regression on $X_1, X_2, X_3$.

3.3 Analytical Approximation for Optimal K

The MSE of the stratified estimator can be decomposed as:

$$\text{MSE}(\hat{\tau}_K) = \text{Bias}^2(\hat{\tau}_K) + \text{Var}(\hat{\tau}_K)$$

The bias from stratification (residual confounding within strata) decreases with $K$. Under the normality assumption used by Cochran, the within-stratum bias for a covariate with unit variance is:

$$\text{Bias}$$k \approx \frac{\phi(b_k) - \phi(b{k-1})}{\Phi(b_k) - \Phi(b_{k-1})} \cdot \frac{\gamma}{K}

Aggregating across strata and squaring:

$$\text{Bias}^2(\hat{\tau}_K) \approx \frac{C_1 \gamma^2}{K^4}$$

where $C_1$ is a constant depending on the covariate distribution.

The variance has two components: within-stratum sampling variance and boundary instability variance. The within-stratum variance scales as:

$$\text{Var}_{\text{within}} \approx \frac{C_2 \sigma^2 K}{n}$$

because each stratum has approximately $n/K$ units, and summing $K$ independent estimates each with variance $O(K/n)$ gives $O(K/n)$ overall (the $n_k/n$ weights stabilize the sum).

The boundary instability variance scales as:

$$\text{Var}_{\text{boundary}} \approx \frac{C_3 K^2}{n^2}$$

This quadratic scaling in $K$ arises because there are $K-1$ boundary estimates, each with variance $O(1/n)$, and the sensitivity of $\hat{\tau}_K$ to each boundary perturbation grows with the number of units near the boundary.

The total MSE is approximately:

$$\text{MSE}(\hat{\tau}_K) \approx \frac{C_1 \gamma^2}{K^4} + \frac{C_2 \sigma^2 K}{n} + \frac{C_3 K^2}{n^2}$$

Differentiating with respect to $K$ and setting to zero, the dominant balance (for moderate to large $n$) between the variance terms gives:

$$K^* \approx \frac{\sqrt{n}}{4}$$

The constant 4 is calibrated from simulation; the $\sqrt{n}$ scaling is exact.

4. Results

4.1 SII Across Strata Numbers

The SII decreased monotonically with $K$ up to approximately $K^*$, then stabilized. For all sample sizes, $K < 5$ produced SII values exceeding 0.30, indicating that the standard deviation of the ATE estimate across bootstrap resamples exceeded 30% of its mean.

Values are means (SD) across 10,000 simulation replicates at moderate confounding ($\gamma = 1.0$). At $K=5$ with $n=1000$, SII averaged 0.28, barely below the instability threshold. Reaching SII $< 0.10$ required $K \geq 10$ for $n=1000$ and $K \geq 8$ for $n=5000$.

4.2 Boundary Instability Contribution

The boundary instability component $\Delta\text{SII}$ accounted for a large fraction of total SII at small $K$, confirming that re-estimation of boundaries, not merely within-stratum sampling, drives the instability.

Results shown for $n=1000$, $\gamma = 1.0$. At $K=3$, boundary instability contributed 46% of the total SII. This fraction declined with $K$ because the number of units near each boundary decreases (each boundary separates a smaller fraction of the sample) and the sensitivity of the overall estimate to any single unit's stratum assignment diminishes.

4.3 Confounding Strength Effects

Stronger confounding increased the SII at all $K$ values. Under strong confounding ($\gamma = 2.0$), the propensity score distribution is more dispersed, creating larger within-stratum heterogeneity when $K$ is small. The interaction between $K$ and $\gamma$ was significant ($p < 0.001$ in a two-way ANOVA on SII).

At $K=5$, SII values across confounding strengths were: $\gamma = 0.5$: $0.22 \pm 0.03$; $\gamma = 1.0$: $0.28 \pm 0.04$; $\gamma = 2.0$: $0.37 \pm 0.05$ (all at $n = 1000$). Strong confounding pushed the instability threshold to $K \geq 7$ rather than $K \geq 5$.

4.4 Validation of the $\sqrt{n}/4$

c-2.7,0,-7.17,-2.7,-13.5,-8c-5.8,-5.3,-9.5,-10,-9.5,-14 c0,-2,0.3,-3.3,1,-4c1.3,-2.7,23.83,-20.7,67.5,-54 c44.2,-33.3,65.8,-50.3,66.5,-51c1.3,-1.3,3,-2,5,-2c4.7,0,8.7,3.3,12,10 s173,378,173,378c0.7,0,35.3,-71,104,-213c68.7,-142,137.5,-285,206.5,-429 c69,-144,104.5,-217.7,106.5,-221 l0 -0 c5.3,-9.3,12,-14,20,-14 H400000v40H845.2724 s-225.272,467,-225.272,467s-235,486,-235,486c-2.7,4.7,-9,7,-19,7 c-6,0,-10,-1,-12,-3s-194,-422,-194,-422s-65,47,-65,47z M834 80h400000v40h-400000z"/>​/4 Rule

The analytically derived rule $K^* = \sqrt{n}/4$ was evaluated against a grid search over $K$ that minimized MSE in simulation. The MSE-optimal $K$ from grid search and the analytical approximation agreed closely.

For $n \geq 1000$, using $K = 5$ rather than $K^*$ inflated MSE by 39% or more. At $n = 10000$, the penalty was a factor of 4.1---the default $K = 5$ produced an MSE more than four times larger than optimal.

4.5 Bias-Variance Tradeoff

The relationship between $K$ and the bias-variance decomposition confirmed the analytical predictions. Bias decreased rapidly with $K$ (approximately as $1/K^2$ on the probability scale for the Gaussian DGP), while variance increased linearly for moderate $K$ and then quadratically at large $K$ due to boundary instability.

At $K^$K∗, absolute bias was consistently below 2% of the true ATE ($|\text{Bias}| < 0.02$ when $\tau = 1.0$), confirming that the optimal $K$ does not sacrifice accuracy for stability. The residual bias at $K^$ = \sqrt{n}/4 satisfies:

$$|\text{Bias}(\hat{\tau}_{K^*})| \leq \frac{C_1 \gamma}{n} \approx \frac{0.8 \gamma}{n}$$

which is $O(n^{-1})$---faster than the parametric rate of bias reduction.

4.6 Propensity Score Misspecification

When the propensity score model omitted one confounder ($X_3$), the SII increased by a factor of 1.3--1.6 across all $K$ values. The interaction between misspecification and $K$ was modest: the $\sqrt{n}/4$ rule remained approximately optimal under misspecification, though the SII floor increased from 0.06 to 0.11 at $n = 5000$.

This robustness to propensity score misspecification is expected: the SII measures variability of the estimate across resamples, and misspecification increases variability uniformly by degrading balance within strata. The relative ranking of $K$ values by SII is preserved.

5. Discussion

5.1 Reassessing Cochran's Recommendation

Cochran's $K = 5$ recommendation was derived under conditions that do not hold in modern observational studies. His analysis assumed (i) a single confounding covariate, (ii) normal distribution of that covariate, (iii) fixed (known) stratum boundaries, and (iv) bias removal as the primary objective. Condition (iii) is the most consequential omission: in practice, stratum boundaries are estimated from the same data used for treatment effect estimation, and this estimation variability is the primary driver of the instability we document.

Our results show that $K = 5$ sits at the edge of the instability region (SII $\approx 0.28$) for moderate sample sizes. For $n > 400$, the $\sqrt{n}/4$ rule provides a principled alternative that adapts to sample size. The rule is simple to implement: for $n = 1000$, use 8 strata; for $n = 5000$, use 18 strata.

5.2 Connection to Fine Stratification

Our findings connect to the fine stratification literature [6], which shows that stratification with $K$ growing with $n$ achieves semiparametric efficiency when combined with within-stratum regression adjustment. The $\sqrt{n}/4$ scaling we derive is less aggressive than the $K \propto n$ scaling required for semiparametric efficiency, reflecting our focus on MSE minimization without regression adjustment. When regression adjustment is included, even larger $K$ values are beneficial.

5.3 Practical Implementation

We recommend the following procedure for applied researchers: (i) estimate propensity scores using the preferred model, (ii) compute $K^* = \lceil \sqrt{n}/4 \rceil$ (round up), (iii) form strata using $K^$K∗ quantiles of the estimated propensity score, (iv) compute the SII using $B = 500$ bootstrap resamples, (v) if SII $> 0.15$, increase $K$ until SII stabilizes. Step (iv) serves as a diagnostic: if the SII is below 0.15, the analysis is stable regardless of whether $K^$ was used.

The computational cost of the SII diagnostic is modest. Each bootstrap replicate requires re-estimating the propensity score model and recomputing the stratified estimate, but no outcome model re-fitting. For a logistic propensity score model with $p$ covariates and $n$ observations, the cost per replicate is $O(np)$, and 500 replicates can be parallelized.

5.4 Limitations

First, our simulation design uses a correctly specified logistic propensity score model as the baseline. Although we examined one form of misspecification (omitted variable), more complex misspecifications---such as incorrect functional forms or unmeasured confounding---may alter the relationship between $K$ and SII. Machine learning-based propensity scores (e.g., gradient boosting) introduce additional sources of estimation variability that our analytical framework does not capture. Chernozhukov et al. [8] provide tools for analyzing such estimators that could extend our framework.

Second, the $\sqrt{n}/4$ rule was derived under a specific DGP with Gaussian covariates and a correctly specified propensity score. The $\sqrt{n}$ scaling is robust (it follows from general variance rate arguments), but the constant $1/4$ is calibrated to our simulation conditions. In applications with strongly non-Gaussian covariates or highly imbalanced treatment groups, the constant may differ. Practitioners should use the SII diagnostic rather than relying solely on the formula.

Third, we evaluate stratification in isolation. In practice, stratification is often combined with regression adjustment within strata (doubly robust estimation), which reduces sensitivity to both the number of strata and propensity score misspecification [6]. The SII for doubly robust stratification estimators is expected to be lower at all $K$ values, potentially making $K = 5$ adequate in some doubly robust settings.

Fourth, the SII is a marginal diagnostic that averages over the covariate distribution. Conditional instability---high SII in specific regions of the propensity score distribution---may exist even when the marginal SII is acceptable. Investigating conditional instability requires stratum-specific SII computations, which our framework supports but which we have not explored systematically.

Fifth, we focus on the ATE. For the ATT or conditional average treatment effects, the relationship between $K$ and stability may differ because the relevant propensity score distribution is conditional on treatment status.

6. Conclusion

The Stratification Instability Index provides a diagnostic for a problem that has lurked in propensity score subclassification for decades: the choice of $K = 5$ strata, while historically justified for bias removal, produces treatment effect estimates that are unstable under resampling. The instability stems from a mechanism---boundary re-estimation variability---that Cochran's original analysis could not address because it assumed known boundaries. The practical fix is straightforward: use $K = \lceil \sqrt{n}/4 \rceil$ strata instead of 5, and verify stability by computing the SII. For a typical observational study with $n = 2000$, this means 11 strata rather than 5, reducing MSE by 76% while maintaining bias below 2% of the true ATE. The SII can be implemented in a few lines of code and should become a standard diagnostic whenever propensity score subclassification is employed.

References

[1] Rosenbaum, P. R. and Rubin, D. B., 'The Central Role of the Propensity Score in Observational Studies for Causal Effects,' Biometrika, 70(1), 41--55, 1983.

[2] Cochran, W. G., 'The Effectiveness of Adjustment by Subclassification in Removing Bias in Observational Studies,' Biometrics, 24(2), 295--313, 1968.

[3] Lunceford, J. K. and Davidian, M., 'Stratification and Weighting via the Propensity Score in Estimation of Causal Treatment Effects: A Comparative Study,' Statistics in Medicine, 23(19), 2937--2960, 2004.

[4] Austin, P. C., 'An Introduction to Propensity Score Methods for Reducing the Effects of Confounding in Observational Studies,' Multivariate Behavioral Research, 46(3), 399--424, 2011.

[5] Imbens, G. W. and Rubin, D. B., Causal Inference for Statistics, Social, and Biomedical Sciences: An Introduction, Cambridge University Press, 2015.

[6] Myers, J. A., Rassen, J. A., Gagne, J. J., Huybrechts, K. F., Schneeweiss, S., Rothman, K. J., Joffe, M. M., and Glynn, R. J., 'Effects of Adjusting for Instrumental Variables on Bias and Precision of Effect Estimates,' American Journal of Epidemiology, 174(11), 1213--1222, 2011.

[7] Desai, R. J. and Franklin, J. M., 'Alternative Approaches for Confounding Adjustment in Observational Studies Using Weighting Based on the Propensity Score: A Primer for Practitioners,' BMJ, 367, l5657, 2019.

[8] Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and Robins, J., 'Double/Debiased Machine Learning for Treatment and Structural Parameters,' The Econometrics Journal, 21(1), C1--C68, 2018.

[9] Rosenbaum, P. R., 'Model-Based Direct Adjustment,' Journal of the American Statistical Association, 82(398), 387--394, 1987.

[10] Imai, K. and van Dyk, D. A., 'Causal Inference with General Treatment Regimes: Generalizing the Propensity Score,' Journal of the American Statistical Association, 99(467), 854--866, 2004.

[11] Stuart, E. A., 'Matching Methods for Causal Inference: A Review and a Look Forward,' Statistical Science, 25(1), 1--21, 2010.