Transformed Fay–Herriot Models in Python with svy-sae

Log, arcsine, and logit transformations with back-transformed estimates

Fit Fay-Herriot small area models on a transformed scale in Python with svy-sae: log for skewed outcomes, arcsine and logit for proportions, back-transformations, and Slud-Maiti or bootstrap MSE on the original scale.

Keywords

transformed Fay-Herriot model, log transformation small area estimation, arcsine transformation proportions, Slud-Maiti MSE, back-transformation bias correction

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 assumes normal area effects and sampling errors on the scale of the direct estimates. Skewed outcomes such as expenditure or income break that assumption, and so can proportions, which are bounded by 0 and 1. Fitting the model on a transformed scale, where the assumptions hold better, and transforming the estimates back is the usual remedy.

fh() does both steps: it moves the direct estimates and their sampling variances to the model scale, fits there, and returns estimates and MSEs on the original scale.

This tutorial demonstrates the log transformation on district mean expenditure, the choice of back-transformation and MSE, and the arcsine and logit transformations for proportions, applied to the 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 use the same survey, districts, and census as 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")
)

The outcome here is the mean per-capita expenditure of each district, with the true district means from the census:

expenditure = (
    hld_sample.estimation.mean(y="pc_exp", by="geo2")
    .to_polars()
    .select("geo2", "est", "se", "cv", "lci", "uci", "df", "n")
    .with_columns(
        (pl.col("se") ** 2).alias("var_direct"),
        pl.col("n").log().alias("log_n"),
        pl.col("est").log().alias("log_est"),
    )
)

true_expenditure = hld_census.group_by("geo2").agg(pl.col("pc_exp").mean().alias("true_mean"))

The variance of a mean falls with the sample size and grows with the level of the outcome, so the variance model of Smoothing variances uses both:

expenditure = sae.smooth_variance(
    expenditure, x=["log_n", "log_est"], variance="var_direct", df="df"
)

model = sae.AreaLevel(expenditure)

A skewed outcome

Fit the model on the original scale first, with mobile phone ownership as the covariate:

fit_raw = model.fh(
    y="est", x=["cell"], area="geo2", variance="variance_smoothed", pop_data=district_x
)

fit_raw.to_polars("findings").select("code", "title", "detail")

shape: (3, 3)

code	title	detail
str	str	str
"SAE_NONNORMAL_RESIDUALS"	"Residuals not normal"	"Shapiro-Wilk W = 0.959, p = 0.…
"SAE_NONNORMAL_AREA_EFFECTS"	"Area effects not normal"	"Shapiro-Wilk W = 0.886, p = 4.…
"SAE_INFLUENTIAL_AREAS"	"3 influential areas"	"Cook's distance above 0.0678: …

Both normality checks fail: the area effects and residuals are skewed, as expenditure is.

The log transformation

Pass transformation="log". There is no analytical MSE for these estimates with districts outside the sample, so request the bootstrap:

fit_log = model.fh(
    y="est",
    x=["cell"],
    area="geo2",
    variance="variance_smoothed",
    pop_data=district_x,
    transformation="log",
    mse="bootstrap",
    n_reps=1000,
    rstate=12345,
)

fit_log.to_polars("findings").select("code", "title")

shape: (1, 2)

code	title
str	str
"SAE_INFLUENTIAL_AREAS"	"3 influential areas"

The normality findings are gone. The model is fitted to the log of the direct estimates, with sampling variances \(D_i / \hat\theta_i^2\) by the delta method, and the estimates and MSEs come back in currency units:

fit_log.to_polars("predictions").select("area", "pred", "mse", "cv").head()

shape: (5, 4)

area	pred	mse	cv
str	f64	f64	f64
"geo_01_01"	4601.745855	87699.371756	0.064354
"geo_01_02"	4805.161994	119466.849775	0.071931
"geo_01_03"	4551.053808	78561.347154	0.061587
"geo_01_04"	4386.12568	94649.88017	0.070142
"geo_01_05"	4798.70828	104188.358667	0.067264

The per-district diagnostics keep both scales: the direct estimate and variance on the original scale, and on the model scale with a _t suffix.

Mobile phone ownership was chosen on the log scale: stepwise() takes the same transformation= argument as fh() and compares the candidates on that scale.

Back-transformation

Transforming the model-scale estimate back with \(\exp\) gives the median, not the mean, of the district’s expenditure given the data. The backtransform option of LogTransform sets how to return to the original scale:

Table 1: Back-transformations for LogTransform

backtransform	Estimate	MSE
"bc" (default)	\(\exp(\tilde\theta_i + s^2_i / 2)\), the mean given the data	Slud–Maiti with ML (Slud and Maiti 2006), or bootstrap
"naive"	\(\exp(\tilde\theta_i)\), the median	Bootstrap
"crude"	\(\exp(\tilde\theta_i + \text{mse}_i / 2)\)	Analytical only

