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"),
  bootstrap = TRUE,
  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.021150604
#> 2    Area_02 0.04273984 0.04000000      -0.2194713     0.033486455
#> 3    Area_03 6.61571316 0.04166667       4.6159626     0.013043208
#> 4    Area_04 6.23637659 0.03571429       5.2761556     0.024534637
#> 5    Area_05 3.90330728 0.08333333       5.6670604     0.037904679
#> 6    Area_06 4.18112746 0.02380952       5.0320352     0.016834623
#> 7    Area_07 5.11069626 0.02380952       4.9811352     0.014053099
#> 8    Area_08 1.26974831 0.05000000       1.5226914     0.026653680
#> 9    Area_09 2.53503718 0.02222222       3.5927593     0.016580653
#> 10   Area_10 6.37934945 0.05263158       5.0739790     0.019103575
#> 11   Area_11 1.94237636 0.03030303       3.1567750     0.015026321
#> 12   Area_12 6.16401327 0.02702703       6.0470741     0.014352974
#> 13   Area_13 3.21091033 0.02083333       2.8855341     0.010788376
#> 14   Area_14 2.64910648 0.03703704       2.1242307     0.030220000
#> 15   Area_15 1.96870917 0.03846154       3.3526977     0.018626148
#> 16   Area_16 4.02261620 0.03448276       4.2250664     0.029810333
#> 17   Area_17 6.62251124 0.05555556       5.8135386     0.025658514
#> 18   Area_18 1.78206408 0.05000000       1.8398888     0.031689509
#> 19   Area_19 3.67762650 0.05882353       4.5615501     0.026043168
#> 20   Area_20 4.53683605 0.07142857       3.8141823     0.032142336
#> 21   Area_21 6.54020334 0.02631579       6.4190328     0.015278202
#> 22   Area_22 4.17529738 0.06666667       3.3563247     0.023170666
#> 23   Area_23 3.61164478 0.05263158       5.0499082     0.018167386
#> 24   Area_24 3.26737137 0.02857143       3.4429183     0.013133334
#> 25   Area_25 4.31371129 0.05555556       2.8398809     0.016326903
#> 26   Area_26 4.29320168 0.05000000       4.4880609     0.027711808
#> 27   Area_27 4.79593009 0.06666667       5.6118172     0.040046917
#> 28   Area_28 1.39586453 0.02439024       1.3447014     0.020863228
#> 29   Area_29 1.75747607 0.04761905       2.0704647     0.013748396
#> 30   Area_30 7.03054359 0.03571429       5.2990814     0.023153900
#> 31   Area_31 5.71163406 0.02777778       5.5768790     0.022340726
#> 32   Area_32 4.26725955 0.02439024       2.9293680     0.015602263
#> 33   Area_33 4.06595513 0.02777778       5.5599192     0.014903131
#> 34   Area_34 5.35081995 0.02564103       5.4551451     0.017430229
#> 35   Area_35 1.61017708 0.07692308       2.2189151     0.034535550
#> 36   Area_36 5.85728459 0.03030303       5.0806900     0.018968088
#> 37   Area_37 4.01020285 0.02702703       4.3947065     0.025568021
#> 38   Area_38 4.60056997 0.08333333       4.6976748     0.045191711
#> 39   Area_39 4.25093185 0.02040816       4.6713304     0.012532009
#> 40   Area_40 4.61888304 0.06666667       6.4351168     0.018502707
#> 41   Area_41 2.05759723 0.02222222       2.2493178     0.015380445
#> 42   Area_42 4.02118521 0.02564103       3.1208296     0.015681542
#> 43   Area_43 6.88024386 0.05263158       6.1245959     0.072910354
#> 44   Area_44 1.56221280 0.02941176       2.7819417     0.014616352
#> 45   Area_45 5.98058025 0.02127660       5.7648721     0.008091791
#> 46   Area_46 5.96758672 0.02222222       4.6617281     0.015106086
#> 47   Area_47 8.32962309 0.02000000       7.8744444     0.013644239
#> 48   Area_48 1.43097959 0.02941176       2.7227587     0.014208174
#> 49   Area_49 1.01745484 0.07692308       2.8736794     0.034466187
#> 50   Area_50 6.88856129 0.04166667       4.3946852     0.024273463
#> 51   Area_51 3.68527035 0.03225806       4.0760400     0.021092538
#> 52   Area_52 4.67601131 0.03030303       5.0455598     0.015186498
#> 53   Area_53 2.50373745 0.02083333       3.9137410     0.012668616
#> 54   Area_54 2.63136965 0.02173913       3.2616952     0.013312870
#> 55   Area_55 3.18783384 0.04166667       4.0419625     0.023615521
#> 56   Area_56 7.87355569 0.04000000       7.5428429     0.011129045
#> 57   Area_57 7.12741366 0.03571429       6.6575681     0.019162342
#> 58   Area_58 5.65875629 0.05882353       4.6237529     0.025963958
#> 59   Area_59 2.30968302 0.02500000       1.9851610     0.015702140
#> 60   Area_60 4.47627250 0.02173913       4.7240893     0.014292655
#> 61   Area_61 2.38392852 0.02083333       2.1114553     0.017386113
#> 62   Area_62 2.53208726 0.04000000       2.4854768     0.015512306
#> 63   Area_63 5.56855882 0.05555556       6.6214332     0.019895426
#> 64   Area_64 0.92735353 0.05000000       1.1653985     0.020448694
#> 65   Area_65 4.88855184 0.02564103       3.7606831     0.015044426
#> 66   Area_66 3.31214984 0.05263158       5.3823841     0.010856887
#> 67   Area_67 1.52726950 0.02325581       2.3133485     0.016795905
#> 68   Area_68 4.40878855 0.02127660       4.4734425     0.016129507
#> 69   Area_69 2.22273147 0.03448276       0.6820827     0.014236756
#> 70   Area_70 6.42408535 0.04347826       6.0177516     0.012162971
#> 71   Area_71 3.35656716 0.06666667       3.5862996     0.024148744
#> 72   Area_72 5.46313040 0.03703704       5.1316058     0.023269390
#> 73   Area_73 4.72021401 0.09090909       5.2895530     0.030655253
#> 74   Area_74 6.53798206 0.05882353       4.7040736     0.026067722
#> 75   Area_75 6.27414527 0.05000000       4.6479569     0.029302346
#> 76   Area_76 6.20887050 0.03333333       6.0260522     0.029446345
#> 77   Area_77 5.71789021 0.06250000       4.9534653     0.018619757
#> 78   Area_78 3.05226922 0.03125000       2.3827959     0.014070677
#> 79   Area_79 2.02207807 0.05882353       2.7679932     0.027160222
#> 80   Area_80 8.60147781 0.05263158       8.0247737     0.019701001
#> 81   Area_81 3.09721018 0.02173913       3.4336344     0.023372044
#> 82   Area_82 7.41660997 0.03846154       6.6645883     0.013760574
#> 83   Area_83 2.28788908 0.06666667       3.2020414     0.023625796
#> 84   Area_84 7.22179102 0.02040816       7.8235886     0.012488577
#> 85   Area_85 5.76386565 0.06666667       6.1391598     0.032806274
#> 86   Area_86 4.67584433 0.02702703       5.5000523     0.020980709
#> 87   Area_87 1.65961121 0.02439024       1.5172285     0.014866415
#> 88   Area_88 5.01709843 0.02272727       5.1823634     0.025795404
#> 89   Area_89 1.75334208 0.02702703       2.6106353     0.020045075
#> 90   Area_90 5.33088028 0.03448276       2.7160922     0.011662084
#> 91   Area_91 7.01441155 0.10000000       8.4245453     0.020021971
#> 92   Area_92 2.78415535 0.07142857       3.9432778     0.019951082
#> 93   Area_93 4.74526632 0.05000000       4.0209464     0.026326195
#> 94   Area_94 7.77533877 0.03448276       7.1847548     0.021310251
#> 95   Area_95 1.47097284 0.09090909       3.4589957     0.013389205
#> 96   Area_96 3.61976837 0.05263158       2.9538557     0.030311746
#> 97   Area_97 0.40577186 0.06250000       2.2767197     0.025488004
#> 98   Area_98 2.25785409 0.02127660       3.0580487     0.016067858
#> 99   Area_99 4.84838972 0.02222222       3.6453263     0.011528485
#> 100 Area_100 8.27389818 0.03125000       7.8241605     0.016480845
#> 
#> * 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 (Bootstrap):
#>              Min       Mean        Max
#> MSE  0.008091791 0.02096572 0.07291035
#> RMSE 0.089954384 0.14214669 0.27001917

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         Fay-Herriot 0.04142231 0.1986204
#> 2 Spatial Fay-Herriot 0.04143614 0.1986506
#> 3     Geoadditive SAE 0.02096572 0.1421467

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.