geoaddSAE2-intro

Background

geosae fits the area-level Geoadditive Small Area Estimation (Geoadditive SAE) model: a semiparametric extension of the Fay-Herriot model in which

linear covariate effects are modelled as usual,
nonlinear covariate effects are modelled with penalized splines (P-splines), and
spatial variation is modelled with a bivariate thin plate regression spline,

all represented jointly as a linear mixed model and fitted by Restricted Maximum Likelihood (REML). The Mean Squared Error (MSE) of the resulting small area predictor is obtained by parametric bootstrap.

The classical Fay-Herriot (FH) and Spatial Fay-Herriot (SFH) models are not re-implemented: geosae calls the existing implementations by the sae package so that the three approaches can be fitted and compared consistently.

A minimal example

library(geoaddSAE2)
data(simulated_sae)
head(simulated_sae)
#>      area          y x1       x2         x3        lat      lon     vardir
#> 1 Area_01 5.75116528  1 4.138301  3.0539692  -6.111182 122.5995 0.06666667
#> 2 Area_02 0.04273984  0 1.037720 -2.2803723   2.401187 110.3099 0.04000000
#> 3 Area_03 6.61571316  1 4.116264  2.5466352  -4.047392 117.4762 0.04166667
#> 4 Area_04 6.23637659  1 3.917563  0.4794186   4.011296 138.9058 0.03571429
#> 5 Area_05 3.90330728  1 3.520527 -0.6569142   4.987944 117.2135 0.08333333
#> 6 Area_06 4.18112746  1 2.923643 -0.3154003 -10.225540 135.9561 0.02380952
#>          theta
#> 1  5.899357370
#> 2 -0.001369921
#> 3  6.786057967
#> 4  6.261011210
#> 5  3.739533377
#> 6  3.994656977

simulated_sae has a direct estimator y, a known sampling variance vardir, a linear covariate x1, a nonlinear covariate x2, and spatial coordinates lat/lon.

Geoadditive SAE only

fit <- geosae(
  data      = simulated_sae,
  formula   = y ~ x1,
  vardir    = vardir,
  nonlinear = "x2",
  spatial   = c("lat", "lon"),
  B         = 50,
  seed      = 1
)

