Introduction to geoaddSAE2

Background

The geoaddSAE2 package provides the geosae() function for fitting an area-level Geoadditive Small Area Estimation (Geoadditive SAE) model.

The Geoadditive SAE model extends the Fay-Herriot model by allowing the relationship between the direct estimator and auxiliary variables to contain linear, nonlinear, and spatial components. Specifically:

These components are jointly represented within a linear mixed model framework. Model parameters and smoothing parameters are estimated using Restricted Maximum Likelihood (REML).

The uncertainty of the resulting small area estimates is quantified using the Mean Squared Error (MSE) estimated by parametric bootstrap.

The classical Fay-Herriot (FH) and Spatial Fay-Herriot (SFH) models are not reimplemented in geoaddSAE2. When model comparison is requested, geosae() uses the existing implementations provided by the sae package. This allows Geoadditive SAE, FH, and SFH estimates to be compared within a consistent workflow.

Example data

The package includes the simulated_sae dataset, which can be used to demonstrate the main features of geosae().

library(geoaddSAE2)

data(simulated_sae)

head(simulated_sae)
#>       area        y x1       x2         x3        lat      lon     vardir
#> 1 Area_001 4.826302  0 4.138301  3.0539692  -6.111182 122.5995 0.06666667
#> 2 Area_002 1.117877  1 1.037720 -2.2803723   2.401187 110.3099 0.04000000
#> 3 Area_003 6.690850  1 4.116264  2.5466352  -4.047392 117.4762 0.04166667
#> 4 Area_004 6.311514  1 3.917563  0.4794186   4.011296 138.9058 0.03571429
#> 5 Area_005 2.978444  0 3.520527 -0.6569142   4.987944 117.2135 0.08333333
#> 6 Area_006 4.256265  1 2.923643 -0.3154003 -10.225540 135.9561 0.02380952
#>      theta
#> 1 4.974495
#> 2 1.073767
#> 3 6.861195
#> 4 6.336148
#> 5 2.814671
#> 6 4.069794

The dataset contains:

Exploratory identification before modelling

The geosae() function requires the user to explicitly specify the variables included in the nonlinear component through the nonlinear argument and the coordinates used for the spatial component through the spatial argument.

Before fitting the model, two exploratory checks can be useful:

  1. examining whether continuous auxiliary variables show evidence of nonlinear relationships with the direct estimate; and
  2. examining whether residual spatial autocorrelation remains after accounting for auxiliary variables.

These procedures should be treated as exploratory diagnostics rather than strict decision rules. The final model specification should also consider the substantive context, theoretical justification, and characteristics of the data.

Identifying nonlinear relationships

A continuous auxiliary variable can first be examined using a weighted generalized additive model.

Because direct estimates may have different sampling variances, the model can be fitted using inverse sampling variance weights:

\[w_i = \frac{1}{D_i},\]

where \(D_i\) denotes the sampling variance of the direct estimate for area \(i\).

The following example compares a linear model with a penalized spline model for each candidate continuous covariate.

library(mgcv)
d <- simulated_sae
w <- 1 / d$vardir

nonlinear_candidates <- c("x2", "x3")

smooth_models <- lapply(
  nonlinear_candidates,
  function(v) {
    gam(
      as.formula(sprintf("y ~ s(%s, bs = 'ps')", v)),
      data = d,
      weights = w,
      method = "REML"
    )
  }
)

names(smooth_models) <- nonlinear_candidates

op <- par(mfrow = c(1, length(smooth_models)))

for (v in names(smooth_models)) {
  plot(
    smooth_models[[v]],
    residuals = TRUE,
    pch = 16,
    cex = 0.6,
    shade = TRUE,
    xlab = v,
    ylab = sprintf("s(%s)", v)
  )
}


par(op)

If a variable shows a clear nonlinear pattern, it can be specified in geosae() through the nonlinear argument.

For example:

nonlinear = "x3"

Multiple nonlinear covariates can be supplied as a character vector when supported by the model specification.

Identifying spatial autocorrelation

Spatial autocorrelation can be examined before fitting the Geoadditive SAE model by testing whether residuals from a nonspatial model still exhibit spatial structure.

One common diagnostic is Moran’s I.

In this example, a spatial neighborhood structure is constructed using the four nearest neighbors based on the area coordinates.

For polygon-based regional data, a contiguity-based neighborhood structure, such as queen contiguity, may be more appropriate.

coords <- cbind(d$lon, d$lat)

