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

August 6, 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", source="bundled")
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 wealth quintile. Note the two categorical variables: urbrur uses the default reference (first alphabetically), while quint_nat specifies an explicit reference with ref:

lin_model = hld_sample.glm.fit(
    y="tot_exp",
    x=[
        "hhsize",
        "rooms",
        svy.Cat("urbrur"),
        svy.Cat("quint_nat", ref="poorest"),
    ],
    family="gaussian",
)

print(lin_model)

╭─────────────────────────────────── GLM: Gaussian (identity) ────────────────────────────────────╮
│ Modeling: tot_exp                                                                               │
│                                                                                                 │
│ Observations         825  AIC           16088.7063                                              │
│ DF Residuals          19  BIC                    -                                              │
│ Deviance      1.3468e+10  Scale         5.1801e+08                                              │
│ R-squared        0.72978  R-sq (adj)       0.72746                                              │
│                           Iterations             2                                              │
│ F-stat (adj)    42.89804  Prob (F-adj)      <0.001                                              │
│                                                                                                 │
│                                                                                                 │
│  Term                      Coef.     Std.Err.          t    P>|t|         [0.025        0.975]  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  _intercept_         -9617.12355   1228.91295   -7.82572   <0.001   -12189.26793   -7044.97918  │
│  hhsize               2394.70155    172.17753   13.90833   <0.001     2034.32983    2755.07327  │
│  rooms                1110.54046    155.90909    7.12300   <0.001      784.21899    1436.86194  │
│  urbrur_Urban          247.12744    614.28772    0.40230   0.6920    -1038.59153    1532.84641  │
│  quint_nat_average    6001.53950    766.99616    7.82473   <0.001     4396.19808    7606.88091  │
│  quint_nat_poor       3497.05840    700.86701    4.98962   <0.001     2030.12689    4963.98991  │
│  quint_nat_rich       9521.99935    797.47609   11.94017   <0.001     7852.86272   11191.13598  │
│  quint_nat_richest   16170.30178   1063.42595   15.20586   <0.001    13944.52570   18396.07787  │
│                                                                                                 │
╰─────────────────────────────────────────────────────────────────────────────────────────────────╯

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: (8, 8)

term	estimate	std_err	conf_low	conf_high	statistic	p_value	df
str	f64	f64	f64	f64	f64	f64	i64
"_intercept_"	-9617.123552	1228.912954	-12189.267926	-7044.979177	-7.825716	2.3209e-7	19
"hhsize"	2394.70155	172.177532	2034.329833	2755.073266	13.908328	2.0633e-11	19
"rooms"	1110.540462	155.90909	784.218987	1436.861937	7.123	8.9891e-7	19
"urbrur_Urban"	247.127441	614.287717	-1038.591527	1532.846409	0.402299	0.691954	19
"quint_nat_average"	6001.539496	766.996162	4396.19808	7606.880913	7.824732	2.3252e-7	19
"quint_nat_poor"	3497.058403	700.867009	2030.126894	4963.989913	4.989618	0.000081	19
"quint_nat_rich"	9521.999353	797.476086	7852.862722	11191.135984	11.940169	2.8190e-10	19
"quint_nat_richest"	16170.301784	1063.425946	13944.525699	18396.077868	15.205856	4.3296e-12	19

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 wealth quintile as categorical predictors:

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

print(logit_model)

╭───────────────────────────────── GLM: Binomial (logit) ─────────────────────────────────╮
│ Modeling: pov_status                                                                    │
│                                                                                         │
│ Observations         825  AIC           126.9280                                        │
│ DF Residuals          19  BIC                  -                                        │
│ Deviance        120.7610  Scale           1.0000                                        │
│ R-squared        0.85627  R-sq (adj)     0.85503                                        │
│                           Iterations          20                                        │
│ F-stat (adj)  3615.01748  Prob (F-adj)    <0.001                                        │
│                                                                                         │
│                                                                                         │
│  Term                    Coef.   Std.Err.           t    P>|t|      [0.025      0.975]  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  _intercept_         -23.30984    1.35363   -17.22020   <0.001   -26.14303   -20.47666  │
│  hhsize                0.22773    0.14016     1.62478   0.1207    -0.06563     0.52109  │
│  rooms                 0.01661    0.16771     0.09903   0.9222    -0.33441     0.36762  │
│  urbrur_Rural          1.15890    0.64411     1.79923   0.0879    -0.18924     2.50704  │
│  quint_nat_poor       22.60247    0.85662    26.38578   <0.001    20.80956    24.39539  │
│  quint_nat_poorest    42.84121    0.79260    54.05115   <0.001    41.18227    44.50015  │
│  quint_nat_rich        0.37281    1.25830     0.29628   0.7702    -2.26084     3.00647  │
│  quint_nat_richest     0.82403    0.95051     0.86693   0.3968    -1.16541     2.81347  │
│                                                                                         │
╰─────────────────────────────────────────────────────────────────────────────────────────╯

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: (8, 8)

