Small Area Estimation Using a Projection Estimator with a Multilevel Regression Model

Introduction

The saeproj.multilevel package provides tools for small area estimation using a projection estimator with a linear multilevel regression working model.

The method is designed for a two-survey setting:

a smaller model survey (data_model) containing the response variable and auxiliary predictors;
a larger projection survey (data_proj) containing auxiliary predictors and survey design information, but not the response variable.

The main function is:

sae_ml_linear()

The function fits a multilevel regression model using the model survey, predicts the response variable for units in the projection survey, aggregates predicted values by domain, and applies a design-based residual correction.

This vignette demonstrates the complete workflow using the simulated datasets included in the package.

Method Overview

Let $i$ denote a unit and $d$ denote a domain. A random-intercept multilevel regression working model can be written as:

\[ Y_{id} = \beta_0 + \beta_1 X_{1id} + \beta_2 X_{2id} + \cdots + \beta_p X_{pid} + u_d + e_{id}, \]

where:

$Y_{id}$ is the response variable for unit $i$ in domain $d$;
$X_{1id}, \ldots, X_{pid}$ are auxiliary variables;
$\beta_0, \beta_1, \ldots, \beta_p$ are regression coefficients;
$u_d$ is the domain-specific random effect;
$e_{id}$ is the unit-level residual.

The multilevel regression model is fitted using the smaller model survey. Predictions are then generated for all units in the larger projection survey.

For each domain $d$, the synthetic projection estimate is denoted by:

\[ \hat{Y}^{SYN}_d. \]

A design-based residual correction is calculated from the model survey:

\[ \hat{B}_d. \]

The final projection estimator is:

\[ \hat{Y}^{PR}_d = \hat{Y}^{SYN}_d + \hat{B}_d. \]

The plug-in variance estimator is:

\[ \widehat{\mathrm{Var}} \left( \hat{Y}^{PR}_d \right) = \widehat{\mathrm{Var}} \left( \hat{Y}^{SYN}_d \right) + \widehat{\mathrm{Var}} \left( \hat{B}_d \right). \]

The reported plug-in variance is approximate and does not fully account for uncertainty in the estimated multilevel regression parameters.

Load the Package and Example Data

library(saeproj.multilevel)

data("saeml_modelsvy")
data("saeml_projsvy")

The package includes two simulated survey datasets:

saeml_modelsvy is the smaller model-survey dataset and contains the target variable Y;
saeml_projsvy is the larger projection-survey dataset and does not contain the target variable Y.

dim(saeml_modelsvy)
#> [1] 250  11
dim(saeml_projsvy)
#> [1] 15000    10

head(saeml_modelsvy)
#> # A tibble: 6 × 11
#>   prov  kab_kota id_individu     Z1     Z2     X1    X2     X3     X4     Y
#>   <chr> <chr>          <int>  <dbl>  <dbl>  <dbl> <int>  <dbl>  <dbl> <dbl>
#> 1 35    1                 14 -0.282 -0.117 -0.495     0  1.10   1.94   40.4
#> 2 35    1                 98 -0.282 -0.117  0.925     0 -1.24   1.40  146. 
#> 3 35    1                134 -0.282 -0.117 -0.734     0 -0.402  1.13   68.9
#> 4 35    1                256 -0.282 -0.117 -0.122     0  1.39   0.276  29.1
#> 5 35    1                620 -0.282 -0.117  0.719     1 -1.55  -1.76   84.6
#> 6 35    2               1019 -1.31  -0.645 -0.175     1 -0.434 -1.77  105. 
#> # ℹ 1 more variable: WEIND <dbl>
head(saeml_projsvy)
#> # A tibble: 6 × 10
#>   prov  kab_kota id_individu     Z1     Z2     X1    X2     X3     X4 WEIND
#>   <chr> <chr>          <int>  <dbl>  <dbl>  <dbl> <int>  <dbl>  <dbl> <dbl>
#> 1 35    1                  1 -0.282 -0.117  0.995     1 -0.215  1.18   3.33
#> 2 35    1                  4 -0.282 -0.117  0.761     0 -0.355 -0.777  3.33
#> 3 35    1                 18 -0.282 -0.117  0.872     1 -0.910  0.405  3.33
#> 4 35    1                 21 -0.282 -0.117 -0.615     0 -0.656 -0.361  3.33
#> 5 35    1                 24 -0.282 -0.117  0.200     0 -0.408  0.203  3.33
#> 6 35    1                 28 -0.282 -0.117  0.423     1 -1.08   1.43   3.33