nb <- spdep::knn2nb(spdep::knearneigh(coords, k = 4))
lw <- spdep::nb2listw(nb, style = "W")

fit_nonspatial <- mgcv::gam(y ~ x1 + x2 + s(x3, bs = "ps"),
                            data = d,weights = 1 / vardir,method = "REML")

res_std <- residuals(fit_nonspatial,type = "pearson")
spdep::moran.test(res_std,lw)
#> 
#>  Moran I test under randomisation
#> 
#> data:  res_std  
#> weights: lw    
#> 
#> Moran I statistic standard deviate = 2.1994, p-value = 0.01393
#> alternative hypothesis: greater
#> sample estimates:
#> Moran I statistic       Expectation          Variance 
#>       0.132778232      -0.010101010       0.004220286

set.seed(1)
spdep::moran.mc(res_std,lw,nsim = 999)
#> 
#>  Monte-Carlo simulation of Moran I
#> 
#> data:  res_std 
#> weights: lw  
#> number of simulations + 1: 1000 
#> 
#> statistic = 0.13278, observed rank = 977, p-value = 0.023
#> alternative hypothesis: greater

The test is performed on Pearson residuals from the weighted linear model. This adjustment is useful because the direct estimates may have different sampling variances across areas.

Evidence of positive and statistically significant residual spatial autocorrelation indicates that observations from nearby areas tend to retain similar residual values after accounting for the auxiliary variables. In this situation, incorporating a spatial component may improve the model specification.

For regional polygon data, spatial neighbors can instead be defined using contiguity. For example:

nb <- spdep::poly2nb(regions, queen = TRUE)

lw <- spdep::nb2listw(
  nb,
  style = "W",
  zero.policy = TRUE
)

Isolated regions with no contiguous neighbors require additional consideration. Possible alternatives include nearest-neighbor or distance-based neighborhood structures.

The choice of spatial weights should reflect the geographical structure and substantive interpretation of spatial dependence in the application.

Fitting the Geoadditive SAE model

After identifying the model components, the Geoadditive SAE model can be fitted using geosae().

The following example specifies:

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

