Smoothing Sampling Variances in Python with svy-sae

Generalized variance functions for Fay–Herriot inputs

Smooth direct sampling variances with a generalized variance function in Python before fitting a Fay-Herriot model: Gamma and log-linear fits, precision weights, leverage, and leave-one-out checks with svy-sae.

Keywords

generalized variance function Python, variance smoothing small area estimation, sampling variance Fay-Herriot, GVF survey Python, smooth_variance svy-sae

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.

The Fay–Herriot model treats the sampling variance \(D_i\) of each direct estimate as known (Rao and Molina 2015). In practice it is itself an estimate, and a poor one where the sample is small: a district sampled in a few EAs has a variance estimate with few degrees of freedom, and a district sampled in one EA, or with no poor household in the sample, gets a variance of zero.

Those errors carry straight into the model. Each estimate is a weighted average of the direct estimate and the regression prediction, with weight \(\gamma_i = \sigma^2_u / (\sigma^2_u + D_i)\) on the direct estimate. A variance that is too small pulls the estimate toward the direct one; a variance of zero keeps it there, with an MSE of zero.

A generalized variance function (GVF) replaces the direct variances with values predicted from a model of how the variance depends on domain characteristics, such as the sample size. Every district then gets a usable variance, including those whose direct variance cannot be used.

This tutorial demonstrates how to smooth the district variances with smooth_variance(), check the variance model, and see what smoothing changes in a Fay–Herriot fit.

Setting up the inputs

We continue with the district poverty rates and census covariates built in Preparing inputs:

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")
)

A variance model

The variance of an estimated proportion is roughly \(\text{deff} \cdot p(1-p)/n\), so on the log scale it falls linearly with \(\log n\). Add the log of the number of sampled households:

direct = direct.with_columns(pl.col("n").log().alias("log_n"))

Fit the variance model and add the smoothed variances:

smoothed = sae.smooth_variance(direct, x="log_n", variance="var_direct", df="df")

smoothed.select("geo2", "n", "df", "var_direct", "variance_smoothed").head()

shape: (5, 5)

geo2	n	df	var_direct	variance_smoothed
str	i64	f64	f64	f64
"geo_01_01"	250	9.0	0.001872	0.00636
"geo_01_02"	50	1.0	0.00324	0.024909
"geo_01_03"	300	11.0	0.000216	0.005448
"geo_01_04"	75	2.0	0.005641	0.01766
"geo_01_05"	100	3.0	0.006911	0.013836

smooth_variance() takes any frame with one row per domain and returns it with one more column, variance_smoothed (set by name=). Pass se= instead of variance= to smooth the square of a standard error column.

The model is fitted on the districts whose direct variance carries information: a positive variance with positive df. The others are left out of the fit and predicted, which is what the four districts from Preparing inputs need:

smoothed.filter(pl.col("var_direct") < 1e-12).select(
    "geo2", "est", "n", "df", "var_direct", "variance_smoothed"
)

shape: (4, 6)

geo2	est	n	df	var_direct	variance_smoothed
str	f64	i64	f64	f64	f64
"geo_02_01"	0.12	25	0.0	1.2037e-35	0.044846
"geo_09_01"	0.0	25	0.0	0.0	0.044846
"geo_10_01"	0.0	50	1.0	0.0	0.024909
"geo_10_07"	0.44	25	0.0	8.7750e-33	0.044846

Their smoothed variances are those of any district with the same sample size.

Why df matters

df= weights each district by its design degrees of freedom, the precision of its variance estimate, and leaves out the districts with df of zero. Without it, the two single-EA districts enter the fit with variances of about \(10^{-33}\), which are rounding error, not zero:

unweighted = sae.smooth_variance(
    direct, x="log_n", variance="var_direct", method="loglinear"
)

unweighted.select(
    pl.col("variance_smoothed").min().alias("smallest"),
    pl.col("variance_smoothed").max().alias("largest"),
)

shape: (1, 2)

smallest	largest
f64	f64
6.2917e-12	3862.474767

On the log scale they sit about 70 units below the other districts and drag the fit with them.

WarningAlways pass df

svy’s estimation output has a df column. Pass it as df= so that districts whose variance is not estimable are predicted rather than fitted.

Choosing a method

The method argument sets how the model is fitted:

Table 1: Variance smoothing methods