The model survey contains 250 observations from 50 domains, with five sampled units in each domain.

model_domain_count <- table(saeml_modelsvy$kab_kota)

c(
  n_domains = length(model_domain_count),
  min_units_per_domain = min(model_domain_count),
  max_units_per_domain = max(model_domain_count)
)
#>            n_domains min_units_per_domain max_units_per_domain 
#>                   50                    5                    5

The projection survey contains 15,000 observations from the same 50 domains, with 300 sampled units in each domain.

proj_domain_count <- table(saeml_projsvy$kab_kota)

c(
  n_domains = length(proj_domain_count),
  min_units_per_domain = min(proj_domain_count),
  max_units_per_domain = max(proj_domain_count)
)
#>            n_domains min_units_per_domain max_units_per_domain 
#>                   50                  300                  300

The response variable Y is available only in the model survey.

"Y" %in% names(saeml_modelsvy)
#> [1] TRUE

"Y" %in% names(saeml_projsvy)
#> [1] FALSE

Fit the Multilevel Projection Estimator

The following specification fits a random-intercept multilevel regression model with kab_kota as both the random-effect grouping variable and the domain variable.

result <- sae_ml_linear(
  formula = Y ~ X1 + X2 + X3 + X4 + Z1 + Z2 + (1 | kab_kota),
  data_model = saeml_modelsvy,
  data_proj = saeml_projsvy,
  domain = "kab_kota",
  cluster_ids = ~1,
  weight = "WEIND",
  strata = "kab_kota",
  summary_function = "mean"
)

In this example:

Y is the target variable;
X1, X2, X3, and X4 are unit-level auxiliary variables;
Z1 and Z2 are area-level auxiliary variables;
kab_kota identifies the small area domain and the random-intercept group;
WEIND is the survey weight;
cluster_ids = ~1 indicates that the simulated example is treated as having no separate PSU clustering structure;
strata = "kab_kota" specifies the strata variable in the simulated survey design;
summary_function = "mean" requests domain mean estimation.

A concise summary of the fitted estimator can be displayed with:

summary(result)
#> SAE Projection Estimator using Linear Multilevel Model
#> -------------------------------------------------------
#> Formula   : Y ~ X1 + X2 + X3 + X4 + Z1 + Z2 + (1 | kab_kota) 
#> Estimator : bias_corrected 
#> Domains   : 50 
#> 
#> Model diagnostics:
#>   nobs        : 250 
#>   sigma       : 9.6444 
#>   ICC         : 0.9048 
#>   singular    : FALSE 
#>   convergence : OK 
#> 
#> Estimates:
#>  kab_kota  estimate variance       se       rse
#>         1  63.63811 28.51940 5.340356  8.391758
#>         2 123.57033 22.92302 4.787799  3.874554
#>         3  72.21099 24.03748 4.902803  6.789553
#>         4  89.15406 25.01544 5.001543  5.610001
#>         5 160.68935 12.41104 3.522931  2.192386
#>         6  27.48805 28.16499 5.307070 19.306825

Final Domain-Level Estimates

The final domain-level estimates are stored in:

head(result$estimates)
#>   kab_kota  estimate variance       se       rse
#> 1        1  63.63811 28.51940 5.340356  8.391758
#> 2        2 123.57033 22.92302 4.787799  3.874554
#> 3        3  72.21099 24.03748 4.902803  6.789553
#> 4        4  89.15406 25.01544 5.001543  5.610001
#> 5        5 160.68935 12.41104 3.522931  2.192386
#> 6        6  27.48805 28.16499 5.307070 19.306825