fit
#> === Geoadditive Small Area Estimation Model ===
#> 
#> Formula (linear) : y ~ x1 
#> Nonlinear terms  : x2 
#> Spatial terms    : lat, lon 
#> Number of Areas  : 100 
#> 
#>         area direct_est direct_mse geoadditive_est geoadditive_mse
#> 1    Area_01 5.75116528 0.06666667       6.5584718     0.021507138
#> 2    Area_02 0.04273984 0.04000000      -0.2194713     0.027056796
#> 3    Area_03 6.61571316 0.04166667       4.6159626     0.016272767
#> 4    Area_04 6.23637659 0.03571429       5.2761556     0.018154033
#> 5    Area_05 3.90330728 0.08333333       5.6670604     0.031721017
#> 6    Area_06 4.18112746 0.02380952       5.0320352     0.016113246
#> 7    Area_07 5.11069626 0.02380952       4.9811352     0.017078133
#> 8    Area_08 1.26974831 0.05000000       1.5226914     0.025793983
#> 9    Area_09 2.53503718 0.02222222       3.5927593     0.019069102
#> 10   Area_10 6.37934945 0.05263158       5.0739790     0.018862764
#> 11   Area_11 1.94237636 0.03030303       3.1567750     0.020523867
#> 12   Area_12 6.16401327 0.02702703       6.0470741     0.013692276
#> 13   Area_13 3.21091033 0.02083333       2.8855341     0.015969708
#> 14   Area_14 2.64910648 0.03703704       2.1242307     0.034773287
#> 15   Area_15 1.96870917 0.03846154       3.3526977     0.022162359
#> 16   Area_16 4.02261620 0.03448276       4.2250664     0.028777049
#> 17   Area_17 6.62251124 0.05555556       5.8135386     0.021529082
#> 18   Area_18 1.78206408 0.05000000       1.8398888     0.042404317
#> 19   Area_19 3.67762650 0.05882353       4.5615501     0.023133610
#> 20   Area_20 4.53683605 0.07142857       3.8141823     0.041533565
#> 21   Area_21 6.54020334 0.02631579       6.4190328     0.015254173
#> 22   Area_22 4.17529738 0.06666667       3.3563247     0.022375091
#> 23   Area_23 3.61164478 0.05263158       5.0499082     0.023353784
#> 24   Area_24 3.26737137 0.02857143       3.4429183     0.022572405
#> 25   Area_25 4.31371129 0.05555556       2.8398809     0.012530196
#> 26   Area_26 4.29320168 0.05000000       4.4880609     0.036623890
#> 27   Area_27 4.79593009 0.06666667       5.6118172     0.053183809
#> 28   Area_28 1.39586453 0.02439024       1.3447014     0.020906532
#> 29   Area_29 1.75747607 0.04761905       2.0704647     0.016888943
#> 30   Area_30 7.03054359 0.03571429       5.2990814     0.023559425
#> 31   Area_31 5.71163406 0.02777778       5.5768790     0.013523816
#> 32   Area_32 4.26725955 0.02439024       2.9293680     0.013310307
#> 33   Area_33 4.06595513 0.02777778       5.5599192     0.019003137
#> 34   Area_34 5.35081995 0.02564103       5.4551451     0.015787673
#> 35   Area_35 1.61017708 0.07692308       2.2189151     0.043365623
#> 36   Area_36 5.85728459 0.03030303       5.0806900     0.021185932
#> 37   Area_37 4.01020285 0.02702703       4.3947065     0.026027452
#> 38   Area_38 4.60056997 0.08333333       4.6976748     0.044418646
#> 39   Area_39 4.25093185 0.02040816       4.6713304     0.016613418
#> 40   Area_40 4.61888304 0.06666667       6.4351168     0.020398292
#> 41   Area_41 2.05759723 0.02222222       2.2493178     0.021959039
#> 42   Area_42 4.02118521 0.02564103       3.1208296     0.010592489
#> 43   Area_43 6.88024386 0.05263158       6.1245959     0.056845280
#> 44   Area_44 1.56221280 0.02941176       2.7819417     0.011989826
#> 45   Area_45 5.98058025 0.02127660       5.7648721     0.010667476
#> 46   Area_46 5.96758672 0.02222222       4.6617281     0.008909519
#> 47   Area_47 8.32962309 0.02000000       7.8744444     0.011747158
#> 48   Area_48 1.43097959 0.02941176       2.7227587     0.019769323
#> 49   Area_49 1.01745484 0.07692308       2.8736794     0.026301599
#> 50   Area_50 6.88856129 0.04166667       4.3946852     0.019989585
#> 51   Area_51 3.68527035 0.03225806       4.0760400     0.025988901
#> 52   Area_52 4.67601131 0.03030303       5.0455598     0.019095078
#> 53   Area_53 2.50373745 0.02083333       3.9137410     0.007286870
#> 54   Area_54 2.63136965 0.02173913       3.2616952     0.015090068
#> 55   Area_55 3.18783384 0.04166667       4.0419625     0.022261292
#> 56   Area_56 7.87355569 0.04000000       7.5428429     0.018283139
#> 57   Area_57 7.12741366 0.03571429       6.6575681     0.035425589
#> 58   Area_58 5.65875629 0.05882353       4.6237529     0.020822989
#> 59   Area_59 2.30968302 0.02500000       1.9851610     0.015979246
#> 60   Area_60 4.47627250 0.02173913       4.7240893     0.019235812
#> 61   Area_61 2.38392852 0.02083333       2.1114553     0.018373663
#> 62   Area_62 2.53208726 0.04000000       2.4854768     0.016750046
#> 63   Area_63 5.56855882 0.05555556       6.6214332     0.023519343
#> 64   Area_64 0.92735353 0.05000000       1.1653985     0.016368407
#> 65   Area_65 4.88855184 0.02564103       3.7606831     0.019966752
#> 66   Area_66 3.31214984 0.05263158       5.3823841     0.013930276
#> 67   Area_67 1.52726950 0.02325581       2.3133485     0.016404268
#> 68   Area_68 4.40878855 0.02127660       4.4734425     0.014926272
#> 69   Area_69 2.22273147 0.03448276       0.6820827     0.015245959
#> 70   Area_70 6.42408535 0.04347826       6.0177516     0.016022447
#> 71   Area_71 3.35656716 0.06666667       3.5862996     0.016581429
#> 72   Area_72 5.46313040 0.03703704       5.1316058     0.015774639
#> 73   Area_73 4.72021401 0.09090909       5.2895530     0.048032836
#> 74   Area_74 6.53798206 0.05882353       4.7040736     0.032786770
#> 75   Area_75 6.27414527 0.05000000       4.6479569     0.023467529
#> 76   Area_76 6.20887050 0.03333333       6.0260522     0.018904108
#> 77   Area_77 5.71789021 0.06250000       4.9534653     0.024554384
#> 78   Area_78 3.05226922 0.03125000       2.3827959     0.014814459
#> 79   Area_79 2.02207807 0.05882353       2.7679932     0.023619817
#> 80   Area_80 8.60147781 0.05263158       8.0247737     0.017236775
#> 81   Area_81 3.09721018 0.02173913       3.4336344     0.019804194
#> 82   Area_82 7.41660997 0.03846154       6.6645883     0.017538824
#> 83   Area_83 2.28788908 0.06666667       3.2020414     0.029333122
#> 84   Area_84 7.22179102 0.02040816       7.8235886     0.011028152
#> 85   Area_85 5.76386565 0.06666667       6.1391598     0.032499964
#> 86   Area_86 4.67584433 0.02702703       5.5000523     0.018748214
#> 87   Area_87 1.65961121 0.02439024       1.5172285     0.018917163
#> 88   Area_88 5.01709843 0.02272727       5.1823634     0.018209458
#> 89   Area_89 1.75334208 0.02702703       2.6106353     0.016110202
#> 90   Area_90 5.33088028 0.03448276       2.7160922     0.011258881
#> 91   Area_91 7.01441155 0.10000000       8.4245453     0.016594209
#> 92   Area_92 2.78415535 0.07142857       3.9432778     0.024295958
#> 93   Area_93 4.74526632 0.05000000       4.0209464     0.025945634
#> 94   Area_94 7.77533877 0.03448276       7.1847548     0.014407159
#> 95   Area_95 1.47097284 0.09090909       3.4589957     0.019977457
#> 96   Area_96 3.61976837 0.05263158       2.9538557     0.026825693
#> 97   Area_97 0.40577186 0.06250000       2.2767197     0.030972634
#> 98   Area_98 2.25785409 0.02127660       3.0580487     0.011413899
#> 99   Area_99 4.84838972 0.02222222       3.6453263     0.016277527
#> 100 Area_100 8.27389818 0.03125000       7.8241605     0.014522433
#> 
#> * Note: Use summary() for model diagnostics or extract $estimation for full results.

