geoaddSAE2-intro

The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

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.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.