Generalized Linear Models (GLMs) for Complex Surveys in Python

Linear, logistic, and count regression with design-adjusted standard errors

Tutorials

GLM

Regression

Python

Fit linear, logistic, and Poisson regression models to complex survey data in Python. Learn design-adjusted GLMs with proper standard errors, categorical predictors, predictions, and marginal effects using the svy library.

Author

Mamadou S. Diallo, Ph.D.

Published

January 18, 2026

Modified

April 18, 2026

Keywords

survey weighted regression Python, logistic regression survey weights Python, linear regression complex survey Python, GLM complex survey Python, design-adjusted standard errors Python, Poisson regression survey data Python, survey regression analysis Python, binomial GLM survey Python, categorical predictors survey regression, Taylor linearization GLM Python, marginal effects survey Python, predictive margins survey Python

Generalized Linear Models (GLMs) extend ordinary linear regression to accommodate response variables with non-normal error distributions, such as binary, categorical, or count data.

In complex survey analysis, fitting these models requires special attention. While point estimates (coefficients) are computed using weighted estimating equations, standard variance estimation methods (like OLS) generally underestimate uncertainty because they assume observations are independent and identically distributed (i.i.d.).

The svy GLM module fits common regression models — linear (Gaussian), logistic (Binomial), Poisson, and Gamma — while correctly estimating standard errors using the survey design information (stratification, clustering, and weighting).

Setting Up the Sample

We’ll use the World Bank (2023) synthetic sample data:

import polars as pl
import svy

# Load data and define design
hld_data = svy.datasets.load(name="hld_sample_wb_2023")
hld_design = svy.Design(stratum=("geo1", "urbrur"), psu="ea", wgt="hhweight")
hld_sample = svy.Sample(data=hld_data, design=hld_design)

# Create a binary poverty status variable for logistic regression examples
hld_sample = hld_sample.wrangling.mutate(
    {
        "hhpovline": svy.col("hhsize") * 1800,
        "pov_status": svy.when(svy.col("tot_exp") < svy.col("hhpovline")).then(1).otherwise(0),
    }
)

Linear Regression

Linear regression is used when the outcome variable is continuous. In survey analysis, this is equivalent to solving weighted least squares, but with variance estimates that account for the complex design.

Estimate a model predicting total household expenditure from household size, number of rooms, area type, and tenure status. Note the two categorical variables: urbrur uses the default reference (first alphabetically), while statocc specifies an explicit reference with ref:

lin_model = hld_sample.glm.fit(
    y="tot_exp",
    x=[
        "hhsize",
        "rooms",
        svy.Cat("urbrur"),
        svy.Cat("statocc", ref="Occupied for free"),
    ],
    family="gaussian",
)

print(lin_model)

╭─────────────────────────────── GLM: Gaussian (identity) ────────────────────────────────╮
│ Modeling: tot_exp                                                                       │
│                                                                                         │
│ Observations        8000  AIC           2.5590e+11                                      │
│ DF Residuals         301  BIC           2.5590e+11                                      │
│ Deviance      2.5590e+11  Scale         8.5017e+08                                      │
│ R-squared        0.40869  R-sq (adj)       0.40832                                      │
│                           Iterations             2                                      │
│ F-stat (adj)   229.38407  Prob (F-adj)      <0.001                                      │
│                                                                                         │
│                                                                                         │
│  Term                  Coef.    Std.Err.          t    P>|t|       [0.025       0.975]  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  _intercept_        99.46746   387.26665    0.25684   0.7975   -662.62549    861.56041  │
│  hhsize            818.04073    55.05433   14.85879   <0.001    709.70061    926.38086  │
│  rooms            1957.89110   146.07358   13.40346   <0.001   1670.43634   2245.34587  │
│  urbrur_Urban     4833.23715   291.01783   16.60804   <0.001   4260.54999   5405.92432  │
│  statocc_Owned     555.87593   277.89565    2.00030   0.0464      9.01161   1102.74026  │
│  statocc_Rented    281.24103   276.31716    1.01782   0.3096   -262.51703    824.99909  │
│                                                                                         │
╰─────────────────────────────────────────────────────────────────────────────────────────╯

The output includes estimated coefficients, design-based standard errors, t-statistics, and confidence intervals. The model-level statistics include the residual deviance, AIC, and a Wald F-test for overall significance.

Exporting Results

Export coefficients to a Polars DataFrame:

lin_model.to_polars()

shape: (6, 8)

term	estimate	std_err	conf_low	conf_high	statistic	p_value	df
str	f64	f64	f64	f64	f64	f64	i64
"_intercept_"	99.467458	387.266648	-662.62549	861.560407	0.256845	0.797474	301
"hhsize"	818.040733	55.054329	709.70061	926.380855	14.85879	7.7458e-38	301
"rooms"	1957.891105	146.073579	1670.436336	2245.345873	13.403458	1.9154e-32	301
"urbrur_Urban"	4833.237153	291.017833	4260.549989	5405.924317	16.608045	2.0316e-44	301
"statocc_Owned"	555.875933	277.89565	9.011606	1102.740261	2.000305	0.046365	301
"statocc_Rented"	281.241032	276.317163	-262.517027	824.999092	1.01782	0.309581	301