The output contains the following columns:

Column	Description
`kab_kota`	Domain identifier
`estimate`	Final bias-corrected projection estimate
`variance`	Plug-in variance of the final estimate
`se`	Standard error of the final estimate
`rse`	Relative standard error in percent

The estimates can also be extracted as an ordinary data frame.

estimates <- as.data.frame(result)

head(estimates)
#>   kab_kota  estimate variance       se       rse
#> 1        1  63.63811 28.51940 5.340356  8.391758
#> 2        2 123.57033 22.92302 4.787799  3.874554
#> 3        3  72.21099 24.03748 4.902803  6.789553
#> 4        4  89.15406 25.01544 5.001543  5.610001
#> 5        5 160.68935 12.41104 3.522931  2.192386
#> 6        6  27.48805 28.16499 5.307070 19.306825

Estimation Components

Detailed estimation components are stored in:

head(result$estimation_details)
#>   kab_kota estimate_synthetic variance_synthetic  correction
#> 1        1           64.51902           6.651059 -0.88090302
#> 2        2          123.47261           6.399579  0.09772008
#> 3        3           72.45427           5.914642 -0.24328336
#> 4        4           89.31018           6.307633 -0.15612183
#> 5        5          159.60585           6.707629  1.08349965
#> 6        6           29.12853           5.924376 -1.64048198
#>   variance_correction estimate_final variance_final se_final rse_final n_model
#> 1           21.868346       63.63811       28.51940 5.340356  8.391758       5
#> 2           16.523438      123.57033       22.92302 4.787799  3.874554       5
#> 3           18.122836       72.21099       24.03748 4.902803  6.789553       5
#> 4           18.707803       89.15406       25.01544 5.001543  5.610001       5
#> 5            5.703415      160.68935       12.41104 3.522931  2.192386       5
#> 6           22.240615       27.48805       28.16499 5.307070 19.306825       5
#>   n_proj
#> 1    300
#> 2    300
#> 3    300
#> 4    300
#> 5    300
#> 6    300

The table contains the following components:

Column	Description
`estimate_synthetic`	Synthetic projection estimate
`variance_synthetic`	Variance of the synthetic projection estimate
`correction`	Design-based residual correction
`variance_correction`	Variance of the residual correction
`estimate_final`	Final bias-corrected projection estimate
`variance_final`	Plug-in variance of the final estimate
`se_final`	Standard error of the final estimate
`rse_final`	Relative standard error of the final estimate
`n_model`	Number of model-survey observations in the domain
`n_proj`	Number of projection-survey observations in the domain

The synthetic component can be inspected separately.

head(
  result$estimation_details[, c(
    "kab_kota",
    "estimate_synthetic",
    "variance_synthetic"
  )]
)
#>   kab_kota estimate_synthetic variance_synthetic
#> 1        1           64.51902           6.651059
#> 2        2          123.47261           6.399579
#> 3        3           72.45427           5.914642
#> 4        4           89.31018           6.307633
#> 5        5          159.60585           6.707629
#> 6        6           29.12853           5.924376

The design-based residual correction can also be inspected separately.

head(
  result$estimation_details[, c(
    "kab_kota",
    "correction",
    "variance_correction"
  )]
)
#>   kab_kota  correction variance_correction
#> 1        1 -0.88090302           21.868346
#> 2        2  0.09772008           16.523438
#> 3        3 -0.24328336           18.122836
#> 4        4 -0.15612183           18.707803
#> 5        5  1.08349965            5.703415
#> 6        6 -1.64048198           22.240615

Model Parameters

Estimated model parameters are stored in result$model_parameters.

# Fixed-effect estimates
result$model_parameters$fixed_effects
#> (Intercept)          X1          X2          X3          X4          Z1 
#>   92.430911   25.107730   18.802686  -22.530931   20.487619   -9.648481 
#>          Z2 
#>   -6.124119

