Fay-Herriot Model with Spectral Clustering

• FH-SC is a fully Bayesian small area estimation method that integrates spectral clustering-derived priors into the classical Fay-Herriot framework to capture complex covariate-driven heterogeneity.
• It employs a three-level hierarchical model with spectral clustering to form data-driven clusters, enabling improved precision via Rao-Blackwellized estimates and rigorous uncertainty quantification.
• The approach supports benchmarked estimation and introduces the novel Conditional Posterior Mean Square Error (CPMSE) metric, demonstrating significant gains over traditional SAE methods.

The Fay-Herriot Model with Spectral Clustering (FH-SC) is a fully Bayesian methodology for small area estimation (SAE) that integrates spectral clustering-derived random-effect priors into the classical Fay-Herriot (FH) model framework. Unlike traditional spatial or geographic-based SAE models, FH-SC leverages external covariates to induce data-driven clusters, enhances precision by borrowing strength within clusters of similar areas, and supports rigorous benchmarking and uncertainty quantification, including closed-form Rao-Blackwellized estimators and a novel Conditional Posterior Mean Square Error (CPMSE) metric (Fúquene-Patiño, 17 Dec 2025).

1. Hierarchical Model Specification and Priors

The FH-SC model posits the following three-level hierarchical structure for $m$ small areas partitioned into $C$ clusters via spectral clustering:

• Sampling Level: For each cluster $c$ (with $n_c$ areas),

$y_c \mid \theta_c, D_c \sim N_{n_c}(\theta_c, D_c)$

where $\mathbf{y}_c = (y_{1,c}, ..., y_{n_c,c})'$ denotes direct survey estimates with known sampling variances $D_c = \mathrm{diag}(D_{1,c},...,D_{n_c,c})$, and $\theta_c = (\theta_{1,c},..., \theta_{n_c,c})'$ are the true area effects.

• Linking Level: Linking model for true parameters,

$$\theta_c = A_{\rho,c}^{-1} \eta_c , \quad \eta_c \mid \beta, u_c = X_c \beta_c + Z_c u_c, \quad u_c \sim N_{h_c}(0, G_{\phi,c})$$

where $X_c$ is the design matrix, $C$0 the regression coefficients, $C$1 the matrix for random effects $C$2, and $C$3 their covariance.

The cluster regularization operator

$C$4

uses the cluster Laplacian $C$5 to induce within-cluster smoothness and regularization.

• Priors: Typical prior assignments are
  • $C$6 flat or Normal,
  • $C$7 Inverse-Gamma or Gamma (on precisions),
  • $C$8, with $C$9.

This construction enables flexible clustering effects, with cluster-wise or global $c$0 and $c$1.

2. Spectral Clustering for Cluster Geometry

Prior to model fitting, spectral clustering is performed using external covariates $c$2 (e.g., poverty or educational indices):

1. Similarity Matrix: Construct $c$3 by $c$4.
2. Adjacency Matrix: Build $c$5 (e.g., $c$6-nearest-neighbor or $c$7-threshold) with $c$8 if $c$9, $n_c$0 are neighbors, else 0.
3. Laplacian: Form unnormalized graph Laplacian $n_c$1, with $n_c$2, $n_c$3.
4. Eigenvector Embedding: Extract the first $n_c$4 eigenvectors $n_c$5 of $n_c$6, stack rows into $n_c$7.
5. Clustering: Apply $n_c$8-means to rows of $n_c$9 to assign areas to clusters $y_c \mid \theta_c, D_c \sim N_{n_c}(\theta_c, D_c)$0.
6. Block Laplacian and Regularizer: For each cluster, set $y_c \mid \theta_c, D_c \sim N_{n_c}(\theta_c, D_c)$1, assemble $y_c \mid \theta_c, D_c \sim N_{n_c}(\theta_c, D_c)$2.
7. Final Operator: $y_c \mid \theta_c, D_c \sim N_{n_c}(\theta_c, D_c)$3; $y_c \mid \theta_c, D_c \sim N_{n_c}(\theta_c, D_c)$4.