Here \(s^2_i\) is the variance of the district’s true value given the data, on the log scale. Compare the bias-corrected and naive estimates:

fit_naive = model.fh(
    y="est",
    x=["cell"],
    area="geo2",
    variance="variance_smoothed",
    pop_data=district_x,
    transformation=sae.LogTransform(backtransform="naive"),
    mse="bootstrap",
    n_reps=1000,
    rstate=12345,
)

fit_log.to_polars("predictions").select("area", "pred").join(
    fit_naive.to_polars("predictions").select("area", pl.col("pred").alias("naive")), on="area"
).select(((pl.col("pred") / pl.col("naive") - 1) * 100).alias("bc above naive (%)")).describe()

shape: (9, 2)

statistic	bc above naive (%)
str	f64
"count"	61.0
"null_count"	0.0
"mean"	0.139803
"std"	0.00626
"min"	0.118898
"25%"	0.137154
"50%"	0.141239
"75%"	0.143587
"max"	0.15157

The correction is small here because the log-scale variances are small.

Analytical MSE

The Slud–Maiti MSE of the bias-corrected estimate needs ML estimation, and covers sampled areas only. Fit the sampled districts by ML without pop_data, so the covariate must be in the frame of direct estimates:

sampled = sae.AreaLevel(expenditure.join(district_x.select("geo2", "cell"), on="geo2"))

fit_sm = sampled.fh(
    y="est",
    x=["cell"],
    area="geo2",
    variance="variance_smoothed",
    transformation="log",
    method="ml",
)

fit_sm_boot = sampled.fh(
    y="est",
    x=["cell"],
    area="geo2",
    variance="variance_smoothed",
    transformation="log",
    method="ml",
    mse="bootstrap",
    n_reps=1000,
    rstate=12345,
)

(
    fit_sm.to_polars("predictions").select("area", "pred", pl.col("mse").alias("Slud–Maiti"))
    .join(fit_sm_boot.to_polars("predictions").select("area", pl.col("mse").alias("bootstrap")), on="area")
    .join(true_expenditure, left_on="area", right_on="geo2")
    .select(
        pl.col("Slud–Maiti", "bootstrap").mean(),
        ((pl.col("pred") - pl.col("true_mean")) ** 2).mean().alias("squared error"),
    )
)

shape: (1, 3)

Slud–Maiti	bootstrap	squared error
f64	f64	f64
84483.109203	37638.163098	44557.43785

Here the two disagree, and the bootstrap is the closer to the actual squared error. The Slud–Maiti MSE is a second-order approximation, and its term for the uncertainty of \(\hat\sigma^2_u\) is large when \(\hat\sigma^2_u\) is poorly determined, as here:

print(fit_sm.component("area"))

VarianceComponent(name='area', kind=variance, est=0.00147477, se=0.00678549)

Its standard error is several times the estimate. With few areas or a weak area variance, compare the two MSEs before choosing.

Each transformation has its own MSE options:

Table 2: Transformations for fh()

Transformation	Domain	Sampling variance on the model scale	Analytical MSE
"log", LogTransform()	\(\hat\theta > 0\)	\(D / \hat\theta^2\)	"bc" with ML, sampled areas only; "crude"
"arcsin", ArcsinTransform()	\(0 \le \hat\theta \le 1\)	\(D / (4\hat\theta(1-\hat\theta))\), or \(1/(4\tilde n)\) with eff_size=	None: bootstrap
"logit", LogitTransform()	\(0 < \hat\theta < 1\)	\(D / (\hat\theta(1-\hat\theta))^2\)	None: bootstrap

A name gives the transformation with its defaults; the classes carry its options.

How close did we get?

def scores(name, frame):
    scored = frame.join(true_expenditure, left_on="area", right_on="geo2")
    covered = (pl.col("true_mean") >= pl.col("lci")) & (pl.col("true_mean") <= pl.col("uci"))
    return scored.select(
        pl.lit(name).alias("estimates"),
        ((pl.col("pred") - pl.col("true_mean")) ** 2).mean().sqrt().alias("rmse"),
        pl.col("mse").mean().alias("mean mse"),
        ((pl.col("pred") - pl.col("true_mean")) ** 2).mean().alias("mean squared error"),
        covered.sum().alias("interval covers truth"),
    )


direct_pred = expenditure.select(
    pl.col("geo2").alias("area"), pl.col("est").alias("pred"),
    pl.col("var_direct").alias("mse"), "lci", "uci",
)

pl.concat([
    scores("direct", direct_pred),
    scores("FH, original scale", fit_raw.to_polars("predictions")),
    scores("FH, log scale", fit_log.to_polars("predictions")),
])

shape: (3, 5)