# Random-effect and residual variance components
result$model_parameters$variance_components
#>        grp        var1 var2     vcov     sdcor
#> 1 kab_kota (Intercept) <NA> 884.1652 29.734916
#> 2 Residual        <NA> <NA>  93.0136  9.644356

# Residual variance
result$model_parameters$residual_variance
#> [1] 93.0136

The estimated domain random effects can be accessed as follows.

head(result$model_parameters$random_effects$kab_kota)
#>   (Intercept)
#> 1  -41.868278
#> 2    4.644520
#> 3  -11.562970
#> 4   -7.420286
#> 5   51.497456
#> 6  -77.970168

Model Diagnostics

Model diagnostics are stored in:

result$diagnostics
#> $sigma
#> [1] 9.644356
#> 
#> $residual_variance
#> [1] 93.0136
#> 
#> $random_effects
#>        grp        var1 var2     vcov     sdcor
#> 1 kab_kota (Intercept) <NA> 884.1652 29.734916
#> 2 Residual        <NA> <NA>  93.0136  9.644356
#> 
#> $random_effect_groups
#> [1] "kab_kota"
#> 
#> $random_effect_dims
#> kab_kota 
#>        1 
#> 
#> $is_random_intercept_only
#> [1] TRUE
#> 
#> $icc
#> [1] 0.9048142
#> 
#> $icc_note
#> [1] "ICC computed for random-intercept structure."
#> 
#> $singular_fit
#> [1] FALSE
#> 
#> $convergence
#> [1] "OK"
#> 
#> $nobs
#> [1] 250
#> 
#> $REML
#> [1] TRUE
#> 
#> $logLik
#> [1] -1004.968
#> 
#> $AIC
#> [1] 2027.936
#> 
#> $BIC
#> [1] 2059.629

The intraclass correlation coefficient is reported only for a pure random-intercept structure. For random-slope or more complex random-effect structures, the simple ICC is returned as NA.

data.frame(
  icc = result$diagnostics$icc,
  singular_fit = result$diagnostics$singular_fit,
  convergence = result$diagnostics$convergence,
  sigma = result$diagnostics$sigma,
  residual_variance = result$diagnostics$residual_variance,
  REML = result$diagnostics$REML,
  AIC = result$diagnostics$AIC,
  BIC = result$diagnostics$BIC
)
#>         icc singular_fit convergence    sigma residual_variance REML      AIC
#> 1 0.9048142        FALSE          OK 9.644356           93.0136 TRUE 2027.936
#>        BIC
#> 1 2059.629

The fitted lmerMod object can be accessed directly.

fit <- result$fitted_model

summary(fit)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: Y ~ X1 + X2 + X3 + X4 + Z1 + Z2 + (1 | kab_kota)
#>    Data: data
#> Control: control
#> 
#> REML criterion at convergence: 2009.9
#> 
#> Scaled residuals: 
#>      Min       1Q   Median       3Q      Max 
#> -2.16963 -0.59140  0.04972  0.52613  2.44401 
#> 
#> Random effects:
#>  Groups   Name        Variance Std.Dev.
#>  kab_kota (Intercept) 884.17   29.735  
#>  Residual              93.01    9.644  
#> Number of obs: 250, groups:  kab_kota, 50
#> 
#> Fixed effects:
#>             Estimate Std. Error t value
#> (Intercept)  92.4309     4.4555  20.745
#> X1           25.1077     0.6709  37.422
#> X2           18.8027     1.3788  13.637
#> X3          -22.5309     0.6800 -33.133
#> X4           20.4876     0.5781  35.439
#> Z1           -9.6485     4.3391  -2.224
#> Z2           -6.1241     4.8901  -1.252
#> 
#> Correlation of Fixed Effects:
#>    (Intr) X1     X2     X3     X4     Z1    
#> X1 -0.003                                   
#> X2 -0.155 -0.023                            
#> X3 -0.019 -0.007  0.010                     
#> X4  0.005  0.003  0.051  0.109              
#> Z1  0.253 -0.018 -0.004 -0.007  0.006       
#> Z2 -0.054 -0.017  0.003 -0.012  0.024 -0.024

