Finite Population Simulation: Validating fastsae against Ground Truth

Introduction: Why Finite Population Simulation?

In Small Area Estimation (SAE), simulation studies serve as the definitive benchmark for validating statistical models before deploying them to official statistics (such as poverty mapping or regional health indicators). There are two fundamental paradigms:

1. Model-Based Simulation: Data are generated directly from the parametric distribution assumed by the model (e.g., drawing \(y_d \sim \text{Beta}(\dots)\)). While useful for verifying algorithmic correctness and parameter recovery, it evaluates the model under its own idealized assumptions.
2. Design-Based Simulation on a Finite Population (Gold Standard): A realistic finite population of discrete individuals or households (\(N = 33{,}511\) across \(D = 60\) domains) is generated. The true area-level quantities \(P_d\) are calculated directly by enumerating all units in each domain: \[P_d = \frac{1}{N_d} \sum_{i=1}^{N_d} y_{di} \quad (\text{Exact Non-Parametric Ground Truth})\] Survey samples of small sizes (\(n_d \ll N_d\), e.g., \(n_d = 20 - 55\)) are drawn repeatedly. Direct estimators and design-based sampling variances are computed from the survey sample and fed into Small Area models.

This design-based simulation rigorously tests whether: - Area-level continuous Beta models (hb_area and tipsae) remain robust when the underlying data consist of discrete individual binary outcomes (\(y_{di} \in \{0, 1\}\)). - Design-based sampling variances \(\widehat{\text{Var}}(\hat{p}_d)\) with finite population corrections are accurately handled by the linking model. - Model credible intervals maintain exact nominal 95% coverage under repeated survey sampling.

Simulation Architecture

The simulation architecture mirrors a national household survey (such as Susenas or EU-SILC):

+-------------------------------------------------------------------------------+
| 1. SYNTHETIC FINITE POPULATION (N = 33,511 units, D = 60 domains)            |
|    - Unit covariates: x_1di ~ N(1.6, 0.6^2), x_2di ~ Unif(0, 2.2)            |
|    - Spatial dependency: Besag ICAR random effects (W adjacency matrix)       |
|    - Unit outcome: y_di ~ Bernoulli(pi_di) in {0, 1}                          |
|    - Exact Ground Truth: P_d = sum(y_di) / N_d (Mean P_d = 27.4%)            |
+-------------------------------------------------------------------------------+
                                      |
                                      v
+-------------------------------------------------------------------------------+
| 2. SURVEY SAMPLING (R = 100 Monte Carlo Replications)                         |
|    - SRS without replacement per domain (n_d in [22, 55], mean n_d = 40)      |
|    - Sampling fraction f_d = n_d / N_d ~ 7.1%                                 |
|    - Direct Estimator: p_hat_d = mean(y_sample_d)                             |
|    - Design Sampling Variance: V_hat_d = (1 - f_d) * p(1 - p) / (n_d - 1)     |
+-------------------------------------------------------------------------------+
                                      |
                   +------------------+------------------+
                   |                                     |
                   v                                     v
+------------------------------------+ +------------------------------------+
| 3. fastsae::hb_area()             | | 4. tipsae::fit_sae()               |
|    - Bayesian INLA (Full Laplace)  | |    - Hamiltonian Monte Carlo (Stan)|
|    - Family: "beta"                | |    - Likelihood: "beta"            |
|    - Spatial: "besag"              | |    - Spatial error: TRUE           |
|    - Strategy: "laplace"           | |    - Execution: ~10-37 sec / fit   |
|    - Calibrated PC-Prior (1, 0.05) | |                                    |
|    - Execution: ~3-4 sec / fit     | |                                    |
+------------------------------------+ +------------------------------------+
                   |                                     |
                   +------------------+------------------+
                                      |
                                      v
+-------------------------------------------------------------------------------+
| 5. EMPIRICAL EVALUATION AGAINST FIXED GROUND TRUTH P_d                        |
|    - Empirical RMSE, Empirical RRMSE (%), and Empirical Relative Bias (%)     |
|    - 95% Credible Interval Coverage Rate across 100 replications              |
|    - Efficiency Gain vs Direct Estimator (MSE Ratio)                          |
+-------------------------------------------------------------------------------+