Logistic Regression

Logistic regression is used when the outcome variable is binary (0/1), such as whether a household is below the poverty line. It models the log-odds of the outcome as a linear combination of the predictors.

Model the likelihood of a household being poor using household size and rooms as continuous predictors, and urban/rural and tenure as categorical predictors:

logit_model = hld_sample.glm.fit(
    y="pov_status",
    x=[
        "hhsize",
        "rooms",
        svy.Cat("urbrur", ref="Urban"),
        svy.Cat("statocc"),
    ],
    family="binomial",
    link="logit",
)

print(logit_model)

╭────────────────────────────── GLM: Binomial (logit) ──────────────────────────────╮
│ Modeling: pov_status                                                              │
│                                                                                   │
│ Observations       8000  AIC           823.6660                                   │
│ DF Residuals        301  BIC           865.5892                                   │
│ Deviance       811.6660  Scale           1.0000                                   │
│ R-squared       0.42494  R-sq (adj)     0.42458                                   │
│                          Iterations           6                                   │
│ F-stat (adj)  115.17864  Prob (F-adj)    <0.001                                   │
│                                                                                   │
│                                                                                   │
│  Term                Coef.   Std.Err.           t    P>|t|     [0.025     0.975]  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  _intercept_      -4.14685    0.26206   -15.82421   <0.001   -4.66255   -3.63115  │
│  hhsize            0.72939    0.03994    18.26359   <0.001    0.65080    0.80798  │
│  rooms            -0.68916    0.05974   -11.53561   <0.001   -0.80672   -0.57159  │
│  urbrur_Rural      1.99915    0.13429    14.88630   <0.001    1.73487    2.26342  │
│  statocc_Owned     0.23676    0.14548     1.62745   0.1047   -0.04953    0.52304  │
│  statocc_Rented   -0.09581    0.18105    -0.52918   0.5971   -0.45210    0.26048  │
│                                                                                   │
╰───────────────────────────────────────────────────────────────────────────────────╯

Note

Interpretation: The coefficients are on the log-odds scale. A positive coefficient indicates that the predictor increases the probability of the outcome. To interpret them as odds ratios (OR), exponentiate the coefficients (\(e^\beta\)).

logit_model.to_polars()

shape: (6, 8)

term	estimate	std_err	conf_low	conf_high	statistic	p_value	df
str	f64	f64	f64	f64	f64	f64	i64
"_intercept_"	-4.146852	0.262057	-4.662549	-3.631155	-15.824206	1.8402e-41	301
"hhsize"	0.729387	0.039937	0.650797	0.807978	18.263587	1.1296e-50	301
"rooms"	-0.689155	0.059742	-0.806719	-0.571591	-11.535613	9.7120e-26	301
"urbrur_Rural"	1.999147	0.134294	1.734872	2.263422	14.8863	6.1121e-38	301
"statocc_Owned"	0.23676	0.145479	-0.049525	0.523045	1.627447	0.104689	301
"statocc_Rented"	-0.09581	0.181053	-0.4521	0.26048	-0.529184	0.597068	301

Prediction

The predict() method computes fitted values with confidence intervals on the response scale. For logistic models, predictions are on the probability scale (the inverse-logit is applied automatically):

preds = logit_model.predict(hld_sample.data, y_col="pov_status")

print(preds)

╭─────── GLM Predictions (95% CI) ───────╮
│ n              8000  DF          301.0 │
│ Mean ŷ       0.2293  Mean SE    0.0140 │
│ Min ŷ        0.0000  Max ŷ      1.0000 │
│ Mean resid  -0.0002  Std resid  0.3183 │
│                                        │
│ Use .to_polars() for full results      │
╰────────────────────────────────────────╯

Export predictions to a DataFrame:

preds.to_polars().head(10)

shape: (10, 5)

yhat	se	lci	uci	residuals
f64	f64	f64	f64	f64
0.020435	0.003774	0.01419	0.029346	-0.020435
0.014739	0.002888	0.010013	0.021646	-0.014739
0.005443	0.000997	0.003795	0.007802	-0.005443
0.041468	0.006453	0.030474	0.0562	-0.041468
0.003756	0.000827	0.002434	0.005791	-0.003756
0.010785	0.001661	0.007961	0.014595	-0.010785
0.008073	0.001644	0.005404	0.012044	-0.008073
0.010785	0.001661	0.007961	0.014595	-0.010785
0.014739	0.002888	0.010013	0.021646	-0.014739
0.035682	0.005905	0.025718	0.049311	-0.035682

Prediction on New Data

You can also predict on new or counterfactual data. The new data must contain all predictor columns used in the model:

# Counterfactual: what if all households were urban with 4 rooms?
new_data = hld_sample.data.with_columns(
    pl.lit("Urban").alias("urbrur"),
    pl.lit(4).alias("rooms"),
)