Residual diagnostics can be inspected using standard model diagnostic plots.

plot(
  fitted(fit),
  resid(fit),
  xlab = "Fitted values",
  ylab = "Residuals",
  main = "Residuals versus Fitted Values"
)

abline(h = 0, lty = 2)

qqnorm(resid(fit))
qqline(resid(fit))

The estimated random effects can also be inspected directly.

lme4::ranef(fit)
#> $kab_kota
#>     (Intercept)
#> 1  -41.86827839
#> 2    4.64451965
#> 3  -11.56297019
#> 4   -7.42028578
#> 5   51.49745648
#> 6  -77.97016757
#> 7    2.72074285
#> 8   29.35519933
#> 9  -42.62230470
#> 10 -27.79642104
#> 11   2.32279301
#> 12 -26.41673557
#> 13  -0.43544244
#> 14  -9.35890943
#> 15  12.60856337
#> 16 -52.39002472
#> 17  27.65965064
#> 18  42.23075634
#> 19 -26.95422097
#> 20  11.67685373
#> 21 -20.34097431
#> 22  15.81966705
#> 23 -59.78128098
#> 24 -45.21979483
#> 25  38.58000279
#> 26  24.74995717
#> 27  11.64353691
#> 28  -5.04989563
#> 29  28.08578562
#> 30   9.41398453
#> 31  13.18078825
#> 32  13.90819327
#> 33  12.24598079
#> 34  -6.91596274
#> 35 -31.06755394
#> 36   0.09122059
#> 37 -27.66702893
#> 38 -11.37982376
#> 39  -5.45123152
#> 40  33.12828614
#> 41  16.13556004
#> 42  54.29484285
#> 43 -12.48045973
#> 44  32.63175812
#> 45  18.41890554
#> 46   2.40170965
#> 47  10.37144283
#> 48  10.60677527
#> 49 -15.23594635
#> 50  34.96078070
#> 
#> with conditional variances for "kab_kota"

Retaining Unit-Level Predictions and Residuals

Set keep_unit = TRUE when unit-level predictions and model residuals are needed.

The following code is shown for illustration and is not evaluated when the vignette is built because it refits the same model.

result_unit <- sae_ml_linear(
  formula = Y ~ X1 + X2 + X3 + X4 + Z1 + Z2 + (1 | kab_kota),
  data_model = saeml_modelsvy,
  data_proj = saeml_projsvy,
  domain = "kab_kota",
  cluster_ids = ~1,
  weight = "WEIND",
  strata = "kab_kota",
  summary_function = "mean",
  keep_unit = TRUE
)

head(result_unit$unit_projection)

head(result_unit$unit_model_residual)

When keep_unit = TRUE:

unit_projection contains the projection-survey data with an additional .prediction column;
unit_model_residual contains the model-survey data with additional .fitted_model and .model_residual columns.

Direct Estimation

Set return_direct = TRUE to calculate direct design-based estimates from the model survey.

The direct estimator is returned separately and does not replace the final projection estimator.

result_direct <- sae_ml_linear(
  formula = Y ~ X1 + X2 + X3 + X4 + Z1 + Z2 + (1 | kab_kota),
  data_model = saeml_modelsvy,
  data_proj = saeml_projsvy,
  domain = "kab_kota",
  cluster_ids = ~1,
  weight = "WEIND",
  strata = "kab_kota",
  summary_function = "mean",
  return_direct = TRUE
)

head(result_direct$direct_estimator)

Multiple Domain Variables

The domain argument can be supplied as a character scalar, a character vector, or a one-sided formula.

The package data include both prov and kab_kota. Therefore, both variables can be used jointly as domain identifiers.

result_multi <- sae_ml_linear(
  formula = Y ~ X1 + X2 + X3 + X4 + Z1 + Z2 + (1 | kab_kota),
  data_model = saeml_modelsvy,
  data_proj = saeml_projsvy,
  domain = c("prov", "kab_kota"),
  cluster_ids = ~1,
  weight = "WEIND",
  strata = "kab_kota",
  summary_function = "mean"
)

