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),
}
)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.
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:
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 1.6485e+07 │ │ 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 access to electricity as categorical predictors:
logit_model = hld_sample.glm.fit(
y="pov_status",
x=[
"hhsize",
"rooms",
svy.Cat("urbrur", ref="Urban"),
svy.Cat("electricity"),
],
family="binomial",
link="logit",
)
print(logit_model)╭────────────────────────────── GLM: Binomial (logit) ──────────────────────────────╮ │ Modeling: pov_status │ │ │ │ Observations 825 AIC 481.5853 │ │ DF Residuals 22 BIC - │ │ Deviance 467.2386 Scale 0.7535 │ │ R-squared 0.44387 R-sq (adj) 0.44116 │ │ Iterations 6 │ │ F-stat (adj) 58.51423 Prob (F-adj) <0.001 │ │ │ │ │ │ Term Coef. Std.Err. t P>|t| [0.025 0.975] │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ _intercept_ -1.57232 0.56277 -2.79387 0.0106 -2.73944 -0.40520 │ │ hhsize 0.72408 0.10260 7.05720 <0.001 0.51130 0.93687 │ │ rooms -0.67166 0.09970 -6.73660 <0.001 -0.87843 -0.46489 │ │ urbrur_Rural 1.39844 0.37209 3.75833 0.0011 0.62677 2.17010 │ │ electricity_Yes -2.36184 0.31177 -7.57552 <0.001 -3.00842 -1.71527 │ │ │ ╰───────────────────────────────────────────────────────────────────────────────────╯
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: (5, 8)
| term | estimate | std_err | conf_low | conf_high | statistic | p_value | df |
|---|---|---|---|---|---|---|---|
| str | f64 | f64 | f64 | f64 | f64 | f64 | i64 |
| "_intercept_" | -1.57232 | 0.562775 | -2.739444 | -0.405196 | -2.793869 | 0.010582 | 22 |
| "hhsize" | 0.724083 | 0.102602 | 0.5113 | 0.936867 | 7.057199 | 4.4316e-7 | 22 |
| "rooms" | -0.671661 | 0.099703 | -0.878433 | -0.464889 | -6.736601 | 9.0345e-7 | 22 |
| "urbrur_Rural" | 1.398437 | 0.37209 | 0.626769 | 2.170104 | 3.758327 | 0.001085 | 22 |
| "electricity_Yes" | -2.361843 | 0.311773 | -3.008421 | -1.715265 | -7.575518 | 1.4417e-7 | 22 |
When a Predictor Determines the Outcome (Separation)
The wealth quintile is built from the same consumption that defines pov_status, so it predicts poverty almost perfectly: no household in the average, rich or richest quintiles is poor, and every household in the poorest quintile is. This is quasi-complete separation. The estimates for those levels have no finite value: the fit pushes them toward \(\pm\infty\) and stops once the deviance settles, so the numbers it reaches depend on the convergence tolerance, not on the data.
svy detects this. The coefficients that are not identified keep the value where the iterations stopped, so predict() still works, but their standard errors, tests and intervals are NaN, as are the model Wald test and AIC. A GLM_SEPARATION warning names them, the combination that is still identified, and the levels on which the outcome never varies:
sep_model = hld_sample.glm.fit(
y="pov_status",
x=["hhsize", "rooms", svy.Cat("urbrur", ref="Urban"), svy.Cat("quint_nat")],
family="binomial",
)
print(sep_model)/var/folders/ld/zzh65kn12yn3sykchz7mcbw40000gn/T/ipykernel_11581/1207918023.py:1: SvyUserWarning: [GLM_SEPARATION] Outcome separated: some coefficients are not identified: 'pov_status' is predicted perfectly on 702 of 825 rows (quasi-complete separation), so the estimates of '_intercept_', 'quint_nat_poor', 'quint_nat_poorest', 'quint_nat_rich', 'quint_nat_richest' have no finite value: they drift toward +/-infinity and the values shown are where the iterations stopped. Their standard errors, tests and intervals are NaN, as are the AIC and the model Wald test; the other coefficients are identified. The combination _intercept_ + quint_nat_poor is identified: -0.707372.
sep_model = hld_sample.glm.fit(╭────────────────────────────── GLM: Binomial (logit) ───────────────────────────────╮ │ Modeling: pov_status │ │ │ │ Observations 825 AIC nan │ │ DF Residuals 19 BIC - │ │ Deviance 120.7610 Scale 0.1366 │ │ R-squared 0.85627 R-sq (adj) 0.85503 │ │ Iterations 20 │ │ F-stat (adj) nan Prob (F-adj) nan │ │ │ │ │ │ Term Coef. Std.Err. t P>|t| [0.025 0.975] │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ _intercept_ -23.30884 nan nan nan nan nan │ │ 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.60147 nan nan nan nan nan │ │ quint_nat_poorest 42.78578 nan nan nan nan nan │ │ quint_nat_rich 0.35737 nan nan nan nan nan │ │ quint_nat_richest 0.76412 nan nan nan nan nan │ │ │ ╰────────────────────────────────────────────────────────────────────────────────────╯
The other coefficients (hhsize, rooms, urbrur) are identified and their results stand. R’s svyglm prints finite but meaningless standard errors for the separated terms without a warning; Stata’s svy: logit drops them and reports the same estimates for the rest. The fix is in the model: collapse the levels named in the warning, or leave out a predictor that encodes the outcome.
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 │ │ where: electricity == "No" │ │ │ │ Observations 115 AIC 97.0879 │ │ DF Residuals 8 BIC - │ │ Deviance 90.4262 Scale 1.4452 │ │ 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 │ │ │ ╰────────────────────────────────────────────────────────────────────────────────╯
Fitting with Replicate Weights
glm.fit() uses Taylor linearization unless you pass method="replication", as in estimation and the categorical tests: replicate weights on the design are not used implicitly. With method="replication" the model is refitted with each replicate weight, and the degrees of freedom are the replicate df less the number of non-intercept coefficients:
jk_sample = hld_sample.weighting.create_jk_wgts()
logit_model_jk = jk_sample.glm.fit(
y="pov_status",
x=["hhsize", "rooms", svy.Cat("urbrur", ref="Urban"), svy.Cat("electricity")],
family="binomial",
method="replication",
)
print(logit_model_jk)╭────────────────────────────── GLM: Binomial (logit) ──────────────────────────────╮ │ Modeling: pov_status │ │ │ │ Observations 825 AIC 484.9368 │ │ DF Residuals 22 BIC - │ │ Deviance 467.2386 Scale 0.7535 │ │ R-squared 0.44387 R-sq (adj) 0.44116 │ │ Iterations 6 │ │ F-stat (adj) 53.76662 Prob (F-adj) <0.001 │ │ │ │ │ │ Term Coef. Std.Err. t P>|t| [0.025 0.975] │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ _intercept_ -1.57232 0.64564 -2.43527 0.0234 -2.91131 -0.23333 │ │ hhsize 0.72408 0.11628 6.22684 <0.001 0.48292 0.96524 │ │ rooms -0.67166 0.11266 -5.96180 <0.001 -0.90531 -0.43802 │ │ urbrur_Rural 1.39844 0.41071 3.40490 0.0025 0.54667 2.25020 │ │ electricity_Yes -2.36184 0.32047 -7.36989 <0.001 -3.02646 -1.69723 │ │ │ ╰───────────────────────────────────────────────────────────────────────────────────╯
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 22.0 │ │ Mean ŷ 0.2364 Mean SE 0.0339 │ │ Min ŷ 0.0000 Max ŷ 0.9983 │ │ Mean resid -0.0000 Std resid 0.3160 │ │ │ │ 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.010976 | 0.003874 | 0.005267 | 0.022733 | -0.010976 |
| 0.001477 | 0.000734 | 0.000527 | 0.004135 | -0.001477 |
| 0.01156 | 0.003621 | 0.006025 | 0.022069 | -0.01156 |
| 0.005638 | 0.002153 | 0.00255 | 0.012418 | -0.005638 |
| 0.002741 | 0.001269 | 0.001048 | 0.007148 | -0.002741 |
| 0.005939 | 0.002133 | 0.002817 | 0.01248 | -0.005939 |
| 0.005638 | 0.002153 | 0.00255 | 0.012418 | -0.005638 |
| 0.022381 | 0.006575 | 0.012128 | 0.040944 | -0.022381 |
| 0.005638 | 0.002153 | 0.00255 | 0.012418 | -0.005638 |
| 0.005638 | 0.002153 | 0.00255 | 0.012418 | -0.005638 |
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.0933
Mean predicted probability (actual): 0.2364Marginal 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.051946 0.016124 [0.018507, 0.085385] │ │ 3.00 0.128538 0.017978 [0.091255, 0.165822] │ │ 5.00 0.264870 0.025731 [0.211507, 0.318232] │ │ 7.00 0.473313 0.056696 [0.355733, 0.590892] │ │ 10.00 0.807549 0.071214 [0.659860, 0.955239] │ │ │ ╰──────────────────────────────────────────────────────╯
pred_margins.to_polars()
shape: (5, 6)
| term | margin | se | lci | uci | value |
|---|---|---|---|---|---|
| str | f64 | f64 | f64 | f64 | i64 |
| "hhsize" | 0.051946 | 0.016124 | 0.018507 | 0.085385 | 1 |
| "hhsize" | 0.128538 | 0.017978 | 0.091255 | 0.165822 | 3 |
| "hhsize" | 0.26487 | 0.025731 | 0.211507 | 0.318232 | 5 |
| "hhsize" | 0.473313 | 0.056696 | 0.355733 | 0.590892 | 7 |
| "hhsize" | 0.807549 | 0.071214 | 0.65986 | 0.955239 | 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.062172 0.008076 [0.045424, 0.078920] │ │ │ ╰──────────────────────────────────────────────╯
╭─────── GLM Margins: rooms (ame, 95% CI) ────────╮ │ │ │ Margin SE 95% CI │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ -0.057671 0.010035 [-0.078483, -0.036859] │ │ │ ╰─────────────────────────────────────────────────╯
╭───────────── GLM Margins: urbrur (ame, 95% CI) ──────────────╮ │ │ │ Contrast Margin SE 95% CI │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ Rural - Urban 0.134749 0.037592 [0.056788, 0.212710] │ │ │ ╰──────────────────────────────────────────────────────────────╯
╭────────── GLM Margins: electricity (ame, 95% CI) ──────────╮ │ │ │ Contrast Margin SE 95% CI │ │ ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ │ │ Yes - No -0.292671 0.055480 [-0.407729, -0.177613] │ │ │ ╰────────────────────────────────────────────────────────────╯
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 |
"inverse_gaussian" |
Positive continuous data with a heavier right tail than Gamma |
"negative_binomial" |
Counts too variable for Poisson (overdispersed) |
Link Function (link)
Connects the linear predictors (Xβ) to the expected mean of the distribution (μ):
| Link | Description | Canonical for | Also admitted by |
|---|---|---|---|
"identity" |
No transformation (μ = Xβ) | Gaussian | every family |
"logit" |
Log-odds: ln(p/(1−p)) = Xβ | Binomial | — |
"probit" |
Inverse normal CDF: Φ⁻¹(p) = Xβ | — | Binomial |
"cauchit" |
Inverse Cauchy CDF: tan(π(p−½)) = Xβ | — | Binomial |
"cloglog" |
Complementary log-log: ln(−ln(1−p)) = Xβ | — | Binomial |
"log" |
Log link: ln(μ) = Xβ | Poisson, Negative Binomial | all but Binomial-only use |
"sqrt" |
Square root: √μ = Xβ | — | Poisson, Negative Binomial |
"inverse" |
Reciprocal: 1/μ = Xβ | Gamma | Gaussian, Inverse Gaussian |
"inverse_squared" |
1/μ² = Xβ | Inverse Gaussian | — |
When link is omitted, the canonical link for the family is used automatically. Each family admits exactly the links R’s family constructors do; an unusable pairing (say binomial with inverse_squared) raises rather than fitting a model that means nothing.
probit and cauchit are the two alternatives to logit for a binary outcome: probit has lighter tails, cauchit much heavier ones, so cauchit is the choice when a few observations sit far out on the linear predictor and should not dominate the fit. sqrt is the variance-stabilising link for counts.
Overdispersed counts
Poisson pins the variance to the mean. When counts are more variable than that — which is most real count data — use the negative binomial, whose variance is μ + μ²/θ:
model = sample.glm.fit(y="visits", x=["age", svy.Cat("region")], family="nb")
model.fitted.stats.theta, model.fitted.stats.theta_seθ is estimated alongside the coefficients, and the standard errors come from the joint (coefficients, θ) design-based sandwich, so they account for having estimated it. stats.theta_se is the design-based standard error of θ itself. Pass theta= instead and it is treated as known: the variance then conditions on it, and stats.theta_se is None.
A quick check on whether you need it: fit Poisson first and look at stats.scale, the Pearson dispersion. Far above 1 means the Poisson variance assumption is wrong.
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.