Choosing Covariates for a Fay–Herriot Model in Python with svy-sae

Stepwise selection by AIC, BIC, or KIC

Choose the covariates of a Fay-Herriot small area model in Python with svy-sae: stepwise selection by AIC, BIC, or KIC, search direction, collinearity, and categorical candidates.

Keywords

covariate selection Fay-Herriot, stepwise selection small area estimation, AIC BIC KIC Python, variance inflation factor, area-level model selection Python

WarningAbout these examples

This tutorial illustrates the svy-sae API; it is not a rigorous small area analysis. The models are kept simple so each step stays readable, and in many places, from the choice of covariates to the variance model, better choices could be made to improve them.

A census offers many candidate covariates, and with a few dozen areas the model can afford only a few of them. Too many covariates fit the noise in the direct estimates, and census shares that measure the same thing, such as household assets, are close to collinear: their coefficients become unstable and the estimates with them.

Information criteria trade fit against the number of parameters. Stepwise selection searches the candidates one term at a time, dropping or adding whichever term lowers the criterion most, and stops when no move does.

This tutorial demonstrates how to choose the covariates of the Fay–Herriot model with AreaLevel.stepwise(), compare criteria and search directions, and check the choice against the true district poverty rates.

NoteNot a search for the best model

This tutorial does not try to find the best model for these data. It shows the options svy-sae offers to evaluate a model and the ways to improve it; the choices made below are examples, not recommendations.

Setting up the inputs

We continue with the district poverty rates, smoothed variances, and census covariates from the previous tutorials:

import polars as pl
import svy

import svy_sae as sae

hld_data = svy.datasets.load("hld_sample_wb_2023").data
hld_design = svy.Design(stratum=("geo1", "urbrur"), psu="ea", wgt="hhweight")
hld_sample = svy.Sample(data=hld_data, design=hld_design).wrangling.mutate(
    {
        "hhpovline": svy.col("hhsize") * 1800,
        "pov_status": svy.when(svy.col("tot_exp") < svy.col("hhpovline")).then(1).otherwise(0),
    }
)

direct = (
    hld_sample.estimation.prop(y="pov_status", by="geo2")
    .to_polars()
    .filter(pl.col("pov_status") == 1)
    .select("geo2", "est", "se", "cv", "lci", "uci", "df", "n")
    .with_columns((pl.col("se") ** 2).alias("var_direct"))
)

hld_census = svy.datasets.load("hld_pop_wb_2023").data
yes_no = ["electricity", "car", "refrigerator", "tv", "cell", "bank"]
district_x = (
    hld_census.group_by("geo1", "geo2")
    .agg(
        (pl.col("urbrur") == "Urban").mean().alias("urban"),
        (pl.col("piped_water") != "No piped water").mean().alias("piped_water"),
        (pl.col("flush_toilet") == "Flush toilet").mean().alias("flush_toilet"),
        *[(pl.col(c) == "Yes").mean().alias(c) for c in yes_no],
        pl.col("hhsize").mean().alias("hhsize"),
        pl.col("rooms").mean().alias("rooms"),
        pl.len().alias("n_hh"),
    )
    .sort("geo2")
)

truth = (
    hld_census.group_by("geo2")
    .agg((pl.col("tot_exp") < pl.col("hhsize") * 1800).mean().alias("true_rate"))
    .sort("geo2")
)

direct = direct.with_columns(pl.col("n").log().alias("log_n"))
smoothed = sae.smooth_variance(direct, x="log_n", variance="var_direct", df="df")

model = sae.AreaLevel(smoothed)

Why not use them all?

The census frame has eleven candidates:

candidates = [
    "urban", "electricity", "piped_water", "flush_toilet", "car",
    "refrigerator", "tv", "cell", "bank", "hhsize", "rooms",
]

Fit them all and look at the variance inflation factors, \(1 / (1 - R^2)\) of each covariate regressed on the others, reported with the model checks:

fit_all = model.fh(
    y="est", x=candidates, area="geo2", variance="variance_smoothed", pop_data=district_x
)

fit_all.to_polars("tests").filter(pl.col("name").str.starts_with("vif_")).select(
    "name", "statistic"
).sort("statistic", descending=True)

shape: (11, 2)

name	statistic
str	f64
"vif_refrigerator"	1827.312025
"vif_tv"	1432.602433
"vif_bank"	1416.866021
"vif_piped_water"	827.799686
"vif_flush_toilet"	679.906474
…	…
"vif_electricity"	353.184066
"vif_cell"	253.767611
"vif_car"	164.864742
"vif_rooms"	33.226871
"vif_hhsize"	8.197904

A VIF of 10 is a common warning level; most of these are in the hundreds. Districts with more refrigerators also have more TVs and bank accounts, so the model cannot tell their effects apart. The fit flags it as SAE_COLLINEARITY (see Checking the model).

Stepwise selection

stepwise() takes the same inputs as fh(), with the candidates as x:

selection = model.stepwise(
    y="est",
    x=candidates,
    area="geo2",
    variance="variance_smoothed",
    pop_data=district_x,
    criterion="bic",
)

print(selection)

Selected ['piped_water'] by BIC = -71.1715 in 10 steps, 59 areas

It returns a Selection. Its path records each move and the criterion after it:

selection.to_polars()

shape: (11, 3)