method	Fit	Smoothed variance
"gamma" (default)	Gamma GLM with log link, weighted by df	\(\exp(z'\hat\alpha)\), the mean of \(\hat v\)
"loglinear"	Least squares on \(\log \hat v\), weighted by df	\(\exp(z'\hat\alpha)\), the median of \(\hat v\)
"loglinear_ratio"	As "loglinear"	Scaled so the smoothed variances sum to the direct ones

The Gamma model fits the mean of the variance directly, which is what Fay–Herriot needs. The log-linear fit back-transforms to the median, which falls below the mean of a skewed quantity such as a variance estimate. Compare the totals over the districts used in the fit:

fitted = (pl.col("var_direct") > 1e-12) & (pl.col("df") > 0)

pl.concat(
    sae.smooth_variance(direct, x="log_n", variance="var_direct", df="df", method=m)
    .filter(fitted)
    .select(
        pl.lit(m).alias("method"),
        (pl.col("variance_smoothed").sum() / pl.col("var_direct").sum()).alias(
            "smoothed / direct"
        ),
    )
    for m in ["gamma", "loglinear", "loglinear_ratio"]
)

shape: (3, 2)

method	smoothed / direct
str	f64
"gamma"	1.167081
"loglinear"	0.477454
"loglinear_ratio"	1.0

The log-linear median sums to less than half of the direct variances; the ratio version restores the total by construction.

Checking the variance model

diagnostics=True adds two columns for each district used in the fit:

• variance_smoothed_leverage: its hat value, how much its own variance shapes its prediction,
• variance_smoothed_loo: its smoothed variance when the model is refitted without it.

smoothed = sae.smooth_variance(
    direct, x="log_n", variance="var_direct", df="df", diagnostics=True
)

smoothed.select(
    "geo2", "n", "var_direct", "variance_smoothed",
    "variance_smoothed_leverage", "variance_smoothed_loo",
).sort("variance_smoothed_leverage", descending=True, nulls_last=True).head()

shape: (5, 6)

geo2	n	var_direct	variance_smoothed	variance_smoothed_leverage	variance_smoothed_loo
str	i64	f64	f64	f64	f64
"geo_04_02"	475	0.003331	0.003689	0.327544	0.003921
"geo_01_03"	300	0.000216	0.005448	0.093362	0.006069
"geo_02_08"	300	0.007071	0.005448	0.084875	0.005276
"geo_03_01"	300	0.007745	0.005448	0.084875	0.005204
"geo_10_08"	300	0.002139	0.005448	0.084875	0.005801

The district with the most leverage is the one with the largest sample, alone at that end of \(\log n\). Its smoothed variance barely changes without it, so the fit does not depend on it. The leave-one-out values also compare variance models out of sample: a model that predicts each district’s direct variance well without that district is not just fitting noise. The districts left out of the fit have no leverage, and their leave-one-out value is their smoothed variance, since they never entered the fit.

What smoothing changes

Fit the same Fay–Herriot model, with one covariate, the share of households with piped water (chosen in Choosing covariates), once with the direct variances and once with the smoothed ones, and compare both with the true rates:

def fh_rmse(variance):
    fit = sae.AreaLevel(smoothed).fh(
        y="est", x=["piped_water"], area="geo2", variance=variance,
        pop_data=district_x,
    )
    pred = fit.to_polars("predictions").join(truth, left_on="area", right_on="geo2")
    return pred, pred.select(((pl.col("pred") - pl.col("true_rate")) ** 2).mean().sqrt()).item()


pred_direct, rmse_direct = fh_rmse("var_direct")
pred_smoothed, rmse_smoothed = fh_rmse("variance_smoothed")

print(f"RMSE against the truth, direct variances:   {rmse_direct:.3f}")
print(f"RMSE against the truth, smoothed variances: {rmse_smoothed:.3f}")

RMSE against the truth, direct variances:   0.141
RMSE against the truth, smoothed variances: 0.033

With the direct variances, the four districts keep their direct estimates, with an MSE of zero:

problem = ["geo_02_01", "geo_09_01", "geo_10_01", "geo_10_07"]

pred_direct.filter(pl.col("area").is_in(problem)).select("area", "pred", "mse", "true_rate")

shape: (4, 4)

area	pred	mse	true_rate
str	f64	f64	f64
"geo_02_01"	0.12	1.2037e-35	0.572237
"geo_09_01"	0.0	0.0	0.412918
"geo_10_01"	0.0	0.0	0.423086
"geo_10_07"	0.44	8.7750e-33	0.238306

With the smoothed variances, the model moves them toward the regression prediction and reports a real MSE:

pred_smoothed.filter(pl.col("area").is_in(problem)).select("area", "pred", "mse", "true_rate")

shape: (4, 4)

area	pred	mse	true_rate
str	f64	f64	f64
"geo_02_01"	0.425325	0.001051	0.572237
"geo_09_01"	0.372003	0.000715	0.412918
"geo_10_01"	0.362376	0.000779	0.423086
"geo_10_07"	0.211066	0.000332	0.238306

Summary

This tutorial smoothed the district sampling variances before the Fay–Herriot fit:

1. A variance model with smooth_variance(): log variance on log sample size, weighted by df
2. Unusable variances (zero, or from a single EA) left out of the fit and predicted
3. Methods: "gamma" for the mean of the variance, "loglinear" for its median, "loglinear_ratio" to keep the total
4. Checks with diagnostics=True: leverage and leave-one-out variances

Next steps

The inputs are ready: a direct estimate and a usable variance for every sampled district, and covariates for every district.

Fit the model
Continue to The Fay–Herriot Model →

References

Rao, J. N. K., and I Molina. 2015. Small Area Estimation, 2nd edn. John Wiley & Sons, Hoboken, New Jersey.