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, $\beta_c$ the regression coefficients, $Z_c$ the matrix for random effects $u_c$, and $G_{\phi,c}$ their covariance.

The cluster regularization operator

$A_{\rho,c} = I_{n_c} + \frac{1-\rho}{\rho} L_c$

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

• Priors: Typical prior assignments are
  • $\beta_c \sim$ flat or Normal,
  • $G_{\phi,c} \sim$ Inverse-Gamma or Gamma (on precisions),
  • $\rho \sim \mathrm{Beta}(a, b)$, with $\rho\in(0, 1)$.

This construction enables flexible clustering effects, with cluster-wise or global $\beta$ and $G_{\phi,c}$.

2. Spectral Clustering for Cluster Geometry

Prior to model fitting, spectral clustering is performed using external covariates $x_i^* \in \mathbb{R}^{p^*}$ (e.g., poverty or educational indices):

1. Similarity Matrix: Construct $S$ by $$s_{ij} = \exp \left\{ -\|x_i^* - x_j^*\|^2/(2\sigma_s^2) \right\}$$.
2. Adjacency Matrix: Build $W$ (e.g., $k$-nearest-neighbor or $\epsilon$-threshold) with $W_{ij} = s_{ij}$ if $i$, $j$ are neighbors, else 0.
3. Laplacian: Form unnormalized graph Laplacian $L_u = D_u - W$, with $D_u = \mathrm{diag}(d_i)$, $d_i = \sum_j W_{ij}$.
4. Eigenvector Embedding: Extract the first $C$ eigenvectors $\{v_1,...,v_C\}$ of $L_u$, stack rows into $V \in \mathbb{R}^{m \times C}$.
5. Clustering: Apply $k$-means to rows of $V$ to assign areas to clusters $\{D_1,...,D_C\}$.
6. Block Laplacian and Regularizer: For each cluster, set $$L_c = n_c I_{n_c} - \mathbf{1}_{n_c} \mathbf{1}_{n_c}^\top$$, assemble $L_{SC} = \mathrm{blockdiag}(L_1,...,L_C)$.
7. Final Operator: $A_{\rho,c} = I_{n_c} + \frac{1-\rho}{\rho} L_c$; $$A_\rho = \mathrm{blockdiag}(A_{\rho,1}, ..., A_{\rho,C})$$.

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 $\rho$, given its nonstandard conditional posterior. The joint posterior is:

$$p(\theta, \beta, G_\phi, \rho | y) \propto \prod_{c=1}^C N(y_c; A_{\rho,c}^{-1} \eta_c, D_c) N(\eta_c; X_c\beta_c, Z_c G_{\phi,c} Z_c') \pi(\beta_c)\pi(G_{\phi,c})\pi(\rho)$$

Key updates in each iteration ($\ell$):

• $\theta_c | -$ is multivariate Normal; mean and variance depend on $A_{\rho,c}$, $D_c$, $X_c$, $Z_c$, $G_{\phi,c}$.
• $\beta_c | -$ follows a conjugate Gaussian conditional.
• $G_{\phi,c}|-$ is Gamma or Inverse-Gamma (if diagonal/identity structures).
• $\rho|-$ is updated via MH with a random walk on $\log(\rho)$.

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

$$\hat\theta_j^{RB} = \frac{1}{L-T} \sum_{\ell=T+1}^L E(\theta_j | \beta^{(\ell)}, G_{\phi}^{(\ell)}, \rho^{(\ell)}, y)$$

4. Benchmarking Through Posterior Projections

FH-SC supports benchmarked estimation via posterior-projection. Given $k$ linear constraints $W\theta = p$ ($W\in \mathbb{R}^{k\times m}$, $\operatorname{rank}(W)=k$), benchmarked area draws are defined as solutions to

$$\theta_B^{(\ell)} = \arg\min_{W\tilde\theta = p} (\tilde\theta - \theta^{(\ell)})' A_{\rho^{(\ell)}} (\tilde\theta - \theta^{(\ell)})$$

with closed-form KKT solution (Proposition 4):

$$\theta_B^{(\ell)} = \theta^{(\ell)} + A_{\rho^{(\ell)}}^{-1} W' [W A_{\rho^{(\ell)}}^{-1} W']^{-1}(p - W \theta^{(\ell)})$$

RB-benchmarked estimates are averages of the conditional expectations of $\theta_B^{(\ell)}$:

$$\hat\theta_{B,j}^{RB} = \frac{1}{L-T} \sum_{\ell=T+1}^L E(\theta_{B,j}^{(\ell)} | \cdots)$$

with $E(\theta_B|\cdots)$ also available in closed form (Definition 9):

$$E(\theta_B | ...) = E(\theta | ...) + A_\rho^{-1} W' [W A_\rho^{-1} W']^{-1}(p - W E(\theta | ...))$$

5. Uncertainty Quantification: Conditional Posterior MSE (CPMSE)

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

$$\mathrm{CPMSE}(\hat\theta_{B,j}^{RB}) = E_{post} [ (\hat\theta_{B,j}^{RB} - \theta_j)^2 ] = (\hat\theta_{B,j}^{RB} - \hat\theta_{j}^{RB})^2 + \mathrm{CPMSE}(\hat\theta_j^{RB})$$

where the latter term is the RB-posterior variance of $\theta_j$. Empirically, CPMSE is estimated by averaging over posterior draws:

$$\mathrm{CPMSE} \approx \frac{1}{L-T} \sum_{\ell=T+1}^L \left[ (E(\theta_j|\ell) - \hat\theta_j^{RB})^2 + \mathrm{Var}(\theta_j|\ell) \right ]$$

plus the squared adjustment from benchmarking.

CPMSE serves as a fully Bayesian, generalizable uncertainty measure, with demonstrated frequentist consistency as $m\to\infty$.

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 $E[\theta_j|y]$	$|$CPMSE–MSE$|\rightarrow 0$ as $m$ increases
Data-based (true FH–SC1)	FH–SC1 yielded uniformly smaller absolute/squared errors than FH, especially as $\rho$ rises	CPMSE remained an accurate proxy
Colombian municipalities	FH–SC2 (common $\beta$, cluster-specific $\sigma_c^2$, free $\rho$) outperformed six competitors (including FH, two FH–C, three FH–SC variants) in DIC and predictive deviance; realized RB coefficient of variation reductions of $\sim$92% (non-benchmarked) and $\sim$85% (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 $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$-conditionals and MH sampling for $\rho$,
• 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).