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:

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.
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 keep close to 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: ~12-40 sec / fit   |
|    - Default priors (hb-priors)    | |                                    |
|    - Execution: ~6-9 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:

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\%$
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\%$
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\%$
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 hb_area() with its default engine and priors (inla_engine(); see vignette("hb-priors")) and tipsae::fit_sae() with type_disp = "var" (one chain of 600 iterations). Runtimes are per fit, with five replicates running in parallel.

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

Estimator / Model	Empirical RMSE	Empirical RRMSE (%)	Empirical MARB (%)	Empirical RB (%)	Empirical 95% CrI Coverage	Mean CrI Width	Efficiency Gain (vs Direct)	Runtime (s/fit)
Direct Estimator (Sample $\hat{p}_d$ )	0.0681	26.27%	20.71%	+2.27%	—	—	1.000x	—
`fastsae::hb_area` (Spatial Besag)	0.04913	19.31%	16.09%	+2.36%	92.13%	0.1855	1.783x	6.18 s
`tipsae::fit_sae` (Spatial Besag)	0.04911	19.16%	15.89%	+1.06%	93.18%	0.1979	1.806x	11.65 s

Key Empirical Findings:

Substantial Variance & Error Reduction: Both fastsae and tipsae reduce the empirical RMSE from 0.0681 to 0.0491, an empirical efficiency gain of about 1.8x relative to direct survey estimation. The two Bayesian fits are practically equally accurate.
Uncertainty Quantification: fastsae gives slightly narrower credible intervals ( $0.186$ vs $0.198$ ) with a coverage of 92.1%, against 93.2% for tipsae; both are a little below the nominal 95%.
Computational Scalability: fastsae::hb_area fits in about half the time of the Stan MCMC run here; its Beta model with vardir is fitted twice (see ?hb_area).

An earlier run of this simulation used a PC prior with param = c(1, 0.05) on the area SD and the precision of the Beta likelihood at the direct estimate; its intervals covered 88.3% of the true proportions (RMSE 0.0485, mean width 0.172). The current defaults widen the intervals and bring the coverage close to that of tipsae.

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).

Estimator / Model	Empirical RMSE	Empirical RRMSE (%)	Empirical MARB (%)	Empirical RB (%)	95% CrI Coverage (%)	Mean CrI Width	Efficiency Gain (vs Direct)	Runtime (s/fit)
Direct Estimator (Sample $\hat{p}_d$ )	0.0719	27.14%	21.75%	+2.96%	—	—	1.000x	—
`fastsae::hb_area` (Spatial Besag)	0.0489	18.85%	15.82%	+2.81%	93.38%	0.1935	1.991x	8.53 s
`tipsae::fit_sae` (Spatial Besag)	0.0487	18.68%	15.66%	+2.00%	93.69%	0.2017	2.014x	40.00 s

Key National Scalability Takeaways:

4.7x Speedup: fastsae fits the full 500-domain spatial model in $8.5$ seconds, compared to $40.0$ seconds for Stan MCMC (one chain of 500 iterations).
Equivalent accuracy: the empirical RMSE of the two Bayesian fits is practically the same ( $0.0489$ vs $0.0487$ , against $0.0719$ for the direct estimator).
Coverage: 93.4% for fastsae and 93.7% for tipsae, with slightly narrower intervals for fastsae ( $0.194$ vs $0.202$ ). The earlier run with a PC prior and the precision at the direct estimate covered 89.0%.

Visualizations

1. Domain-Specific Empirical RMSE

Across all 60 small areas, both fastsae and tipsae consistently halve the estimation variance relative to direct survey estimates:

# 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")

2. Empirical 95% Credible Interval Coverage

The empirical coverage rates of the domains scatter around the nominal 95% dashed reference line, on average 92.1% for fastsae and 93.2% for tipsae.

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.