Bayesian Estimation of Variance under Fine Stratification via Mean-Variance Smoothing

1. Problem Statement

Fine stratification is a survey design where the population is divided into a large number of disjoint strata, with very few Primary Sampling Units (PSUs) selected per stratum (often just one or two). While this design yields unbiased point estimators (e.g., the Horvitz-Thompson estimator) for population means, it creates a significant challenge for variance estimation:

• One PSU per stratum: The standard design-based variance estimator is undefined because the second-order inclusion probability ($\pi_{jk}$) is zero.
• Two or Three PSUs per stratum: The variance estimator exists but suffers from high variability and instability.

Current Limitations:
The traditional solution is collapsed strata, where adjacent strata are paired to create "pseudo-strata" with more units, allowing variance calculation. However, this method introduces positive bias in the variance estimator. The bias increases as the population means of the paired strata differ significantly, leading to unnecessarily wide confidence intervals. Existing non-parametric alternatives (e.g., Breidt et al., 2016) still rely on some form of collapsing or kernel weighting that can be sensitive to boundary issues and bandwidth selection.

2. Methodology

The authors propose a Bayesian approach that avoids collapsing strata entirely. Instead, they model the mean and variance functions simultaneously across all strata using penalized splines (P-splines).

Core Model

The method assumes that the stratum means and variances vary smoothly with respect to an auxiliary variable (collapsing index) $x_h$.

• One PSU per Stratum ($n_h=1$):
  The observation $y_h$ in stratum $h$ is modeled as:
  $$y_h \sim N(m(x_h), s^2(x_h))$$
  where $m(x_h)$ is the mean function and $s^2(x_h)$ is the variance function.
• Multiple PSUs per Stratum ($n_h > 1$):
  The model extends to a multivariate normal distribution with an equicorrelation structure:
  $$\mathbf{y}_h \sim N(m(x_h)\mathbf{1}_{n_h}, s^2(x_h)\mathbf{R}_h(\rho))$$
  where $\mathbf{R}_h(\rho)$ is a correlation matrix with parameter $\rho$.

Smoothing via Penalized Splines

Both the mean function $m(x)$ and the log-variance function $\log s^2(x)$ are modeled using P-splines to ensure flexibility and prevent overfitting:
$$m(x) = \beta_0 + \beta_1 x + \dots + \beta_q x^q + \sum_{l=1}^L \beta_{q+l}(x - \kappa_l)_+^q$$
$$\log s^2(x) = \gamma_0 + \gamma_1 x + \dots + \gamma_q x^q + \sum_{l=1}^L \gamma_{q+l}(x - \kappa_l)_+^q$$

• Random Effects: The higher-order spline coefficients ($\beta_{q+l}, \gamma_{q+l}$) are treated as random effects with normal priors to enforce smoothness.
• Weighted Likelihood: To account for the complex survey design, the authors incorporate sampling weights ($\pi_h^{-1}$) into the likelihood function, creating a pseudo-likelihood (or weighted likelihood) framework. This ensures the posterior distribution respects the design-based inclusion probabilities.

Estimation

• Inference: The joint posterior distribution of all parameters (regression coefficients, variance parameters, and correlation $\rho$) is approximated using Markov Chain Monte Carlo (MCMC) (specifically Metropolis-within-Gibbs sampling).
• Variance Estimator: The Bayesian variance estimator for the Horvitz-Thompson mean is derived by summing the posterior mean of the smoothed stratum variances, weighted by the design weights:
  $$\widehat{var}_{Bayes} = \frac{1}{N^2} \sum_{h=1}^H w_h^2 \hat{s}^2(x_h)$$

3. Key Contributions

1. Elimination of Strata Collapsing: Unlike previous methods, this approach does not require pairing or collapsing strata, thereby avoiding the positive bias associated with differences in stratum means.
2. Joint Mean-Variance Smoothing: The method simultaneously models the mean and variance functions, leveraging the correlation between them to improve estimation efficiency.
3. Design-Aware Bayesian Framework: It successfully integrates survey design weights (inclusion probabilities) into a Bayesian hierarchical model, providing a rigorous framework for informative sampling.
4. Uncertainty Quantification: Beyond point estimates, the method provides full posterior distributions, allowing for better characterization of uncertainty in the variance estimates themselves.

4. Results

The authors evaluated their method through simulation studies and a real-data application, comparing it against the Collapsed Strata (varColl) and Kernel-based (varKer) estimators.

• Simulation Studies:
  • Gaussian and Gamma Populations: The Bayesian estimator ($var_{Bayes}$) consistently demonstrated lower Absolute Bias (AB) and Root Mean Squared Error (RMSE) compared to competing methods, particularly as the number of strata increased or under systematic sampling designs.
  • Coverage Probability: The 95% confidence intervals generated by the Bayesian method were closer to the nominal 95% level than those from collapsed or kernel methods, which often suffered from under-coverage or excessive width.
• Real Data Application (NSFG):
  • Using data from the National Survey of Family Growth (2015–2017), the method was applied to estimate variances for pregnancy age, education, and number of babies.
  • Result: The Bayesian estimator produced significantly smaller variance estimates and shorter confidence intervals compared to the collapsed and kernel methods, without sacrificing coverage accuracy. This suggests the traditional methods were overestimating variance due to the bias of collapsing.

5. Significance

This paper offers a robust alternative to the long-standing practice of collapsing strata in fine stratification surveys. By utilizing penalized splines within a Bayesian framework, the authors provide a method that:

• Reduces Bias: Eliminates the systematic overestimation of variance caused by pairing dissimilar strata.
• Improves Precision: Yields tighter confidence intervals, leading to more powerful statistical inferences.
• Enhances Flexibility: Can be easily adapted to multiple PSUs per stratum and extended to multivariate settings or different correlation structures.
• Practical Utility: The methodology is computationally feasible (via MCMC) and directly applicable to major government surveys (like the US Census Bureau's or NCHS surveys) that utilize fine stratification designs.

The work represents a significant step forward in survey methodology, bridging the gap between design-based inference and model-based smoothing to solve the "variance estimation problem" in sparse sampling designs.