Mathematical Formulations of Empirical Properties

Across \(R = 100\) independent survey sample replications, the empirical properties for each domain \(d = 1, \dots, D\) are defined as:

1. Empirical Mean and Bias: \[\overline{\hat{\theta}}_d = \frac{1}{R}\sum_{r=1}^R \hat{\theta}_d^{(r)}, \quad \text{Bias}_d = \overline{\hat{\theta}}_d - P_d, \quad \text{RB}_d = \frac{\text{Bias}_d}{P_d} \times 100\%\]

2. Empirical Root Mean Squared Error (RMSE): \[\text{RMSE}_d = \sqrt{\frac{1}{R}\sum_{r=1}^R \left(\hat{\theta}_d^{(r)} - P_d\right)^2}, \quad \text{RRMSE}_d = \frac{\text{RMSE}_d}{P_d} \times 100\%\]

3. Empirical 95% Credible Interval Coverage Rate (CR): \[\text{CR}_d = \frac{1}{R}\sum_{r=1}^R \mathbb{I}\left(P_d \in \left[\text{CI}_{\text{lower}, d}^{(r)}, \text{CI}_{\text{upper}, d}^{(r)}\right]\right) \times 100\%\]

4. Empirical Efficiency Gain vs Direct Estimator: \[\text{Eff}_d = \frac{\text{MSE}_d(\text{Direct})}{\text{MSE}_d(\text{Model})}\]

Empirical Benchmark Results (\(R = 100\) Replications, \(D = 60\))

The table below summarizes the macro-averaged empirical metrics across all 60 domains from the full \(R = 100\) Monte Carlo simulation, evaluating fastsae with Full Laplace approximation (inla_engine(strategy = "laplace")) and calibrated PC-Prior (param = c(1, 0.05)):

(Cells highlighted in green indicate superior performance for each evaluation criterion).

National-Scale Benchmark (\(D = 500\) Domains, \(N \approx 250{,}000\) Individuals, \(R = 25\))

To evaluate scalability at the national level (equivalent to the \(\approx 514\) regency/city breakdown in Indonesia or hundreds of counties in other nations), a national synthetic population of \(N = 250{,}824\) individuals was generated across a \(20 \times 25\) spatial grid of \(D = 500\) contiguous domains. Survey samples of \(n_d \in [20, 50]\) per domain (\(n \approx 17{,}500\) survey units) were drawn across \(R = 25\) Monte Carlo replications.

(Cells highlighted in green indicate superior performance).

Visualizations

# Code to replicate domain-wise RMSE visualization
ggplot(df_rmse, aes(x = DomainIdx, y = RMSE, color = Method)) +
  geom_line() + geom_point() +
  labs(title = "Domain-Specific Empirical RMSE (R = 100 Replications)",
       x = "Domain Index", y = "Empirical RMSE vs True P_d")

How to Replicate

The complete finite population benchmark script is included in the package repository:

# Run 100 Monte Carlo replications across 5 CPU cores
Rscript benchmarks/sim_mc_finite_population.R --R=100 --cores=5

References

• Janicki, R. (2020). Properties of the beta regression model for small area estimation of proportions and application to estimation of poverty rates. Communications in Statistics - Theory and Methods, 49(9), 2264–2284.
• De Nicolò, S., & Gardini, A. (2024). The R Package tipsae: Tools for Mapping Proportions and Indicators on the Unit Interval. Journal of Statistical Software, 108(1), 1–36.
• Rao, J. N. K., & Molina, I. (2015). Small Area Estimation (2nd ed.). John Wiley & Sons.
• Tzavidis, N., et al. (2018). From start to finish: a framework for the production of small area official statistics. Journal of the Royal Statistical Society: Series A, 181(4), 927–979.
• Rue, H., Martino, S., & Chopin, N. (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society: Series B, 71(2), 319–392.