estimates	rmse	mean mse	mean squared error	interval covers truth
str	f64	f64	f64	u32
"direct"	836.336742	465721.648082	699459.14565	49
"FH, original scale"	196.543123	116894.399574	38629.199201	61
"FH, log scale"	216.141917	56377.887537	46717.328402	59

Both models cut the error of the direct estimates by about three quarters. On the original scale, the error is slightly lower, but the model’s assumptions fail and its MSE is three times the actual squared error. On the log scale, the assumptions hold and the MSE is close to the actual error: the uncertainty the model reports can be trusted.

Proportions

Poverty rates are bounded by 0 and 1, and their sampling variance depends on the rate itself. Two transformations address that: the arcsine of the square root, which stabilizes the variance, and the logit, which maps \((0, 1)\) to the whole line.

Two districts have a direct poverty rate of 0, which the logit cannot take:

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

try:
    poverty.fh(
        y="est", x=["piped_water"], area="geo2", variance="variance_smoothed",
        pop_data=district_x, transformation="logit", mse="bootstrap",
    )
except Exception as error:
    print(error)

  ❌ Direct estimates outside the transformation's domain [MODEL_ERROR]
  The logit needs direct estimates strictly between 0 and 1; areas ['geo_09_01', 'geo_10_01'] are not.
  Hint: ArcsinTransform(eff_size=...) takes estimates of 0 and 1.
  Docs: https://svylab.com/docs/svy/reference/codes.html#model-error

The arcsine takes 0 and 1, but its delta-method variance is undefined there. ArcsinTransform(eff_size=...) replaces it with \(1/(4\tilde n)\), from a column of effective sample sizes \(\tilde n = n / \text{deff}\). A district with no poor household has no design effect of its own, so use the national one:

deff = (
    hld_sample.estimation.prop(y="pov_status", deff="wr")
    .to_polars()
    .filter(pl.col("pov_status") == 1)["deff"]
    .item()
)
print(f"National design effect: {deff:.2f}")

poverty = sae.AreaLevel(smoothed.with_columns((pl.col("n") / deff).alias("n_eff")))

fit_arcsin = poverty.fh(
    y="est",
    x=["piped_water"],
    area="geo2",
    variance="variance_smoothed",
    pop_data=district_x,
    transformation=sae.ArcsinTransform(eff_size="n_eff"),
    mse="bootstrap",
    n_reps=1000,
    rstate=12345,
)
fit_rate = poverty.fh(
    y="est", x=["piped_water"], area="geo2", variance="variance_smoothed", pop_data=district_x
)

def rate_scores(name, fit):
    scored = fit.to_polars("predictions").join(truth, left_on="area", right_on="geo2")
    covered = (pl.col("true_rate") >= pl.col("lci")) & (pl.col("true_rate") <= pl.col("uci"))
    return scored.select(
        pl.lit(name).alias("estimates"),
        ((pl.col("pred") - pl.col("true_rate")) ** 2).mean().sqrt().alias("rmse"),
        covered.sum().alias("interval covers truth"),
    )

pl.concat([rate_scores("original scale", fit_rate), rate_scores("arcsine", fit_arcsin)])

National design effect: 8.92

shape: (2, 3)

estimates	rmse	interval covers truth
str	f64	u32
"original scale"	0.032911	58
"arcsine"	0.037283	38

For these rates the arcsine does not help. The area variance stays at the boundary, as in The Fay–Herriot model, and only the bootstrap MSE exists for the arcsine, which understates the MSE at the boundary: its intervals miss the truth far more often. The untransformed model, with its analytical MSE, remains the better choice for the poverty rates.

TipWhen to transform

Transform when the checks of the untransformed model fail in a way the transformation addresses, such as skewed residuals, and check that they pass afterwards. A transformation is a change of model, not a refinement that is always worth making.

Summary

This tutorial fitted the Fay–Herriot model on transformed scales:

1. Log for a skewed outcome: transformation="log" fixes the normality findings, and its MSE matches the actual error
2. Back-transformations: "bc" for the mean given the data, "naive" for the median, "crude" with the analytical MSE
3. MSE: Slud–Maiti for the log with ML on sampled areas, which overstates when \(\hat\sigma^2_u\) is poorly determined; otherwise the bootstrap
4. Proportions: the logit needs rates strictly between 0 and 1; the arcsine takes 0 and 1 with eff_size=
5. When not to: for the poverty rates, the untransformed model is better

Next steps

The district estimates are final. Regional figures already published from the survey may need to agree with them.

Benchmark the estimates
Continue to Benchmarking →

References

Slud, Eric V., and Tapabrata Maiti. 2006. “Mean-Squared Error Estimation in Transformed Fay-Herriot Models.” Journal of the Royal Statistical Society: Series B 68 (2): 239–57.