preds_cf = logit_model.predict(new_data)
print(f"Mean predicted probability (counterfactual): {preds_cf.yhat.mean():.4f}")
print(f"Mean predicted probability (actual):         {preds.yhat.mean():.4f}")

Mean predicted probability (counterfactual): 0.0615
Mean predicted probability (actual):         0.2293

Marginal Effects and Predictive Margins

After fitting a model, you often want to understand the practical impact of each predictor — not just whether it’s statistically significant. The margins() method provides two complementary views:

Predictive margins (at): predicted values at specific levels of a variable, averaging over the rest of the covariates
Average marginal effects (variables): the average change in the outcome for a unit change in a predictor

Predictive Margins

Compute the predicted probability of poverty at different household sizes, averaging over all other variables in the model:

pred_margins = logit_model.margins(at={"hhsize": [1, 3, 5, 7, 10]})

print(pred_margins)

╭────── GLM Margins: hhsize (predictive, 95% CI) ──────╮
│                                                      │
│  Value     Margin         SE                 95% CI  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│   1.00   0.034043   0.000488   [0.033084, 0.035003]  │
│   3.00   0.115738   0.001482   [0.112821, 0.118655]  │
│   5.00   0.287209   0.002825   [0.281649, 0.292769]  │
│   7.00   0.519985   0.003330   [0.513433, 0.526538]  │
│  10.00   0.832527   0.002176   [0.828244, 0.836809]  │
│                                                      │
╰──────────────────────────────────────────────────────╯

pred_margins.to_polars()

shape: (5, 6)

term	margin	se	lci	uci	value
str	f64	f64	f64	f64	i64
"hhsize"	0.034043	0.000488	0.033084	0.035003	1
"hhsize"	0.115738	0.001482	0.112821	0.118655	3
"hhsize"	0.287209	0.002825	0.281649	0.292769	5
"hhsize"	0.519985	0.00333	0.513433	0.526538	7
"hhsize"	0.832527	0.002176	0.828244	0.836809	10

Average Marginal Effects (AME)

Compute the average marginal effect of each continuous predictor:

ame = logit_model.margins()

for m in ame:
    print(m)

╭───── GLM Margins: hhsize (ame, 95% CI) ──────╮
│                                              │
│    Margin         SE                 95% CI  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  0.074366   0.004072   [0.066353, 0.082379]  │
│                                              │
╰──────────────────────────────────────────────╯

╭─────── GLM Margins: rooms (ame, 95% CI) ────────╮
│                                                 │
│     Margin         SE                   95% CI  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  -0.070264   0.006091   [-0.082251, -0.058278]  │
│                                                 │
╰─────────────────────────────────────────────────╯

The AME tells you: on average, how much does a one-unit increase in the predictor change the predicted probability? Unlike log-odds coefficients, AMEs are on the probability scale and directly interpretable.

Understanding GLM Parameters

Categorical Variables with `svy.Cat()`

Categorical variables need special handling to create indicator (dummy) variables. Wrap them in svy.Cat():

x=[svy.Cat("urbrur"), svy.Cat("statocc")]

By default, the first category alphabetically serves as the reference. To choose a different reference, use ref:

x=[svy.Cat("urbrur", ref="Urban"), svy.Cat("statocc", ref="Owned")]

The resulting coefficients are interpreted relative to the reference category. Mixing both styles (with and without ref) in the same model is common — use ref only when the default alphabetical reference isn’t the most interpretable baseline.

Distribution Family (`family`)

Family	Use Case
`"gaussian"`	Continuous outcomes (Linear Regression)
`"binomial"`	Binary (0/1) outcomes (Logistic Regression)
`"poisson"`	Count data (Poisson Regression)
`"gamma"`	Positively skewed continuous data

Link Function (`link`)

Connects the linear predictors (Xβ) to the expected mean of the distribution (μ):

Link	Description	Default for
`"identity"`	No transformation (μ = Xβ)	Gaussian
`"logit"`	Log of the odds: ln(p/(1−p)) = Xβ	Binomial
`"log"`	Log link: ln(μ) = Xβ	Poisson
`"inverse"`	Reciprocal: 1/μ = Xβ	Gamma

When link is omitted, the canonical link for the family is used automatically.

Next Steps

Now that you’ve covered weighting, estimation, and modeling, explore how to visualize results or export them for reporting.

Explore more tutorials
Return to Tutorials Overview →

References

World Bank. 2023. “Synthetic Data for an Imaginary Country, Sample, 2023.” World Bank, Development Data Group. https://doi.org/10.48529/MC1F-QH23.

Setting Up the Sample

Linear Regression

Exporting Results

Logistic Regression

Prediction

Prediction on New Data

Marginal Effects and Predictive Margins

Predictive Margins

Average Marginal Effects (AME)

Understanding GLM Parameters

Categorical Variables with svy.Cat()

Distribution Family (family)

Link Function (link)

Next Steps

References

Categorical Variables with `svy.Cat()`

Distribution Family (`family`)

Link Function (`link`)