term	estimate	std_err	conf_low	conf_high	statistic	p_value	df
str	f64	f64	f64	f64	f64	f64	i64
"_intercept_"	-23.309844	1.353633	-26.143031	-20.476657	-17.220204	4.7449e-13	19
"hhsize"	0.227729	0.14016	-0.06563	0.521088	1.624775	0.120684	19
"rooms"	0.016608	0.167708	-0.334408	0.367624	0.099028	0.922154	19
"urbrur_Rural"	1.158901	0.64411	-0.189238	2.507039	1.799227	0.087887	19
"quint_nat_poor"	22.602472	0.856616	20.809555	24.395389	26.385781	1.9642e-16	19
"quint_nat_poorest"	42.841211	0.792605	41.182269	44.500152	54.051149	2.8592e-22	19
"quint_nat_rich"	0.372812	1.258302	-2.260844	3.006469	0.296282	0.770229	19
"quint_nat_richest"	0.824028	0.950509	-1.165411	2.813467	0.866933	0.396795	19

Fitting on a Subpopulation with `where`

Pass where to fit the model on a subpopulation — here, households without electricity. As with subpopulation estimation elsewhere in svy, the full design is retained so standard errors stay correct (unlike pre-filtering the data):

logit_model_domain = hld_sample.glm.fit(
    y="pov_status",
    x=[
        "hhsize",
        "rooms",
        svy.Cat("urbrur", ref="Urban"),
    ],
    where = svy.col("electricity") == "No",
    family="binomial",
    link="logit",
)

print(logit_model_domain)

╭──────────────────────────── GLM: Binomial (logit) ─────────────────────────────╮
│ Modeling: pov_status                                                           │
│                                                                                │
│ Observations      115  AIC           97.0879                                   │
│ DF Residuals        8  BIC                 -                                   │
│ Deviance      90.4262  Scale          1.0000                                   │
│ R-squared     0.31679  R-sq (adj)    0.29832                                   │
│                        Iterations          6                                   │
│ F-stat (adj)  3.82528  Prob (F-adj)   0.0762                                   │
│                                                                                │
│                                                                                │
│  Term              Coef.   Std.Err.          t    P>|t|     [0.025     0.975]  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  _intercept_    -3.28986    1.13204   -2.90614   0.0197   -5.90034   -0.67937  │
│  hhsize          0.89074    0.25804    3.45191   0.0087    0.29569    1.48579  │
│  rooms          -0.18832    0.22723   -0.82874   0.4313   -0.71232    0.33568  │
│  urbrur_Rural    1.83275    0.68713    2.66726   0.0285    0.24823    3.41726  │
│                                                                                │
╰────────────────────────────────────────────────────────────────────────────────╯

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               825  DF           19.0 │
│ Mean ŷ       0.2370  Mean SE    0.0115 │
│ Min ŷ        0.0000  Max ŷ      1.0000 │
│ Mean resid  -0.0006  Std resid  0.1637 │
│                                        │
│ 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
2.8443e-10	1.1072e-10	1.2593e-10	6.4240e-10	-2.8443e-10
2.9896e-10	1.5053e-10	1.0421e-10	8.5764e-10	-2.9896e-10
3.6315e-10	9.0998e-11	2.1494e-10	6.1357e-10	-3.6315e-10
1.8417e-10	1.1729e-10	4.8566e-11	6.9842e-10	-1.8417e-10
2.3030e-10	1.1165e-10	8.3482e-11	6.3530e-10	-2.3030e-10
3.6923e-10	1.1142e-10	1.9633e-10	6.9440e-10	-3.6923e-10
2.8919e-10	1.0398e-10	1.3626e-10	6.1378e-10	-2.8919e-10
0.506409	0.163123	0.207466	0.800839	-0.506409
1.8417e-10	1.1729e-10	4.8566e-11	6.9842e-10	-1.8417e-10
2.8919e-10	1.0398e-10	1.3626e-10	6.1378e-10	-2.8919e-10

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.2149
Mean predicted probability (actual):         0.2370

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.181461   0.016286   [0.147375, 0.215548]  │
│   3.00   0.194937   0.009584   [0.174877, 0.214996]  │
│   5.00   0.206859   0.006826   [0.192572, 0.221147]  │
│   7.00   0.216674   0.008117   [0.199684, 0.233664]  │
│  10.00   0.227244   0.009368   [0.207636, 0.246852]  │
│                                                      │
╰──────────────────────────────────────────────────────╯

pred_margins.to_polars()

shape: (5, 6)

term	margin	se	lci	uci	value
str	f64	f64	f64	f64	i64
"hhsize"	0.181461	0.016286	0.147375	0.215548	1
"hhsize"	0.194937	0.009584	0.174877	0.214996	3
"hhsize"	0.206859	0.006826	0.192572	0.221147	5
"hhsize"	0.216674	0.008117	0.199684	0.233664	7
"hhsize"	0.227244	0.009368	0.207636	0.246852	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.005554   0.003173   [-0.001087, 0.012194]  │
│                                               │
╰───────────────────────────────────────────────╯

╭────── GLM Margins: rooms (ame, 95% CI) ───────╮
│                                               │
│    Margin         SE                  95% CI  │
│  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━  │
│  0.000405   0.004082   [-0.008139, 0.008949]  │
│                                               │
╰───────────────────────────────────────────────╯

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("quint_nat")]

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("quint_nat", ref="poorest")]

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, learn how to hand results downstream — as stable, versioned JSON payloads that report templates, QA rules, and pipelines can bind to.

Ship the results
Continue to Serializing Results →

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

Fitting on a Subpopulation with where

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

Fitting on a Subpopulation with `where`

Categorical Variables with `svy.Cat()`

Distribution Family (`family`)

Link Function (`link`)