This procedure results in clusters that capture complex covariate-driven heterogeneity potentially missed by spatial-only approaches.

3. Bayesian Estimation and Rao-Blackwellization

Bayesian inference in FH-SC is performed via Gibbs sampling with Metropolis–Hastings (MH) updates for the cluster penalty parameter $y_c \mid \theta_c, D_c \sim N_{n_c}(\theta_c, D_c)$5, given its nonstandard conditional posterior. The joint posterior is:

$y_c \mid \theta_c, D_c \sim N_{n_c}(\theta_c, D_c)$6

Key updates in each iteration ($y_c \mid \theta_c, D_c \sim N_{n_c}(\theta_c, D_c)$7):

• $y_c \mid \theta_c, D_c \sim N_{n_c}(\theta_c, D_c)$8 is multivariate Normal; mean and variance depend on $y_c \mid \theta_c, D_c \sim N_{n_c}(\theta_c, D_c)$9, $\mathbf{y}_c = (y_{1,c}, ..., y_{n_c,c})'$0, $\mathbf{y}_c = (y_{1,c}, ..., y_{n_c,c})'$1, $\mathbf{y}_c = (y_{1,c}, ..., y_{n_c,c})'$2, $\mathbf{y}_c = (y_{1,c}, ..., y_{n_c,c})'$3.
• $\mathbf{y}_c = (y_{1,c}, ..., y_{n_c,c})'$4 follows a conjugate Gaussian conditional.
• $\mathbf{y}_c = (y_{1,c}, ..., y_{n_c,c})'$5 is Gamma or Inverse-Gamma (if diagonal/identity structures).
• $\mathbf{y}_c = (y_{1,c}, ..., y_{n_c,c})'$6 is updated via MH with a random walk on $\mathbf{y}_c = (y_{1,c}, ..., y_{n_c,c})'$7.

After MCMC sampling, posterior means (ergodic samples) or Rao-Blackwellized (RB) estimates are computed:

$\mathbf{y}_c = (y_{1,c}, ..., y_{n_c,c})'$8

4. Benchmarking Through Posterior Projections

FH-SC supports benchmarked estimation via posterior-projection. Given $\mathbf{y}_c = (y_{1,c}, ..., y_{n_c,c})'$9 linear constraints $D_c = \mathrm{diag}(D_{1,c},...,D_{n_c,c})$0 ($D_c = \mathrm{diag}(D_{1,c},...,D_{n_c,c})$1, $D_c = \mathrm{diag}(D_{1,c},...,D_{n_c,c})$2), benchmarked area draws are defined as solutions to

$D_c = \mathrm{diag}(D_{1,c},...,D_{n_c,c})$3

with closed-form KKT solution (Proposition 4):

$D_c = \mathrm{diag}(D_{1,c},...,D_{n_c,c})$4

RB-benchmarked estimates are averages of the conditional expectations of $D_c = \mathrm{diag}(D_{1,c},...,D_{n_c,c})$5:

$D_c = \mathrm{diag}(D_{1,c},...,D_{n_c,c})$6

with $D_c = \mathrm{diag}(D_{1,c},...,D_{n_c,c})$7 also available in closed form (Definition 9):

$D_c = \mathrm{diag}(D_{1,c},...,D_{n_c,c})$8

5. Uncertainty Quantification: Conditional Posterior MSE (CPMSE)

FH-SC introduces the Conditional Posterior Mean Square Error (CPMSE) for the RB-benchmarked estimators:

$D_c = \mathrm{diag}(D_{1,c},...,D_{n_c,c})$9

where the latter term is the RB-posterior variance of $\theta_c = (\theta_{1,c},..., \theta_{n_c,c})'$0. Empirically, CPMSE is estimated by averaging over posterior draws:

$\theta_c = (\theta_{1,c},..., \theta_{n_c,c})'$1

plus the squared adjustment from benchmarking.