fit
#> === Geoadditive Small Area Estimation Model ===
#> 
#> Formula (linear) : y ~ x1 + x2 
#> Nonlinear terms  : x3 
#> Spatial terms    : lat, lon 
#> Number of Areas  : 100 
#> 
#>         area direct_est direct_mse geoadditive_est geoadditive_mse
#> 1   Area_001  4.8263024 0.06666667       4.8317411     0.008541429
#> 2   Area_002  1.1178770 0.04000000       1.0512350     0.003407798
#> 3   Area_003  6.6908503 0.04166667       6.9546945     0.004859286
#> 4   Area_004  6.3115137 0.03571429       6.3762786     0.007922885
#> 5   Area_005  2.9784444 0.08333333       2.7842691     0.009481147
#> 6   Area_006  4.2562646 0.02380952       4.0742504     0.005181869
#> 7   Area_007  4.1858334 0.02380952       4.1708974     0.003790017
#> 8   Area_008  0.3448855 0.05000000       0.5360236     0.004761747
#> 9   Area_009  1.6101743 0.02222222       1.6590522     0.004976931
#> 10  Area_010  5.4544866 0.05263158       5.2993304     0.003366953
#> 11  Area_011  1.0175135 0.03030303       1.2381039     0.005604833
#> 12  Area_012  5.2391504 0.02702703       5.1716319     0.003351813
#> 13  Area_013  2.2860475 0.02083333       2.3958023     0.006397475
#> 14  Area_014  2.7242436 0.03703704       2.8272652     0.008640358
#> 15  Area_015  1.0438463 0.03846154       1.1126536     0.005062714
#> 16  Area_016  4.0977534 0.03448276       3.8400553     0.003930846
#> 17  Area_017  5.6976484 0.05555556       5.6271837     0.006380121
#> 18  Area_018  0.8572012 0.05000000       0.8759385     0.006657705
#> 19  Area_019  3.7527637 0.05882353       3.7908691     0.006921371
#> 20  Area_020  5.6119732 0.07142857       5.4704984     0.013151335
#> 21  Area_021  5.6153405 0.02631579       5.4709595     0.009325997
#> 22  Area_022  5.2504345 0.06666667       5.2521847     0.004677019
#> 23  Area_023  3.6867819 0.05263158       4.2091554     0.006457031
#> 24  Area_024  3.3425085 0.02857143       3.4573936     0.005780851
#> 25  Area_025  3.3888484 0.05555556       3.2671164     0.006901163
#> 26  Area_026  3.3683388 0.05000000       3.0902378     0.006476772
#> 27  Area_027  3.8710672 0.06666667       3.5931605     0.007104666
#> 28  Area_028  1.4710017 0.02439024       1.3932164     0.009367307
#> 29  Area_029  1.8326132 0.04761905       2.2890876     0.004034720
#> 30  Area_030  8.1056807 0.03571429       7.9615521     0.005951294
#> 31  Area_031  5.7867712 0.02777778       5.5979052     0.006271434
#> 32  Area_032  3.3423967 0.02439024       3.4083279     0.005895901
#> 33  Area_033  4.1410923 0.02777778       3.7808755     0.004360223
#> 34  Area_034  4.4259571 0.02564103       4.7379904     0.010765111
#> 35  Area_035  0.6853142 0.07692308       0.9649417     0.005110834
#> 36  Area_036  4.9324217 0.03030303       4.8065263     0.005038765
#> 37  Area_037  3.0853400 0.02702703       3.1073020     0.004387124
#> 38  Area_038  4.6757071 0.08333333       5.1036175     0.005919271
#> 39  Area_039  3.3260690 0.02040816       3.1053762     0.006460869
#> 40  Area_040  4.6940202 0.06666667       4.9062152     0.004284432
#> 41  Area_041  2.1327344 0.02222222       2.3393480     0.004831009
#> 42  Area_042  4.0963224 0.02564103       4.0392262     0.007226144
#> 43  Area_043  5.9553810 0.05263158       5.6005354     0.006298414
#> 44  Area_044  1.6373499 0.02941176       1.7904855     0.003337508
#> 45  Area_045  5.0557174 0.02127660       5.1931964     0.005597241
#> 46  Area_046  6.0427239 0.02222222       6.0677074     0.003346478
#> 47  Area_047  7.4047602 0.02000000       7.5009800     0.005889805
#> 48  Area_048  2.5061167 0.02941176       2.5309242     0.004926716
#> 49  Area_049  2.0925920 0.07692308       1.9999334     0.007052792
#> 50  Area_050  6.9636984 0.04166667       6.8113523     0.004749930
#> 51  Area_051  2.7604075 0.03225806       2.5407523     0.011678879
#> 52  Area_052  3.7511485 0.03030303       3.8450014     0.006818838
#> 53  Area_053  2.5788746 0.02083333       2.5988258     0.003890244
#> 54  Area_054  1.7065068 0.02173913       1.5278395     0.005592667
#> 55  Area_055  3.2629710 0.04166667       3.0504882     0.004045454
#> 56  Area_056  7.9486928 0.04000000       7.9649409     0.005112815
#> 57  Area_057  6.2025508 0.03571429       6.1919265     0.004537242
#> 58  Area_058  4.7338934 0.05882353       5.1235831     0.004999437
#> 59  Area_059  1.3848202 0.02500000       1.1255494     0.004157284
#> 60  Area_060  4.5514097 0.02173913       4.5489963     0.004293020
#> 61  Area_061  3.4590657 0.02083333       3.4011259     0.008894127
#> 62  Area_062  2.6072244 0.04000000       2.6755887     0.007496457
#> 63  Area_063  5.6436960 0.05555556       5.7914929     0.006351591
#> 64  Area_064  2.0024907 0.05000000       2.1751083     0.003031910
#> 65  Area_065  4.9636890 0.02564103       4.9262858     0.006195432
#> 66  Area_066  3.3872870 0.05263158       3.7607276     0.003271910
#> 67  Area_067  0.6024067 0.02325581       0.3387339     0.004146864
#> 68  Area_068  3.4839257 0.02127660       3.3806968     0.003564383
#> 69  Area_069  1.2978686 0.03448276       1.1258902     0.004964376
#> 70  Area_070  6.4992225 0.04347826       6.3765617     0.003958574
#> 71  Area_071  3.4317043 0.06666667       3.1580549     0.004491638
#> 72  Area_072  4.5382675 0.03703704       4.6689256     0.004019283
#> 73  Area_073  3.7953512 0.09090909       4.1705211     0.008331701
#> 74  Area_074  5.6131192 0.05882353       5.5170237     0.006438238
#> 75  Area_075  6.3492824 0.05000000       6.4553298     0.004007624
#> 76  Area_076  6.2840076 0.03333333       5.8701900     0.007462594
#> 77  Area_077  6.7930274 0.06250000       7.0991415     0.006541877
#> 78  Area_078  2.1274064 0.03125000       2.1818279     0.028871258
#> 79  Area_079  1.0972152 0.05882353       0.9636694     0.005657068
#> 80  Area_080  8.6766150 0.05263158       8.5235853     0.002640405
#> 81  Area_081  3.1723473 0.02173913       3.1396031     0.005646318
#> 82  Area_082  6.4917471 0.03846154       6.3879896     0.004866521
#> 83  Area_083  1.3630262 0.06666667       1.2141609     0.004803740
#> 84  Area_084  6.2969282 0.02040816       6.3883452     0.003244818
#> 85  Area_085  4.8390028 0.06666667       4.7187231     0.006756251
#> 86  Area_086  3.7509815 0.02702703       3.9294310     0.004207409
#> 87  Area_087  0.7347484 0.02439024       0.7374531     0.003043322
#> 88  Area_088  4.0922356 0.02272727       4.3999587     0.003449429
#> 89  Area_089  0.8284792 0.02702703       1.0587830     0.005876176
#> 90  Area_090  5.4060174 0.03448276       5.3337220     0.006893037
#> 91  Area_091  6.0895487 0.10000000       5.8489799     0.008924722
#> 92  Area_092  1.8592925 0.07142857       1.8258583     0.004372337
#> 93  Area_093  3.8204035 0.05000000       4.0738184     0.006645939
#> 94  Area_094  7.8504759 0.03448276       7.6439270     0.002529567
#> 95  Area_095  1.5461100 0.09090909       2.1062171     0.008802101
#> 96  Area_096  3.6949055 0.05263158       3.5670914     0.015028637
#> 97  Area_097  1.4809090 0.06250000       1.6064617     0.004123376
#> 98  Area_098  1.3329912 0.02127660       1.4314250     0.005992331
#> 99  Area_099  3.9235269 0.02222222       3.8280151     0.003462983
#> 100 Area_100  8.3490353 0.03125000       8.5106176     0.005628432
#> 
#> * Note: Use summary() for model diagnostics or extract $estimation for full results.

