Diagnostic Evaluation of Small Area Estimation Models

Conducts a comprehensive diagnostic evaluation of small area estimation models (HB, EBP, or EBLUP). Evaluates estimation precision, efficiency gains over direct estimators, external calibration and bias diagnostics (Brown et al., 2001), goodness-of-fit tests, and residual spatial autocorrelation (Moran's I).

Usage

diagnose(
  object,
  W = NULL,
  truth = NULL,
  rse_threshold = 25,
  alpha_level = 0.05
)

Arguments

object: A fitted model object of class "fastsae" (e.g., from hb_area, eblup_fh, eblup_sfh, or eblup_stfh).
W: Optional spatial proximity or adjacency matrix of dimension \(D \times D\). If NULL and the fitted model contains a spatial matrix (e.g. eblup_sfh or spatial hb_area), it is automatically extracted from object.
truth: Optional numeric vector containing true parameter values \(\theta_d\) (useful for simulation studies).
rse_threshold: Numeric. Threshold for reliable Relative Standard Error (default is 25%).
alpha_level: Numeric. Significance level for hypothesis tests (default is 0.05).

Value

An object of class c("fastsae_diagnose", "list") containing:

model_info: Basic model characteristics (model name, sample size, unsampled count).
precision: Summary of direct vs SAE RSE, proportion of areas with RSE < threshold, and MSE reduction ratios.
brown_test: Results of the Brown et al. (2001) calibration test (\(H_0: \alpha = 0, \beta = 1\)) and Chi-Square goodness-of-fit statistic \(W\).
spatial_test: Moran's I test for spatial autocorrelation in model residuals (if W is available).
simulation_metrics: Relative bias (RB), Relative RMSE, and coverage rates (if truth is supplied).
df_diag: Data frame combining direct estimates, SAE predictions, MSE, RSE, and residuals for diagnostics.
status: Overall assessment summary string.

References

Brown, G., Chambers, R., Heady, P., & Heasman, D. (2001). Evaluation of small area estimation methods: An application to the British Labour Force Survey. ONS Internal Report.
Rao, J. N. K., & Molina, I. (2015). Small Area Estimation (2nd ed.). John Wiley & Sons.
Cliff, A. D., & Ord, J. K. (1981). Spatial Processes: Models & Applications. Pion London.

Examples

library(fastsae)
data(mys)
data(mys_proxmat)

# 1. Fit Fay-Herriot model
fit_fh <- eblup_fh(y ~ x1 + x2, data = mys, vardir = ~vardir, domain = ~area)
#> 
#> ── Fast Small Area Estimation (fastsae) ────────────────────────────────────────
#> Call:
#> eblup_fh(formula = y ~ x1 + x2, vardir = ~vardir, domain = ~area, data = mys)
#> 
#> ✔ Convergence: Yes (in 7 iterations)
#> Model: Fay-Herriot (Area-level)
#> Method: eblup
#> Random effect variance (sigma2_u): 2.569182 
#> 
#> Fixed Effects Coefficients:
#>                   beta  std.error     zvalue pvalue
#> (Intercept)  3.0689476  0.7631917  4.0212018 0.0001
#> x1          -0.0041073  0.0090768 -0.4525074 0.6509
#> x2           0.0861173  0.0309576  2.7817847 0.0054
#> 
#> EBLUP Estimates (First 6 domains):
#>   domain        y    eblup    vardir random_effect       mse       rse
#> 1      1 8.359527 7.594807 0.6618838    2.96835204 0.5602338  9.855256
#> 2      2 7.599650 6.777056 0.8374691    2.52354833 0.6766501 12.137828
#> 3      3 5.514137 5.247512 0.8822257    0.77645304 0.7048881 15.999508
#> 4      4 3.869326 4.247072 0.6581716   -1.47453645 0.5556741 17.551750
#> 5      5 6.305063 6.353632 1.2788021   -0.09757891 0.9308379 15.185005
#> 6      6 3.926807 4.090117 0.3878004   -1.08192989 0.3497100 14.458337
#> ... and 36 more rows.
#> 
diag_fh <- diagnose(fit_fh)
print(diag_fh)
#> ── fastsae Small Area Estimation Diagnostic Report ─────────────────────────────
#> Model: "FH" | Domains: 42 (Sampled: 32, Unsampled: 10)
#> 
#> ── 1. Precision & Efficiency Gain ──
#> 
#> ! Domains with RSE < 25%: 52.4% (Caution: low precision)
#> ✔ Average RSE reduction: Direct "29.49%" -> SAE "25.59%" (Gain: 3.9%)
#> ✔ Variance reduction in 100% of areas (MSE ratio median: 1.33, max: 4.44)
#> 
#> ── 2. Brown et al. (2001) Calibration Tests ──
#> 
#> ! Bias Regression Test (H0: alpha = 0, beta = 1): F = 3.503, p-value = 0.0429 [Potential Systematic Bias]
#> Estimated parameters: alpha = -0.723, beta = 1.1587
#> ! Goodness-of-Fit Statistic W (Chi-Square): W = 4.57 (df = 32), p-value = 1 [Deviation from Survey Variance]
#> ────────────────────────────────────────────────────────────────────────────────
#> ! Final Assessment: CAUTION: Potential systematic bias detected between direct and model predictions.
#> ────────────────────────────────────────────────────────────────────────────────

# 2. Fit Spatial HB model
fit_hb <- hb_area(y ~ x1 + x2, data = mys, vardir = "vardir", spatial = "bym2", W = mys_proxmat)
#> 
#> ── Fast Small Area Estimation (fastsae) ────────────────────────────────────────
#> Call:
#> hb_area(formula = y ~ x1 + x2, data = mys, spatial = "bym2", W = mys_proxmat,
#> vardir = "vardir")
#> 
#> ✔ Convergence: Yes (in - iterations)
#> Model: HB-GAUSSIAN (BYM2)
#> Random effect variance (sigma2_u): 1.678076 
#> Spatial autocorrelation (rho): 0.4817 
#> Spatial mixing fraction (phi): 0.4817 
#> 
#> Fixed Effects Coefficients:
#>                    beta   std.error      zvalue      pvalue    ci_lower
#> (Intercept)  2.9993e+00  6.7355e-01  4.4529e+00  8.4699e-06  1.6908e+00
#> x1          -3.6465e-03  8.0773e-03 -4.5145e-01  6.5167e-01 -1.9694e-02
#> x2           8.6081e-02  2.7821e-02  3.0942e+00  1.9738e-03  3.1295e-02
#>             ci_upper
#> (Intercept)   4.3416
#> x1            0.0121
#> x2            0.1409
#> 
#> HB Estimates (First 6 domains):
#>   domain        y       hb linear_pred        sd       mse      rse ci_lower
#> 1      1 8.359527 7.336937    7.336937 0.7396177 0.5470343 10.08074 5.902625
#> 2      2 7.599650 6.515022    6.515022 0.7988455 0.6381542 12.26159 4.970633
#> 3      3 5.514137 5.150737    5.150737 0.7839678 0.6146055 15.22050 3.623152
#> 4      4 3.869326 4.362793    4.362793 0.7078676 0.5010765 16.22510 2.964204
#> 5      5 6.305063 6.362973    6.362973 0.8812813 0.7766567 13.85015 4.631067
#> 6      6 3.926807 4.146942    4.146942 0.5685674 0.3232689 13.71052 3.028191
#>   ci_upper random_effect    vardir
#> 1 8.802104    2.72064313 0.6618838
#> 2 8.102565    2.28783016 0.8374691
#> 3 6.699835    0.72250797 0.8822257
#> 4 5.741249   -1.32713961 0.6581716
#> 5 8.092337   -0.08182539 1.2788021
#> 6 5.258483   -0.99782653 0.3878004
#> ... and 36 more rows.
#> 
diag_hb <- diagnose(fit_hb)
print(diag_hb)
#> ── fastsae Small Area Estimation Diagnostic Report ─────────────────────────────
#> Model: "HB-GAUSSIAN (BYM2)" | Domains: 42 (Sampled: 32, Unsampled: 10)
#> 
#> ── 1. Precision & Efficiency Gain ──
#> 
#> ! Domains with RSE < 25%: 61.9% (Caution: low precision)
#> ✔ Average RSE reduction: Direct "29.49%" -> SAE "22.92%" (Gain: 6.57%)
#> ✔ Variance reduction in 100% of areas (MSE ratio median: 1.55, max: 5.98)
#> 
#> ── 2. Brown et al. (2001) Calibration Tests ──
#> 
#> ✔ Bias Regression Test (H0: alpha = 0, beta = 1): F = 2.995, p-value = 0.0652 [Statistically Unbiased]
#> Estimated parameters: alpha = -0.8042, beta = 1.1793
#> ! Goodness-of-Fit Statistic W (Chi-Square): W = 7.34 (df = 32), p-value = 1 [Deviation from Survey Variance]
#> 
#> ── 3. Residual Spatial Autocorrelation ──
#> 
#> ✔ Moran's I on Residuals: I = -0.0382 (Expected: -0.0323, p-value = 0.438) [No Residual Spatial Autocorrelation]
#> 
#> ── 5. Bayesian Information Criteria & Predictive Diagnostics ──
#> 
#> ℹ WAIC: 117.4 (p_eff: 13.98) | DIC: 118.44 (p_eff: 19.95)
#> ✔ PIT Calibration Test vs Uniform(0,1): D = 0.084, p-value = 0.9646 [Well-calibrated predictive distribution]
#> ✔ Leave-One-Out CPO: No numerical approximation issues (min CPO = 0.0087)
#> ────────────────────────────────────────────────────────────────────────────────
#> ! Final Assessment: CAUTION: Model goodness-of-fit indicates notable deviation from survey variance.
#> ────────────────────────────────────────────────────────────────────────────────