summary(fit)
#> === Summary: Geoadditive Small Area Estimation Model ===
#> 
#> Nonlinear Terms & Knots:
#>   - x2 : k = 25
#> 
#> Model Diagnostics:
#>   Convergence     : TRUE
#>   Deviance Expl.  : 78.43%
#>   R-squared (adj) : 0.5310
#> 
#> Fixed-Effect Parameters:
#>         term  estimate  std_error      p_value
#>  (Intercept) 3.7845809 0.08884725 0.000000e+00
#>           x1 0.5661607 0.09451434 2.095703e-09
#> 
#> Smooth Terms:
#>        term     edf   ref_df p_value
#>       s(x2) 23.8658 23.99319       0
#>  s(lat,lon) 28.5986 28.98766       0
#> 
#> MSE Summary across areas (Parametric Bootstrap):
#>               Value
#> Mean MSE 0.02157212
#> RMSE     0.14687450

Comparing Geoadditive SAE, Fay-Herriot, and Spatial Fay-Herriot

Set compare = TRUE to also fit FH (always) and SFH (whenever spatial is supplied), and get a side-by-side comparison table.

fit_cmp <- geosae(
  data      = simulated_sae,
  formula   = y ~ x1,
  vardir    = vardir,
  nonlinear = "x2",
  spatial   = c("lat", "lon"),
  compare   = TRUE,
  B         = 50,
  seed      = 1
)
fit_cmp$comparison
#>                 model        mse      rmse
#> 1     Direct Estimate 0.04193097 0.2047705
#> 2         Fay-Herriot 0.04108135 0.2026853
#> 3 Spatial Fay-Herriot 0.04110830 0.2027518
#> 4     Geoadditive SAE 0.02157212 0.1468745

Inspecting results

A fitted geosae object cleanly separates:

fit$estimation – direct estimate, geoadditive estimate, MSE, RMSE per area;
fit$diagnostics – convergence, EDF of smooth terms, variance components, R-sq, deviance explained;
fit$parameters – fixed-effect (linear) coefficient table;
fit$comparison – model | mse | rmse, when compare = TRUE;
fit$models – the underlying fitted model objects, for advanced use.