A more detailed model summary can be obtained using:

summary(fit)
#> === Summary: Geoadditive Small Area Estimation Model ===
#> 
#> Nonlinear Terms & Knots:
#>   - x3 : k = 25
#> 
#> Model Diagnostics:
#>   Convergence     : TRUE
#>   Deviance Expl.  : 99.17%
#>   R-squared (adj) : 0.9903
#> 
#> Fixed-Effect Parameters:
#>         term estimate  std_error       p_value
#>  (Intercept) 0.446386 0.05463733  3.084390e-16
#>           x1 1.042165 0.04079298 5.830432e-144
#>           x2 1.002540 0.01692491  0.000000e+00
#> 
#> Smooth Terms:
#>        term       edf    ref_df     p_value
#>       s(x3) 10.245274 12.285616 0.000000000
#>  s(lat,lon)  2.000227  2.000448 0.000457049
#> 
#> MSE Summary across areas (Parametric Bootstrap):
#>                Value
#> Mean MSE 0.005999761
#> RMSE     0.077458124

The B argument specifies the number of parametric bootstrap replications used to estimate the MSE.

A small value such as B = 50 is useful for demonstration because it reduces computation time. A substantially larger value should generally be used in empirical applications when stable MSE estimates are required.

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

The compare argument can be used to fit benchmark area-level SAE models together with the Geoadditive SAE model.

Set:

compare = TRUE

to fit the Fay-Herriot model and, when spatial information is supplied, the Spatial Fay-Herriot model.

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

fit_cmp$comparison
#>                 model         mse       rmse
#> 1     Direct Estimate 0.041930967 0.20477052
#> 2         Fay-Herriot 0.039665201 0.19916124
#> 3 Spatial Fay-Herriot 0.039717104 0.19929150
#> 4     Geoadditive SAE 0.005999761 0.07745812

The resulting comparison table provides a concise summary of the estimates of model uncertainty across the fitted SAE approaches.

The purpose of this comparison is not to assume that one model is universally preferable. Instead, it allows users to evaluate whether incorporating nonlinear and spatial components provides meaningful improvements for the data under study.

Inspecting fitted results

A fitted geosae object stores the main results in separate components.

The area-level estimates and uncertainty measures are available from:

fit$estimation

This component contains quantities such as:

Model diagnostics are available from:

fit$diagnostics

Depending on the fitted specification, these diagnostics include information such as:

Estimated coefficients for the linear component are available from:

fit$parameters

When compare = TRUE, model comparison results are available from:

fit$comparison

The underlying fitted model objects are stored in:

fit$models

These objects can be used for more advanced diagnostics, visualization, or further model inspection.

Complete example

A complete Geoadditive SAE workflow can therefore be written as:

library(geoaddSAE2)

data(simulated_sae)

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

summary(fit)
#> === Summary: Geoadditive Small Area Estimation Model ===
#> 
#> Nonlinear Terms & Knots:
#>   - x3 : k = 25
#> 
#> Model Diagnostics:
#>   Convergence     : TRUE
#>   Deviance Expl.  : 99.17%
#>   R-squared (adj) : 0.9903
#> 
#> Fixed-Effect Parameters:
#>         term estimate  std_error       p_value
#>  (Intercept) 0.446386 0.05463733  3.084390e-16
#>           x1 1.042165 0.04079298 5.830432e-144
#>           x2 1.002540 0.01692491  0.000000e+00
#> 
#> Smooth Terms:
#>        term       edf    ref_df     p_value
#>       s(x3) 10.245274 12.285616 0.000000000
#>  s(lat,lon)  2.000227  2.000448 0.000457049
#> 
#> MSE Summary across areas (Parametric Bootstrap):
#>                Value
#> Mean MSE 0.005999761
#> RMSE     0.077458124
#> 
#> Model Comparison:
#>                model         mse       rmse
#>      Direct Estimate 0.041930967 0.20477052
#>          Fay-Herriot 0.039665201 0.19916124
#>  Spatial Fay-Herriot 0.039717104 0.19929150
#>      Geoadditive SAE 0.005999761 0.07745812

head(fit$estimation)
#>       area direct_est direct_mse geoadditive_est geoadditive_mse   fh_est
#> 1 Area_001   4.826302 0.06666667        4.831741     0.008541429 4.965387
#> 2 Area_002   1.117877 0.04000000        1.051235     0.003407798 1.120568
#> 3 Area_003   6.690850 0.04166667        6.954695     0.004859286 6.716339
#> 4 Area_004   6.311514 0.03571429        6.376279     0.007922885 6.283086
#> 5 Area_005   2.978444 0.08333333        2.784269     0.009481147 3.036155
#> 6 Area_006   4.256265 0.02380952        4.074250     0.005181869 4.254804
#>       fh_mse  sfh_est    sfh_mse
#> 1 0.06204274 4.973930 0.06218202
#> 2 0.03829867 1.113805 0.03835078
#> 3 0.03980050 6.725553 0.03984349
#> 4 0.03430318 6.286389 0.03432752
#> 5 0.07600067 3.032792 0.07604981
#> 6 0.02317052 4.255011 0.02318165

fit$parameters
#>          term estimate  std_error       p_value
#> 1 (Intercept) 0.446386 0.05463733  3.084390e-16
#> 2          x1 1.042165 0.04079298 5.830432e-144
#> 3          x2 1.002540 0.01692491  0.000000e+00

fit$diagnostics
#> $method
#> [1] "REML"
#> 
#> $converged
#> [1] TRUE
#> 
#> $n_obs
#> [1] 100
#> 
#> $n_smooth_terms
#> [1] 2
#> 
#> $edf_table
#>         term       edf    ref_df     p_value
#> 1      s(x3) 10.245274 12.285616 0.000000000
#> 2 s(lat,lon)  2.000227  2.000448 0.000457049
#> 
#> $variance_components
#>         term      std_dev
#> 1      s(x3) 0.2092375818
#> 2 s(lat,lon) 0.0001364843
#> 
#> $r_sq_adj
#> [1] 0.9903406
#> 
#> $deviance_explained
#> [1] 0.9917304
#> 
#> $reml_score
#>    REML 
#> 12.4034 
#> attr(,"Dp")
#> [1] 52.46465
#> attr(,"ml.pen")
#> [1] 36.21415
#> 
#> $bootstrap
#> $bootstrap$B
#> [1] 50
#> 
#> $bootstrap$n_failed
#> [1] 0

fit$comparison
#>                 model         mse       rmse
#> 1     Direct Estimate 0.041930967 0.20477052
#> 2         Fay-Herriot 0.039665201 0.19916124
#> 3 Spatial Fay-Herriot 0.039717104 0.19929150
#> 4     Geoadditive SAE 0.005999761 0.07745812

This workflow combines linear, nonlinear, and spatial modelling within an area-level SAE framework while also providing uncertainty estimation through parametric bootstrap and optional comparison with conventional Fay-Herriot models.