How to use this
file
This is one complete, runnable notebook, in order. Run the chunks
from top to bottom. Do not skip a chunk, because later chunks use
objects created earlier — for example, fh() in Step 10
needs combined_data from Step 9, and
combined_data needs district_ratio from Step
6.
Every step has the same four parts, so you always
know where to look:
| 1. What this step does (blue box) |
Before the code |
The purpose, in plain words |
| 2. Code + comments |
Code block |
Runnable code, with a comment on every line |
| 3. Output |
Directly under the code |
What R prints when the code runs |
| 4. Explain the output (green box) |
After the output |
How to read the numbers and columns |
Data: all main examples use synthetic / example data
that come with the emdi package (eusilcA_smp
and eusilcA_popAgg). You do not need to download any file
for Steps 1–12.
Your own survey: only two places must change for
your own data:
- The survey design in Steps 3 and 5 — replace
ids,
strata, weights and data with
your own design. See Step 15.
- The poverty threshold in Step 4 — replace
10885.33 with
your own threshold.
Two small fixes made to the source notebooks, so this file
runs in order:
- The English notebook writes
count(district) and then
uses a column called smp_size. count()
actually creates a column called n, so that code would
fail. Step 8 fixes it with
count(district, name = "smp_size").
- The English notebook jumps straight to
fh() without
ever creating combined_data. Step 9 adds that missing step,
exactly as in the Bengali notebook:
combine_data(pop_data = eusilcA_popAgg, ...).
References from the source notebooks:
# Run this chunk ONLY once, if a package is not installed yet.
# eval = FALSE means this chunk is shown but NOT run when you knit the file,
# so knitting never tries to reinstall packages.
# install.packages("emdi") # Small Area Estimation
# install.packages("survey") # Complex survey analysis
# install.packages("tidyverse") # Data handling and plots (dplyr, ggplot2, ...)
Step 1 — Load the
packages
Load three packages. emdi does Small Area Estimation
(SAE) and gives uncertainty (MSE) by analytic and bootstrap methods.
survey handles complex survey data (stratification,
clustering, weighting). tidyverse is the toolbox for data
handling.
1.1 The
emdi package
# Load the emdi package into this R session.
# After this line, all emdi functions (for example fh() and combine_data())
# and its example data sets become available.
library(emdi)
1.2 The
survey package
survey main functions: svydesign() (make a
design), svymean() (weighted mean), svytotal()
(weighted total), svyglm() (weighted regression).
# Load the survey package into this R session.
# After this line, svydesign(), svyby() and svymean() become available.
library(survey)
1.3 The
tidyverse package
tidyverse includes ggplot2 (graphics),
dplyr (data manipulation), tidyr (reshaping),
readr (importing), purrr (functional
programming), tibble (modern data frames),
stringr (text) and forcats (categories).
# Load the whole tidyverse collection into this R session.
# After this line, mutate(), count(), group_by(), summarize() and the
# pipe operator |> become available.
library(tidyverse)
Loading packages normally prints nothing important. A warning such as
“Paket ‘tidyverse’ wurde unter R Version 4.4.1 erstellt” only says the
package was built under R 4.4.1. It is not an error —
your code still runs normally.
Step 2 — Explore the
data
eusilcA_smp is an emdi example data set
with living-standards data for households in EU countries. It is used
only to show the method. Each row is one household.
# Load the example sample data set from the emdi package into memory.
data(eusilcA_smp)
# Show the first 6 rows, so you can see the columns and the values.
head(eusilcA_smp)
# Show the structure: one line per column, with its type and first values.
str(eusilcA_smp)
#> 'data.frame': 1945 obs. of 18 variables:
#> $ eqIncome : num 23485 22704 25946 16153 22587 ...
#> $ gender : Factor w/ 2 levels "male","female": 1 1 2 1 1 1 2 1 2 1 ...
#> $ eqsize : num 2.5 2 1 2.1 1.5 2 1.5 2 1.8 1.5 ...
#> $ cash : num 33003 20073 0 19972 21503 ...
#> $ self_empl : num 0 0 24736 0 0 ...
#> $ unempl_ben : num 0 0 0 0 0 0 0 0 0 0 ...
#> $ age_ben : num 0 0 0 0 0 ...
#> $ surv_ben : num 0 0 0 0 0 0 0 0 0 0 ...
#> $ sick_ben : num 0 0 0 0 0 ...
#> $ dis_ben : num 0 0 0 0 0 0 0 0 0 0 ...
#> $ rent : num 0 0 0 0 0 0 0 0 0 0 ...
#> $ fam_allow : num 0 0 0 5247 2031 ...
#> $ house_allow: num 0 0 1084 0 0 ...
#> $ cap_inv : num 17158.8 4.7 126.3 195.3 19.9 ...
#> $ tax_adj : num -10973 -941 0 0 0 ...
#> $ state : chr "Burgenland" "Burgenland" "Burgenland" "Burgenland" ...
#> $ district : Factor w/ 70 levels "Neusiedl am See",..: 1 1 1 1 1 1 1 1 1 1 ...
#> $ weight : num 9.69 9.69 9.69 9.69 9.69 ...
eqIncome |
Equivalised income — income adjusted for household
size/structure |
gender |
Gender (male / female) |
eqsize |
Household size (equivalence scale) |
cash |
Cash income |
self_empl |
Self-employment income |
unempl_ben |
Unemployment benefits |
age_ben |
Age-related benefits |
surv_ben |
Survival benefits |
sick_ben |
Sickness benefits |
dis_ben |
Disability benefits |
rent |
Rent |
fam_allow |
Family allowances |
house_allow |
Housing allowances |
cap_inv |
Capital investment / income |
tax_adj |
Tax adjustment |
state |
State (Bundesland) |
district |
District — the small area / domain in this notebook |
weight |
Sampling weight — how many population units one sample unit stands
for |
Example — first row shown in the source notebook (row 213):
eqIncome 23,485.32, gender male,
eqsize 2.5, cash 33,003.39, state
Burgenland, district Neusiedl am See, weight
9.6875.
Step 3 — Uncertainty
measures for direct estimates
Question from the training: how can you easily get
several uncertainty measures — standard error (SE), coefficient of
variation (CV), CV in percent, variance and confidence interval (CI) —
for direct point estimates?
Answer: define the survey design once with
svydesign(), then use svyby() with
vartype = c("se", "cv", "cvpct", "var", "ci"). To get
measures for more variables, add them with + in the
formula. To change the estimate, change FUN (for example,
svytotal instead of svymean).
Adjust the design below to your own survey. This example uses
ids = ~1 (no clustering) only because the example data have
no cluster variable. See Step 15 for strata, clusters and finite
population corrections.
3.1 General pattern
(adapt this to your own data)
# eval = FALSE: this is a template only. It is NOT run, because your_data,
# variable1, variable2 and group_variable do not exist in this notebook.
# Load the survey package (already done in Step 1, repeated here for clarity).
library(survey)
# Define YOUR survey design. Change ids, weights and data to match your survey.
sampling_design <- svydesign(ids = ~1,
weights = ~weight,
data = your_data)
# Compute means by group, with all uncertainty measures.
# formula = the variables you want estimates for.
# by = the grouping variable (for example, district).
# FUN = the function that makes the point estimate (here, the mean).
# deff = TRUE asks for the design effect as well.
# vartype = which uncertainty measures to return.
results <- svyby(formula = ~variable1 + variable2,
by = ~group_variable,
design = sampling_design,
FUN = svymean,
deff = TRUE,
vartype = c("se", "cv", "cvpct", "var", "ci"))
# Print the results table.
print(results)
3.2 Worked example:
mean eqIncome and cash by district
Part A — Define the
design
# Define the survey design for the example data.
# ids = ~1 : each unit is its own cluster (no clustering, no strata).
# weights = ~weight : the 'weight' column holds the sampling weights.
# data = eusilcA_smp : the data set the design is built from.
sampling_design <- svydesign(ids = ~1,
weights = ~weight,
data = eusilcA_smp)
Part B — Estimate
the means with uncertainty
# Estimate the weighted mean of eqIncome and cash, separately for each district.
# formula = ~eqIncome + cash : estimate both variables (the + adds a variable).
# by = ~district : do it separately inside each district.
# design = sampling_design : use the design defined just above.
# FUN = svymean : the point estimate is the (weighted) mean.
# deff = TRUE : also compute the design effect (DEff).
# vartype = c(...) : also return se, cv, cvpct, var and ci.
district_means <- svyby(formula = ~eqIncome + cash,
by = ~district,
design = sampling_design,
FUN = svymean,
deff = TRUE,
vartype = c("se", "cv", "cvpct", "var", "ci"))
# Show the first 6 rows of the result.
head(district_means)
district |
District / area name |
eqIncome, cash |
Weighted mean for the district |
se.eqIncome, se.cash |
Standard error. Small SE = reliable mean; large SE = unreliable
mean |
ci_l.eqIncome, ci_u.eqIncome |
Lower and upper end of the (95%) confidence interval |
cv.eqIncome, cv.cash |
Coefficient of variation = SE / mean. Use it to compare areas with
very different means. Rule of thumb: CV up to about 10% is usually
acceptable |
cvpct.eqIncome, cvpct.cash |
Same CV, shown as a percentage |
var.eqIncome, var.cash |
Variance of the estimate |
DEff.eqIncome, DEff.cash |
Design effect. DEff = 1: design works like a simple random sample;
DEff > 1: design increased the variance; DEff < 1: it reduced
it |
Example from the source notebook — Neusiedl am See:
mean eqIncome 19,692.53, se.eqIncome 1,365.96,
ci_l.eqIncome 17,015.30, ci_u.eqIncome
22,369.77, cv.eqIncome 0.0694 (6.9%, reliable),
cv.cash 0.2105 (21.1%, less reliable),
DEff.eqIncome 1.046 (variance about 4.6% higher than under
a simple random sample).
Step 4 — Make the
poverty variable (for the Fay-Herriot model)
The standard Fay-Herriot (FH) model is made for means. For an
indicator between 0 and 1 — a proportion, such as a poverty rate — a
transformation is used. emdi offers two:
arcsine (arcsine square-root) and
logit. This notebook fits both and compares them in
Steps 10–12.
The arcsine transformation has been used, for example, for literacy
rates in Senegal from mobile-phone data, for poverty mapping in Chilean
comunas, and for regional unemployment in German functional urban
areas.
To show the method, we first make a binary poverty variable: 1 if a
household is poor, 0 otherwise.
4.1 Base R
version
# Set the poverty threshold (same value as in the source notebooks).
threshold <- 10885.33
# Make a new 0/1 column called 'poor' in the sample data.
# ifelse(test, 1, 0) gives 1 when the test is TRUE and 0 when it is FALSE.
# Here the test is: is equivalised income below the threshold?
eusilcA_smp$poor <- ifelse(eusilcA_smp$eqIncome < threshold, 1, 0)
# Illustrative only: the arcsine square-root transformation of that 0/1 variable.
# asin() is arcsine, sqrt() is square root. You do NOT need this column for fh():
# fh(transformation = "arcsin") does the transformation itself inside the model.
eusilcA_smp$transformed_poor <- asin(sqrt(eusilcA_smp$poor))
4.2 Tidyverse version
(run this too — the rest of the notebook uses it)
# Store the same poverty threshold under a short name, as in the source notebook.
npt <- 10885.33
# Take the sample data, then (|>) pass it into mutate().
# mutate() adds (or replaces) a column.
# eqIncome < npt is TRUE/FALSE; as.numeric() turns TRUE into 1 and FALSE into 0.
eusilcA_smp <- eusilcA_smp |>
mutate(poor = as.numeric(eqIncome < npt))
# Show the first 6 rows, including the new 'poor' column at the end.
head(eusilcA_smp)
The new column poor is 1 for a
household whose equivalised income is below the threshold (10,885.33)
and 0 otherwise. Section 4.2 overwrites
poor with the same values, in a cleaner way.
Step 5 — Define the
survey design again (now with poor)
Step 3 defined a design from the original data. The data now have a
poor column, so define the design again from the updated
data before using poor in Step 6.
# Define the survey design again, from the updated data that now contains 'poor'.
# ids = ~1 : each unit is its own cluster (no clustering).
# weights = ~weight : the 'weight' column holds the sampling weights.
# data = eusilcA_smp : the updated sample data, now including 'poor'.
sampling_design <- svydesign(ids = ~1,
weights = ~weight,
data = eusilcA_smp)
No table is printed — this step only creates the object
sampling_design. Without an error message, it worked.
Step 6 — District
proportions, with variance and design effect
The mean of a 0/1 variable is a proportion. So
svymean(poor) by district is the (weighted) poverty
proportion in each district.
# Estimate the poverty proportion by district.
# formula = ~poor : estimate only 'poor' (its mean = proportion poor).
# by = ~district : do it separately inside each district.
# design = sampling_design : use the design from Step 5.
# FUN = svymean : the point estimate is the (weighted) mean = proportion.
# deff = TRUE : also compute the design effect for 'poor'.
# vartype = c("var") : this time return only the variance (not se/cv/ci).
district_ratio <- svyby(formula = ~poor,
by = ~district,
design = sampling_design,
FUN = svymean,
deff = TRUE,
vartype = c("var"))
# Show the first 6 rows: district, proportion poor, variance, design effect.
head(district_ratio)
Reading the first rows: Neusiedl am See has
poor = 0.125 (12.5% poor), var = 0.00684,
DEff.poor = 1.046.
Bruck an der Leitha has poor = 0 and
DEff.poor = NaN. The design effect cannot be
computed when the proportion is 0, because there is no variation to
compare. That row is handled in Step 7.
Step 7 — Handle zero
estimates (out-of-sample domains)
Districts with a 0 estimate (and therefore NaN design
effect) cannot be used as in-sample domains. Drop them here; the FH
model will still give them a model-based (indirect) estimate later, as
out-of-sample domains. That is why the model output in Step 10 reports
out-of-sample domains.
# Keep only rows where DEff.poor is NOT missing (NA / NaN).
# !is.na(...) is TRUE for good rows and FALSE for the zero-estimate rows.
# district_ratio[...] keeps only the TRUE rows, so the bad rows are removed.
district_ratio <- district_ratio[!is.na(district_ratio$DEff.poor), ]
# Alternative used in the Bengali notebook (same idea, drops ANY row with NA):
# district_ratio <- district_ratio |> drop_na()
No table is printed. The rows with NaN design effect
(such as Bruck an der Leitha) are now removed from
district_ratio.
Step 8 — Effective
sample size by district
For the arcsine transformation, the variance on the transformed scale
is approximated with the effective sample size:
\[\text{effective sample size} =
\frac{\text{actual sample size}}{\text{design effect}}\]
8.1 Count the sample
size in each district
# Count how many sample households are in each district.
# count(district) groups by district and counts rows in each group.
# name = "smp_size" names the count column 'smp_size' (not the default 'n').
# FIX: the English notebook used count(district) alone, which makes a column
# called 'n', and then tried to use 'smp_size' — that would fail.
district_smp_size <- eusilcA_smp |>
count(district, name = "smp_size")
# Show the first 6 rows: district and its sample size.
head(district_smp_size)
8.2 Merge with the
proportion table
# Merge the proportion table with the sample-size table, by district.
# merge(..., by = "district") matches rows that share the same district.
# all = TRUE keeps districts that appear in either table (so out-of-sample
# districts are not silently lost at this stage).
combined_ratio <- merge(district_ratio, district_smp_size,
by = "district", all = TRUE)
# Show the first 6 rows of the merged table.
head(combined_ratio)
8.3 Compute the
effective sample size
# Add the effective sample size column to the merged table.
# mutate() creates a new column 'eff_smp_size' = sample size / design effect.
combined_ratio <- combined_ratio |>
mutate(eff_smp_size = smp_size / DEff.poor)
# Show the first 6 rows, now including 'eff_smp_size'.
head(combined_ratio)
# From here on, work with this table under the name used by the FH step.
district_ratio <- combined_ratio
Example: Neusiedl am See has sample size 16 and DEff
1.046, so its effective sample size is 16 / 1.046 =
15.30 — slightly less than 16, because the design is
slightly less efficient than a simple random sample. A smaller effective
sample size means less information, so a less certain direct
estimate.
Step 9 — Combine
sample data with population (aggregate) data
The FH model needs, for every district, both the direct estimate
(from the sample, Steps 6–8) and auxiliary variables from the population
(for example, mean cash and self_empl).
combine_data() joins those two sources by domain.
# Load the aggregated population data from the emdi package into memory.
data(eusilcA_popAgg)
# Combine population and sample data by domain.
# pop_data = eusilcA_popAgg : the population (aggregate) data.
# pop_domains = "Domain" : the domain column in the population data.
# smp_data = district_ratio : the sample estimates from Steps 6-8.
# smp_domains = "district" : the domain column in the sample data.
combined_data <- combine_data(
pop_data = eusilcA_popAgg,
pop_domains = "Domain",
smp_data = district_ratio,
smp_domains = "district"
)
# Show the first 6 rows of the combined data.
head(combined_data)
The message “Non-sampled domains exist” is expected.
It means some districts have population data but no usable direct
estimate (for example, the zero-estimate district dropped in Step 7). In
the table, such a district — for example Bruck an der Leitha — shows
NA for poor and var. The FH model
will estimate these districts indirectly, from the model.
Step 12 — Compare
arcsine and logit results
fh objects store their district estimates in
$ind, with columns Domain (district),
Direct (direct estimate), FH (model estimate)
and Out (TRUE for out-of-sample domains). Merging the two
$ind tables by district gives a clean side-by-side
comparison.
12.1 Side-by-side
table
# Build a small table with only the district name and the arcsine FH estimate.
# fh_arcsin$ind$Domain : district names from the arcsine model.
# fh_arcsin$ind$FH : FH (model-based) estimates from the arcsine model.
arcsin_table <- data.frame(district = fh_arcsin$ind$Domain,
FH_arcsin = fh_arcsin$ind$FH)
# Build the same small table for the logit model.
logit_table <- data.frame(district = fh_logit$ind$Domain,
FH_logit = fh_logit$ind$FH)
# Merge the two tables by district, so each row is one district,
# with both estimates side by side.
comparison <- merge(arcsin_table, logit_table, by = "district")
# Show the first rows of the comparison.
head(comparison)
# Alternative, documented helper (same estimates, different presentation):
# estimators() presents point, MSE and CV estimates for an emdi object.
head(estimators(fh_arcsin))
head(estimators(fh_logit))
Source-notebook examples: Neusiedl am See — arcsine
0.1235, logit 0.1337. Oberwart — arcsine 0.3507, logit
0.3863. Often the two are similar. When they differ, check diagnostics,
AIC/BIC and whether the result makes substantive sense, and prefer the
transformation with the better fit.
12.2 Summary of the
predictions
# Summarise the arcsine FH predictions across all districts.
# Min. / 1st Qu. / Median / Mean / 3rd Qu. / Max. describe the distribution.
summary(comparison$FH_arcsin)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.001128 0.082907 0.213241 0.212077 0.313476 0.726852
# Do the same for the logit predictions, for comparison.
summary(comparison$FH_logit)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 0.004382 0.094778 0.206847 0.228064 0.318746 0.691424
Source-notebook summary of the predictions: Min. 0.0011, 1st quartile
0.0829, Median 0.2132, Mean 0.2121, 3rd quartile 0.3135, Max. 0.7269.
Mean (0.212) close to median (0.213) means the predictions are roughly
symmetric — no strong skew.
Why a model-based estimate differs from the direct
estimate. Example from the source notebook: direct estimate
0.0455, FH estimate 0.1031 for the same district. This is expected, not
a mistake:
- Small samples make direct estimates unstable and highly
variable.
- A direct estimate reflects only that one sample.
- The model borrows strength across districts and uses predictors
(
cash, self_empl).
- EBLUP (the FH prediction) shrinks extreme direct estimates toward
the regression-predicted value, especially for small districts.
Step 13 — Model
comparison with AIC and BIC (illustrative example)
Lower AIC and lower BIC both mean a better trade-off between fit and
complexity. BIC penalises complexity more heavily than AIC, especially
with larger samples. For the two FH models above, read AIC/BIC from
summary(fh_arcsin) and summary(fh_logit) in
Steps 10–11.
The table below is the separate, illustrative example from the source
notebook — it is not computed from the FH models above.
It only shows how to read such a table.
# Make a small example table with three made-up models and their AIC/BIC.
# Model : model name. AIC : Akaike Information Criterion. BIC : Bayesian IC.
model_comparison <- data.frame(
Model = c("Model 1", "Model 2", "Model 3"),
AIC = c(1450, 1405, 1420),
BIC = c(1465, 1420, 1440)
)
# Print the example table.
print(model_comparison)
#> Model AIC BIC
#> 1 Model 1 1450 1465
#> 2 Model 2 1405 1420
#> 3 Model 3 1420 1440
Verdict: Model 2 is best here, because it has the
lowest AIC (1405) and the lowest BIC (1420) of the three.
Step 14 —
Multicollinearity and VIF (optional diagnostic)
VIF (Variance Inflation Factor) checks whether predictors are too
correlated with each other. If they are, their separate effects are hard
to tell apart and coefficients become unstable.
The source notebook only wrote a placeholder here
(vif_values <- ...). The code below computes VIF without
any extra package, from the in-sample rows of combined_data
only (out-of-sample rows have no poor value, so they are
left out of this diagnostic).
# Keep only in-sample rows: rows where 'poor' is not missing.
# is.na(combined_data$poor) is TRUE for out-of-sample rows; ! turns it around.
vif_data <- combined_data[!is.na(combined_data$poor), ]
# List the predictor names you want to check.
predictors <- c("cash", "self_empl")
# For each predictor, regress it on the OTHER predictors, then VIF = 1 / (1 - R-squared).
# sapply() runs the same small calculation once per predictor and collects the results.
vif_values <- sapply(predictors, function(v) {
# other predictors: all predictors except the one being checked now.
others <- setdiff(predictors, v)
# Build a formula such as cash ~ self_empl, then fit it on the in-sample data.
fit_vif <- lm(as.formula(paste(v, "~", paste(others, collapse = " + "))),
data = vif_data)
# R-squared of that fit: how much of this predictor the others explain.
r2 <- summary(fit_vif)$r.squared
# VIF formula. If R-squared is 0, VIF = 1 (no problem). If high, VIF is high.
1 / (1 - r2)
})
# Print the VIF for each predictor.
print(vif_values)
#> cash self_empl
#> 2.052035 2.052035
# If you prefer the car package, the equivalent is:
# library(car); car::vif(lm(poor ~ cash + self_empl, data = vif_data))
| VIF = 1 |
No correlation with the other predictors |
No action |
| VIF > 4 |
Possible multicollinearity |
Investigate; consider dropping a redundant predictor |
| VIF > 10 |
Serious multicollinearity |
Remove or combine correlated predictors |
Other diagnostics for a linear model in general: residual vs. fitted
plot, QQ-plot (normality), scale-location plot, and Cook’s distance
(influential observations). The FH summary in Steps 10–11 already gives
skewness, kurtosis and Shapiro-Wilk tests for residuals and random
effects.
Step 15 — Specifying
the survey design in detail
15.1 Main arguments
of svydesign()
All arguments are formulas naming columns in the data
data frame. Omitting fpc means sampling with replacement.
When clusters are numbered uniquely only within a stratum, use
nest = TRUE. For a multistage design, id and
fpc list the stages in order.
id |
Sampling units (PSUs; optionally later stages) |
Formula, for example ~dnum |
strata |
Strata |
Formula, for example ~stype |
weights |
Sampling weights |
Formula, for example ~pw |
fpc |
Finite population corrections |
Formula, for example ~fpc |
15.2 Example:
California API data (included in survey)
# Load the California API example data sets that come with the survey package.
data(api)
# Stratified sample: one PSU per unit, strata = school type, weights = pw, fpc = fpc.
# id = ~1 : no clustering.
# strata = ~stype : stratify by school type.
# weights = ~pw : sampling weights in column 'pw'.
# data = apistrat : the stratified sample data.
# fpc = ~fpc : finite population correction in column 'fpc'.
dstrat <- svydesign(id = ~1, strata = ~stype, weights = ~pw, data = apistrat, fpc = ~fpc)
# Cluster sample: clusters are districts (dnum), weights = pw, fpc = fpc.
# id = ~dnum : the cluster (PSU) is the district number.
dclus1 <- svydesign(id = ~dnum, weights = ~pw, data = apiclus1, fpc = ~fpc)
# Two-stage cluster sample: first stage dnum, second stage snum; two fpc columns.
# id = ~dnum + snum : stage 1 = district, stage 2 = school.
# fpc = ~fpc1 + fpc2 : one fpc column per stage, in the same order.
dclus2 <- svydesign(id = ~dnum + snum, fpc = ~fpc1 + fpc2, data = apiclus2)
# Print a short summary of the stratified design, to confirm it was built.
print(dstrat)
#> Stratified Independent Sampling design
#> svydesign(id = ~1, strata = ~stype, weights = ~pw, data = apistrat,
#> fpc = ~fpc)
The printout confirms the design object was built (it names the
design type, for example a stratified sampling design, and shows the
call used). The three designs show the three common cases:
stratified (dstrat), one-stage
cluster (dclus1) and two-stage
cluster (dclus2).
15.3 Example: NHANES
2017–2018 (file not supplied — do not run as-is)
This chunk is eval = FALSE on purpose: it needs the file
Data/nhanes1718_rmph.Rdata, which was not supplied with the
notebooks. Run it only after you have that file, in that relative
location.
Never delete rows before making the design. Always
use the full data set for svydesign(). To analyse a
subpopulation, make an indicator (like nomiss below) and
use subset(design, ...). Deleting rows first understates
the standard errors.
# Load the NHANES data file from disk into memory (file not supplied here).
load("Data/nhanes1718_rmph.Rdata")
# Count the rows. The source notebook reports 9,254 rows.
nrow(nhanes)
# Clean and recode, keeping ALL rows (see the warning above the chunk).
# %>% is the tidyverse pipe: it passes the data on the left into the next step.
# mutate() adds or changes columns.
# fct_collapse() merges factor levels: two Hispanic groups become "Hispanic",
# and Asian + Other/Multi become "Non-Hispanic Other".
# replace_na(WTSAF2YR, 0) replaces missing fasting weights with 0.
# complete.cases(...) is TRUE only for rows with no missing value in the listed
# columns; combined with WTSAF2YR > 0, it marks the analysis subpopulation.
nhanes.mod <- nhanes %>%
mutate(
race_eth = fct_collapse(RIDRETH3,
"Hispanic" = c("Mexican American", "Other Hispanic"),
"Non-Hispanic Other" = c("Non-Hispanic Asian", "Other/Multi")),
WTSAF2YR = replace_na(WTSAF2YR, 0),
nomiss = complete.cases(LBDGLUSI, BMXWAIST, smoker, RIDAGEYR,
RIAGENDR, race_eth, income) & WTSAF2YR > 0)
# Tell survey how to handle a stratum with only one PSU ("lonely PSU").
# "adjust" centres that PSU, instead of failing; the second option does the
# same for domain (subpopulation) analysis.
options(survey.lonely.psu = "adjust")
options(survey.adjust.domain.lonely = TRUE)
# Three designs, same strata and PSUs, different weights for different NHANES parts.
# strata = ~SDMVSTRA : NHANES strata. id = ~SDMVPSU : NHANES PSUs.
# nest = TRUE : PSU numbers are unique only within strata.
# weights differ: WTINT2YR = interview, WTMEC2YR = examination, WTSAF2YR = fasting.
design.INT <- svydesign(strata = ~SDMVSTRA, id = ~SDMVPSU, weights = ~WTINT2YR, nest = TRUE, data = nhanes.mod)
design.MEC <- svydesign(strata = ~SDMVSTRA, id = ~SDMVPSU, weights = ~WTMEC2YR, nest = TRUE, data = nhanes.mod)
design.FST <- svydesign(strata = ~SDMVSTRA, id = ~SDMVPSU, weights = ~WTSAF2YR, nest = TRUE, data = nhanes.mod)
# Subset the DESIGN, not the data set.
# subset(design, nomiss) keeps all design information but marks only 'nomiss'
# rows as in the subpopulation. This keeps standard errors correct.
design.FST.nomiss <- subset(design.FST, nomiss)
15.4 Design degrees
of freedom
Confidence intervals and p-values depend on the degrees of freedom.
For complex surveys, the recommended design degrees of freedom are:
\[\text{design DF} = \text{number of PSUs}
- \text{number of strata}\]
Most survey functions do not use design DF by default,
but you can ask for it with degf().
# eval = FALSE: this needs the NHANES designs from Step 15.3, which need the
# file that was not supplied. With those designs built, this line would run.
# degf(design) computes the design degrees of freedom for that design.
# degf(design.FST)
Step 16 — Final
combined results: get the FH estimates out
The FH estimates from Steps 10–11 are in $ind. Export
that table instead of copying numbers by hand: one row per district,
with the direct estimate and both model-based estimates.
# Take the arcsine model's district table, which already has Direct and FH.
# fh_arcsin$ind$Domain : district. $Direct : direct estimate. $FH : FH estimate.
# $Out : TRUE if the district is out-of-sample (no direct estimate).
final_results <- data.frame(
district = fh_arcsin$ind$Domain,
Direct = fh_arcsin$ind$Direct,
FH_arcsin = fh_arcsin$ind$FH,
Out_of_sample = fh_arcsin$ind$Out
)
# Add the logit FH estimate, matched by district (merge is safer than
# assuming the two models list districts in the same order).
final_results <- merge(final_results,
data.frame(district = fh_logit$ind$Domain,
FH_logit = fh_logit$ind$FH),
by = "district")
# Show the first 6 rows of the final table.
head(final_results)
# Optional: write it to a CSV file in the working directory, for use elsewhere.
# write.csv(final_results, "Nairobi_FH_results.csv", row.names = FALSE)
Example rows from the source notebook (for checking, not computed
here):
| Eisenstadt-Umgebung |
0.0140 |
0.0320 |
| Güssing |
0.1926 |
0.1535 |
| Jennersdorf |
0.5493 |
0.3563 |
| Neusiedl am See |
0.125 (Step 6) |
arcsine 0.1235 / logit 0.1337 |
NA in the Direct column means no direct
estimate was available for that district — it was estimated by the model
only (Out_of_sample = TRUE).
Step 17 — Other
resources (no code — reading only)
Poverty mapping (World Bank): Poverty and
inequality on the map in the Gambia; Poverty mapping in Bhutan
— small area estimations: Gewog and town level results;
Hardship at local area council level in Vanuatu. More examples:
UNSD Wiki, or search “poverty mapping with small area estimation”.
Fay-Herriot applications to study: SAIPE poverty
estimates (US Census Bureau); State and Metro Area Employment, Hours
& Earnings (US Bureau of Labor Statistics); estimation of
unemployment (Statistics Canada); prevalence of international migration;
regional CPIH expenditure weights. Institutions use slightly different
FH adjustments — studying them shows why the model is used, not
just how.
Other SAE packages: sae (R, broader
model choice), hbsae (R, Bayesian), Stata sae,
Stata Fay-Herriot extensions fayherriot and
FHSAE. Check the package documentation for the extension
you need.
R Markdown help: RStudio R Markdown Cheat Sheet
(also shows how to suppress warnings/messages in output) and a general
Markdown Guide.
End of notebook. To run it: open it in RStudio, then Knit, or run
the chunks top to bottom in order. If a chunk fails, check first that
every earlier chunk ran — most errors in this notebook mean an earlier
object (sampling_design, district_ratio,
combined_data) is missing.