head(result_multi$estimates)

Survey Design Specification

The arguments cluster_ids, weight, and strata are passed to survey::svydesign().

For the simulated data used in this vignette, the survey design specification is:

cluster_ids = ~1
weight = "WEIND"
strata = "kab_kota"

Use cluster_ids = ~1 when observations are treated as unclustered for the survey-design specification. When the original survey has PSU clustering, supply the actual PSU identifier.

For a real survey with PSU clustering and stratification, supply the actual design variables.

result_clustered <- sae_ml_linear(
  formula = Y ~ X1 + X2 + X3 + X4 + Z1 + Z2 + (1 | kab_kota),
  data_model = data_model,
  data_proj = data_proj,
  domain = "kab_kota",
  cluster_ids = "psu_id",
  weight = "survey_weight",
  strata = "stratum",
  summary_function = "mean",
  nest = TRUE
)

In this specification:

psu_id identifies the primary sampling unit;
survey_weight identifies the sampling weight;
stratum identifies the sampling stratum;
nest = TRUE indicates that PSUs are nested within strata.

Important Data Requirements

Before running sae_ml_linear(), ensure that:

data_model contains the response variable and all required predictors.
data_proj contains the fixed-effect predictors, random-effect grouping variables, domain variables, and survey design variables.
The response variable is not required in data_proj.
Required variables do not contain missing values.
Fixed-effect categorical variables in data_proj do not contain levels that are absent from data_model.
New grouping levels for random effects are allowed in data_proj.
Survey-design variables have compatible names in both datasets.
Survey weights are appropriate for the requested mean or total.

For summary_function = "total", survey weights should be appropriate expansion weights for population totals. When summary_function = "mean", the weights are used to obtain weighted domain means.

When a domain occurs in data_proj but not in data_model, the residual correction is set to zero. Therefore, the final projection estimate for that domain is equal to its synthetic estimate.

Output Structure

The function returns an S3 object of class "sae_ml_linear".

names(result)
#> [1] "call"               "formula"            "estimator"         
#> [4] "fitted_model"       "model_parameters"   "estimates"         
#> [7] "estimation_details" "diagnostics"        "notes"

Important components of the output are:

Component	Description
`call`	Matched function call
`formula`	Final model formula after preprocessing
`estimator`	Estimator type
`fitted_model`	Fitted `lmerMod` object
`model_parameters`	Fixed effects, random effects, and variance components
`estimates`	Final domain-level estimates
`estimation_details`	Synthetic, correction, and final-estimation components
`diagnostics`	Model diagnostics
`notes`	Run-specific notes
`unit_projection`	Unit-level predictions when `keep_unit = TRUE`
`unit_model_residual`	Unit-level residual data when `keep_unit = TRUE`
`direct_estimator`	Direct estimates when `return_direct = TRUE`

References

Bates, D., Maechler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1), 1–48. https://doi.org/10.18637/jss.v067.i01

Finch, W. H., Bolin, J. E., & Kelley, K. (2014). Multilevel Modeling Using R. CRC Press.

Food and Agriculture Organization of the United Nations. (2021). Guidelines on Data Disaggregation for SDG Indicators Using Survey Data (1st ed.). https://doi.org/10.4060/cb3253en

Hox, J. J., Moerbeek, M., & van de Schoot, R. (2018). Multilevel Analysis: Techniques and Applications (3rd ed.). Routledge.

Kim, J. K., & Rao, J. N. K. (2012). Combining data from two independent surveys: A model-assisted approach. Biometrika, 99(1), 85–100. https://doi.org/10.1093/biomet/asr063

Moura, F. A. S., & Holt, D. (1999). Small area estimation using multilevel models. Survey Methodology, 25(1), 73–80. https://www150.statcan.gc.ca/n1/pub/12-001-x/1999001/article/4714-eng.pdf

Rao, J. N. K., & Molina, I. (2015). Small Area Estimation (2nd ed.). Wiley.