CPMSE serves as a fully Bayesian, generalizable uncertainty measure, with demonstrated frequentist consistency as $\theta_c = (\theta_{1,c},..., \theta_{n_c,c})'$2.

6. Simulation Evidence and Empirical Performance

Performance of FH-SC has been assessed through model- and data-based simulations and a real-data study on Colombian municipalities.

Summary of Results:

Empirical Setting	Key Findings
Model-based (true FH)	CPMSE closely tracks empirical MSE of benchmarked $\theta_c = (\theta_{1,c},..., \theta_{n_c,c})'$3	$\theta_c = (\theta_{1,c},..., \theta_{n_c,c})'$4CPMSE–MSE$\theta_c = (\theta_{1,c},..., \theta_{n_c,c})'$5 as $\theta_c = (\theta_{1,c},..., \theta_{n_c,c})'$6 increases
Data-based (true FH–SC1)	FH–SC1 yielded uniformly smaller absolute/squared errors than FH, especially as $\theta_c = (\theta_{1,c},..., \theta_{n_c,c})'$7 rises	CPMSE remained an accurate proxy
Colombian municipalities	FH–SC2 (common $\theta_c = (\theta_{1,c},..., \theta_{n_c,c})'$8, cluster-specific $\theta_c = (\theta_{1,c},..., \theta_{n_c,c})'$9, free $$\theta_c = A_{\rho,c}^{-1} \eta_c , \quad \eta_c \mid \beta, u_c = X_c \beta_c + Z_c u_c, \quad u_c \sim N_{h_c}(0, G_{\phi,c})$$0) outperformed six competitors (including FH, two FH–C, three FH–SC variants) in DIC and predictive deviance; realized RB coefficient of variation reductions of $$\theta_c = A_{\rho,c}^{-1} \eta_c , \quad \eta_c \mid \beta, u_c = X_c \beta_c + Z_c u_c, \quad u_c \sim N_{h_c}(0, G_{\phi,c})$$192% (non-benchmarked) and $$\theta_c = A_{\rho,c}^{-1} \eta_c , \quad \eta_c \mid \beta, u_c = X_c \beta_c + Z_c u_c, \quad u_c \sim N_{h_c}(0, G_{\phi,c})$$285% (benchmarked) relative to FH; clusters induced by covariates captured heterogeneity missed by spatial methods

In the Colombian internet access application, external indices (Multidimensional Poverty Index, Educational Index) were used for clustering, resulting in $$\theta_c = A_{\rho,c}^{-1} \eta_c , \quad \eta_c \mid \beta, u_c = X_c \beta_c + Z_c u_c, \quad u_c \sim N_{h_c}(0, G_{\phi,c})$$3 clusters and demonstrating the approach's flexibility and improved precision.

7. Methodological Significance and Extensions

Key features of FH–SC include:

• Use of spectral clustering on non-geographic external covariates to construct Laplacian-smoothness priors for random effects,
• Maintenance of the fully Bayesian paradigm, with closed-form $$\theta_c = A_{\rho,c}^{-1} \eta_c , \quad \eta_c \mid \beta, u_c = X_c \beta_c + Z_c u_c, \quad u_c \sim N_{h_c}(0, G_{\phi,c})$$4-conditionals and MH sampling for $$\theta_c = A_{\rho,c}^{-1} \eta_c , \quad \eta_c \mid \beta, u_c = X_c \beta_c + Z_c u_c, \quad u_c \sim N_{h_c}(0, G_{\phi,c})$$5,
• Closed-form Rao–Blackwellization for both plain and benchmarked estimation,
• Substantial gains in estimation precision (lower coefficient of variation and MSE) over existing Bayesian and frequentist SAE approaches, particularly in settings where traditional spatial clustering fails to capture underlying heterogeneity.

A plausible implication is that FH–SC generalizes seamlessly to other benchmarking contexts and Bayesian SAE estimators where linear constraints and covariate-driven clustering may be beneficial (Fúquene-Patiño, 17 Dec 2025).