step	change	criterion
i64	str	f64
0	""	-38.012504
1	"-flush_toilet"	-42.051917
2	"-bank"	-46.045849
3	"-refrigerator"	-50.054227
4	"-tv"	-54.072561
…	…	…
6	"-cell"	-61.683902
7	"-car"	-65.146713
8	"-hhsize"	-67.983168
9	"-electricity"	-70.599501
10	"-urban"	-71.171537

With the default direction="both", the search starts from all the candidates; each step drops one term or adds back one already dropped. Here it drops ten in turn and keeps piped water, the covariate used in The Fay–Herriot model.

Candidates are compared by maximum likelihood fits without an MSE, so the search is fast. The criteria therefore differ from those that fh() reports for its REML fit.

Table 1: Options of stepwise()

Argument	Values	Description
criterion	"aic" (default), "bic", "kic"	Penalty per parameter: \(2\), \(\log m\) for \(m\) areas, or \(3\)
direction	"both" (default), "backward", "forward"	Moves allowed: drop and add, drop only, or add only
start	list of terms	Terms of the first model; by default all of x, or none with direction="forward"
vif_max	number	Reject an addition that would put any VIF above it

Comparing criteria

The criteria differ in how much a parameter costs. With 59 districts, BIC charges \(\log 59 \approx 4.1\) per parameter, more than AIC (2) or KIC (3):

pl.DataFrame(
    {
        "criterion": c,
        "terms": ", ".join(model.stepwise(
            y="est", x=candidates, area="geo2", variance="variance_smoothed",
            pop_data=district_x, criterion=c,
        ).terms),
    }
    for c in ["aic", "bic", "kic"]
)

shape: (3, 2)

criterion	terms
str	str
"aic"	"urban, piped_water"
"bic"	"piped_water"
"kic"	"urban, piped_water"

AIC and KIC also keep the share of urban households.

Search direction

A stepwise search is greedy: each step takes the best single move, so different starting points can end in different models. Search forward from no covariates instead:

forward = model.stepwise(
    y="est",
    x=candidates,
    area="geo2",
    variance="variance_smoothed",
    pop_data=district_x,
    criterion="bic",
    direction="forward",
)

print(forward)

Selected ['cell'] by BIC = -71.6490 in 1 steps, 59 areas

The forward search adds mobile phone ownership and stops. Its BIC is slightly lower than piped water’s. The first search never reached it: once down to piped water alone, adding mobile phones does not lower the BIC, and no single move swaps one term for another.

A difference of less than about 2 in an information criterion is commonly read as no real difference. Both covariates measure how well-off a district is; run both directions and choose between near-ties on substance, such as which covariate is measured more reliably in every district. We keep piped water.

Tipvif_max

vif_max=10 stops the search from adding a term that would push any VIF above 10. It only constrains additions: a backward search from a collinear start still reaches the collinear models on its way down.

Categorical candidates

A candidate can be a categorical variable, which enters and leaves the model as one term with all its levels. Add the region as a candidate:

with_region = model.stepwise(
    y="est",
    x=candidates + [svy.Cat("geo1")],
    area="geo2",
    variance="variance_smoothed",
    pop_data=district_x,
    criterion="bic",
)

with_region.to_polars().head(3)

shape: (3, 3)

step	change	criterion
i64	str	f64
0	""	-7.71644
1	"-geo1"	-38.012504
2	"-flush_toilet"	-42.051917

The region, nine dummies, is the first term dropped.

Using the selection

selection.terms is ready for fh(x=...):

fit = model.fh(
    y="est",
    x=selection.terms,
    area="geo2",
    variance="variance_smoothed",
    pop_data=district_x,
)

How close did we get?

Compare each candidate set with the true rates:

def rmse(terms):
    pred = model.fh(
        y="est", x=terms, area="geo2", variance="variance_smoothed", pop_data=district_x
    ).to_polars("predictions")
    scored = pred.join(truth, left_on="area", right_on="geo2")
    return scored.select(((pl.col("pred") - pl.col("true_rate")) ** 2).mean().sqrt()).item()


pl.DataFrame(
    {"covariates": name, "terms": len(terms), "rmse": rmse(terms)}
    for name, terms in [
        ("all eleven", candidates),
        ("AIC: urban, piped water", ["urban", "piped_water"]),
        ("BIC: piped water", ["piped_water"]),
        ("forward BIC: mobile phones", ["cell"]),
    ]
)

shape: (4, 3)

covariates	terms	rmse
str	i64	f64
"all eleven"	11	0.063236
"AIC: urban, piped water"	2	0.035372
"BIC: piped water"	1	0.032911
"forward BIC: mobile phones"	1	0.034543

The full model’s error is about twice that of any selected one. The selected models are close to each other, as their criteria were.

Summary

This tutorial chose the covariates of the Fay–Herriot model:

1. Collinearity in the full model, read from the VIFs in to_polars("tests")
2. Stepwise selection with stepwise(), which returns the chosen terms and the search path
3. Criteria: BIC penalizes parameters more than AIC or KIC with this many areas
4. Direction: a greedy search can end at different near-equivalent models; run more than one
5. Categorical candidates, such as svy.Cat("geo1"), enter and leave as one term

Next steps

The model is chosen. Before using its estimates, check its assumptions.

Check the model
Continue to Checking the Model →