Statsmodels

statsmodels Skill

statsmodels general-purpose statistical modeling library for Python. Covers OLS/WLS/GLS, GLM (logit, probit, Poisson, negative binomial), discrete choice models, time series (ARIMA, SARIMAX, VAR), mixed effects (MixedLM), robust regression, hypothesis tests, and comprehensive diagnostics. Supports R-style formula API. Use when fitting regressions without fixed effects, running GLMs or logit/probit, analyzing time series, or using formula syntax. For fixed effects or DiD, use pyfixest; for panel/IV/system models, use linearmodels.

Comprehensive skill for statistical modeling with statsmodels. Use decision trees below to find the right guidance, then load detailed references.

What is statsmodels?

statsmodels is the general-purpose statistical modeling library for Python:

• Two APIs: Formula API (smf.ols("y ~ x1 + x2", data=df)) for R-style modeling, and array API (sm.OLS(y, X)) for programmatic control
• Broad model coverage: OLS, WLS, GLS, GLM (all families), logit, probit, multinomial, count models, zero-inflated models, quantile regression, robust regression
• Time series: ARIMA, SARIMAX, VAR, exponential smoothing, state space models, unit root tests
• Diagnostics: Heteroskedasticity tests, normality tests, specification tests, VIF, influence measures, residual analysis
• Hypothesis testing: t-tests, F-tests, Wald tests, likelihood ratio tests, multiple comparison corrections

How to Use This Skill

Reference File Structure

Reading Order

1. New to statsmodels? Start with quickstart.md then linear-models.md 2. Need GLM or logit/probit? Read quickstart.md then glm-discrete.md 3. Time series analysis? Read quickstart.md then time-series.md 4. Checking model assumptions? Read diagnostics.md 5. Coming from R? Read quickstart.md (formula API mirrors R syntax)

Related Skills

• pyfixest: Use instead of statsmodels when your model needs absorbed fixed effects, IV with FE, or difference-in-differences. pyfixest is faster for FE models; statsmodels is broader for everything else
• linearmodels: Use for panel data models (FE, RE, between, first difference, Fama-MacBeth), IV/GMM without FE (2SLS, LIML, GMM), system estimation (SUR, 3SLS), and asset pricing. Built on top of statsmodels; extends it for structured data
• svy: Use for survey-weighted regression and estimation with complex survey designs. Important: statsmodels WLS is NOT equivalent to survey-weighted regression — WLS handles heteroscedastic errors but does not account for stratification, clustering, or finite population corrections. If your data comes from a complex probability survey (NHANES, ACS PUMS, CPS, ECLS-K, etc.), load the svy skill instead
• data-scientist: Provides methodology guidance (when to use which model, assumption checking protocol, interpretation). Load alongside statsmodels for the "why"; statsmodels provides the "how"
• polars: Data manipulation before modeling. statsmodels accepts pandas DataFrames; convert with df.to_pandas() if using Polars
• plotnine: Publication-quality visualization of model results and diagnostics

Quick Decision Trees

"I need to fit a regression model"

What kind of regression?
├─ Linear (continuous outcome)
│   ├─ Basic OLS → ./references/linear-models.md
│   ├─ Weighted least squares → ./references/linear-models.md
│   │   (⚠ WLS ≠ survey-weighted regression — for complex surveys, use `svy` skill)
│   ├─ Correlated errors (GLS) → ./references/linear-models.md
│   ├─ Robust to outliers (M-estimator) → ./references/linear-models.md
│   └─ Quantile regression → ./references/linear-models.md
├─ Binary outcome (0/1)
│   ├─ Logit → ./references/glm-discrete.md
│   └─ Probit → ./references/glm-discrete.md
├─ Count outcome (0, 1, 2, ...)
│   ├─ Poisson → ./references/glm-discrete.md
│   ├─ Negative binomial → ./references/glm-discrete.md
│   └─ Zero-inflated → ./references/glm-discrete.md
├─ Multinomial (3+ categories)
│   └─ Multinomial logit → ./references/glm-discrete.md
├─ GLM (custom family/link)
│   └─ GLM framework → ./references/glm-discrete.md
└─ Need fixed effects?
    └─ Use pyfixest instead (faster FE absorption)

"I need to analyze time series"

What time series task?
├─ Forecast a single series
│   ├─ ARIMA / SARIMAX → ./references/time-series.md
│   └─ Exponential smoothing → ./references/time-series.md
├─ Multiple interrelated series
│   └─ VAR / VECM → ./references/time-series.md
├─ Test for stationarity
│   ├─ ADF test → ./references/time-series.md
│   └─ KPSS test → ./references/time-series.md
├─ Examine autocorrelation
│   └─ ACF / PACF → ./references/time-series.md
└─ Structural time series
    └─ Unobserved components → ./references/time-series.md

"I need to check model assumptions"

What assumption to check?
├─ Heteroskedasticity → ./references/diagnostics.md
│   ├─ Breusch-Pagan test
│   └─ White test
├─ Normality of residuals → ./references/diagnostics.md
│   ├─ Jarque-Bera test
│   └─ Shapiro-Wilk test
├─ Specification / functional form → ./references/diagnostics.md
│   └─ RESET test
├─ Multicollinearity → ./references/diagnostics.md
│   ├─ VIF
│   └─ Condition number
├─ Influential observations → ./references/diagnostics.md
│   ├─ Cook's distance
│   └─ Leverage / DFFITS
├─ Serial correlation → ./references/diagnostics.md
│   └─ Durbin-Watson / Breusch-Godfrey
└─ All of the above → ./references/diagnostics.md

"I need to test hypotheses"

What kind of test?
├─ Single coefficient significance → ./references/hypothesis-testing.md
├─ Joint significance (F-test) → ./references/hypothesis-testing.md
├─ Linear restrictions (Wald) → ./references/hypothesis-testing.md
├─ Compare nested models (LR test) → ./references/hypothesis-testing.md
├─ Multiple comparisons correction → ./references/hypothesis-testing.md
└─ Chi-squared test → ./references/hypothesis-testing.md

"Something isn't working"

Common issues?
├─ Missing constant / intercept → ./references/gotchas.md
├─ Convergence warnings → ./references/gotchas.md
├─ predict() errors → ./references/gotchas.md
├─ Formula parsing issues → ./references/gotchas.md
├─ summary() formatting → ./references/gotchas.md
├─ statsmodels vs pyfixest → ./references/gotchas.md
└─ General troubleshooting → ./references/gotchas.md

File-First Execution in Research Workflows

Important: In data research pipelines (see CLAUDE.md), statsmodels analyses are executed through script files, not interactively. This ensures auditability and reproducibility.

The pattern: 1. Write model code to scripts/stage8_analysis/{step}_{model-name}.py 2. Execute via Bash with automatic output capture wrapper script 3. Validation results get automatically embedded in scripts as comments 4. If failed, create versioned copy for fixes

Closely read agent_reference/SCRIPT_EXECUTION_REFERENCE.md for the mandatory file-first execution protocol covering complete code file writing, output capture, and file versioning rules.

See:

• agent_reference/SCRIPT_EXECUTION_REFERENCE.md — Script execution protocol and format with validation

The examples below show statsmodels syntax. In research workflows, wrap them in scripts following the file-first pattern.

---

Quick Reference

Essential Imports

import statsmodels.api as sm           # Array API
import statsmodels.formula.api as smf  # Formula API (R-style)

Core Operations

Formula Syntax

# Additive terms
"y ~ x1 + x2 + x3"

# Interaction (with main effects)
"y ~ x1 * x2"           # equivalent to x1 + x2 + x1:x2

# Interaction only (no main effects)
"y ~ x1 : x2"

# Categorical variable
"y ~ C(region)"          # treatment coding (default)
"y ~ C(region, Treatment(reference='West'))"  # explicit reference

# Suppress intercept
"y ~ x1 + x2 - 1"

# Polynomial
"y ~ x1 + I(x1**2)"     # I() protects Python operators

Topic Index

Citation

When this library is used as a primary analytical tool, include in the report's Software & Tools references:

Seabold, S. & Perktold, J. (2010). "Statsmodels: Econometric and Statistical Modeling with Python." Proceedings of the 9th Python in Science Conference.

Cite when: statsmodels is used for GLM estimation, time series modeling, or statistical hypothesis testing central to the analysis. Do not cite when: Only used for post-estimation diagnostics supporting another library's primary estimation.

For method-specific citations (e.g., individual estimators or techniques), consult the reference files in this skill and agent_reference/CITATION_REFERENCE.md.

statsmodels Diagnostics Reference

Regression diagnostics for OLS and GLM models using statsmodels v0.14.6. All examples assume a fitted OLS results object unless otherwise noted.

---

Contents

1. Heteroskedasticity Tests 2. Normality Tests 3. Specification Tests 4. Multicollinearity 5. Influence Measures 6. Residual Analysis 7. Serial Correlation 8. Comprehensive Diagnostic Checklist 9. References and Further Reading

---

Heteroskedasticity Tests

Breusch-Pagan Test

from statsmodels.stats.diagnostic import het_breuschpagan

# Fit OLS first
results = smf.ols("y ~ x1 + x2", data=df).fit()

# Run test (uses residuals and regressors)
lm_stat, lm_pvalue, f_stat, f_pvalue = het_breuschpagan(
    results.resid, results.model.exog
)

• H0: homoskedasticity (constant variance)
• Reject if p-value < 0.05 → heteroskedasticity present
• Returns four values: LM statistic, LM p-value, F-statistic, F p-value
• The LM test assumes normal errors; the F-statistic is more robust to

non-normality. Use the F variant (3rd and 4th return values) as the default

White Test

from statsmodels.stats.diagnostic import het_white

lm_stat, lm_pvalue, f_stat, f_pvalue = het_white(
    results.resid, results.model.exog
)

• More general than Breusch-Pagan (doesn't assume specific functional form)
• Tests for heteroskedasticity AND specification error simultaneously
• Adds squared terms and cross-products of regressors to the auxiliary regression
• Can lose power with many regressors because the auxiliary regression grows large

What to Do If Heteroskedasticity Is Detected

• Use robust standard errors: .fit(cov_type='HC1') or 'HC3' (HC3 preferred

in small samples)

• Use WLS (sm.WLS) if the variance structure is known or estimable
• Consider transforming the dependent variable (log, sqrt) if the relationship

is multiplicative

---

Normality Tests

Jarque-Bera

Tests skewness and kurtosis jointly:

from statsmodels.stats.stattools import jarque_bera

jb_stat, jb_pvalue, skew, kurtosis = jarque_bera(results.resid)

• H0: residuals are normally distributed
• Also available directly in results.summary() diagnostics section
• Returns four values: JB statistic, p-value, skewness, kurtosis

Omnibus Test (D'Agostino-Pearson)

from statsmodels.stats.stattools import omni_normtest

stat, pvalue = omni_normtest(results.resid)

• Also reported in results.summary() output
• Combines skewness and kurtosis tests using chi-squared distribution

Shapiro-Wilk (scipy)

from scipy.stats import shapiro

stat, pvalue = shapiro(results.resid)

• More powerful than Jarque-Bera for small samples
• For N > 5000, scipy issues a warning that p-values may not be accurate
• H0: sample is drawn from a normal distribution

Anderson-Darling (scipy)

from scipy.stats import anderson

result = anderson(results.resid, dist='norm')
# result.statistic — the test statistic
# result.critical_values — array of critical values
# result.significance_level — corresponding significance levels [15, 10, 5, 2.5, 1]

• Returns critical values at multiple significance levels rather than a single

p-value

• Reject H0 when result.statistic > result.critical_values[idx] for the

chosen significance level

Practical Note on Large Samples

In large samples, normality tests almost always reject because any minor departure from normality becomes statistically detectable. Focus on the degree of non-normality — examine QQ plots, skewness, and kurtosis values numerically — rather than relying solely on p-values. OLS inference is robust to moderate non-normality when N is large.

---

Specification Tests

RESET Test (Ramsey)

Ramsey's Regression Equation Specification Error Test:

from statsmodels.stats.diagnostic import linear_reset

reset_result = linear_reset(results, power=3, use_f=True)
print(reset_result.fvalue)   # F-statistic
print(reset_result.pvalue)   # p-value
print(reset_result.df_denom) # denominator df

• H0: model is correctly specified
• power=3 adds ŷ² and ŷ³ as auxiliary regressors (default behavior)
• use_f=True returns F-statistic; use_f=False returns chi-squared
• Rejection suggests omitted variables, wrong functional form, or omitted

interactions — but does not identify which

Harvey-Collier Test (Linearity)

from statsmodels.stats.diagnostic import linear_harvey_collier

t_stat, pvalue = linear_harvey_collier(results)

• H0: relationship is linear
• Tests using recursive residuals — sensitive to structural breaks as well as

nonlinearity

---

Multicollinearity

Variance Inflation Factor (VIF)

from statsmodels.stats.outliers_influence import variance_inflation_factor

# X must include a constant column
X = sm.add_constant(df[["x1", "x2", "x3"]])

vif_data = pd.DataFrame({
    "Variable": X.columns,
    "VIF": [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
})
print(vif_data)

Interpretation thresholds:

• VIF = 1: no collinearity
• VIF > 5: moderate concern
• VIF > 10: serious multicollinearity

The constant's VIF is meaningless — exclude it from display or interpretation. variance_inflation_factor takes the full design matrix (as a numpy array) and a column index; loop over column indices to compute VIF for each variable.

Condition Number

# Available directly from results object
print(results.condition_number)

# Or compute manually
import numpy as np
cond = np.linalg.cond(results.model.exog)

• Condition number > 30: potential multicollinearity concern
• Also reported in the bottom section of results.summary() output

---

Influence Measures

OLSInfluence (Comprehensive)

from statsmodels.stats.outliers_influence import OLSInfluence

influence = OLSInfluence(results)

# Cook's distance
cooks_d, pvalues = influence.cooks_distance
# Rule of thumb: Cook's D > 4/n is influential

# Leverage (hat matrix diagonal)
leverage = influence.hat_matrix_diag
# Rule of thumb: leverage > 2*(k+1)/n is high leverage

# DFFITS
dffits = influence.dffits[0]
# Rule of thumb: |DFFITS| > 2*sqrt((k+1)/n)

# DFBETAS (influence on each coefficient separately)
dfbetas = influence.dfbetas
# Shape: (n_obs, n_params)
# Rule of thumb: |DFBETAS| > 2/sqrt(n)

# Externally studentized residuals
student_resid = influence.resid_studentized_external
# |studentized resid| > 2 is a potential outlier

Where k = number of predictors (excluding constant) and n = number of observations.

Summary Frame (All Metrics at Once)

summary = influence.summary_frame()
print(summary.columns.tolist())
# Includes: dfb_Intercept, dfb_x1, ..., cooks_d, standard_resid,
#           hat_diag, dffits_internal, student_resid

The summary frame is the most efficient way to flag observations that exceed multiple thresholds simultaneously. Filter on cooks_d and hat_diag in combination to distinguish high-influence from high-leverage points.

---

Residual Analysis

Residual Types

results.resid                              # Raw residuals (y - ŷ)
results.resid_pearson                      # Pearson residuals (for GLM)
influence.resid_studentized_external       # Externally studentized
influence.resid_studentized_internal       # Internally studentized

Externally studentized residuals delete observation i before estimating variance — they follow a t-distribution under H0 and are preferred for outlier detection.

Partial Regression Plots (Added Variable Plots)

from statsmodels.graphics.regressionplots import plot_partregress_grid

fig = plot_partregress_grid(results)
fig.tight_layout()
fig.savefig("partial_regression.png", dpi=300, bbox_inches='tight')
plt.close()

Each panel shows the marginal relationship between one predictor and the outcome after removing the linear effect of all other predictors. Useful for spotting nonlinearity and influential points for individual coefficients.

Residual-vs-Fitted Plot

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(8, 6))
ax.scatter(results.fittedvalues, results.resid, alpha=0.5)
ax.axhline(y=0, color='r', linestyle='--')
ax.set_xlabel("Fitted Values")
ax.set_ylabel("Residuals")
ax.set_title("Residuals vs Fitted")
fig.savefig("resid_vs_fitted.png", dpi=300, bbox_inches='tight')
plt.close()

Look for: random scatter around zero (good), funnel pattern (heteroskedasticity), curved pattern (nonlinearity).

QQ Plot

from statsmodels.graphics.gofplots import qqplot

fig = qqplot(results.resid, line='45', fit=True)
fig.savefig("qq_plot.png", dpi=300, bbox_inches='tight')
plt.close()

line='45' draws the 45-degree reference line. fit=True fits location and scale parameters before plotting. Heavy tails appear as S-shaped departure from the line; skewness appears as systematic curvature.

---

Serial Correlation

Durbin-Watson Statistic

from statsmodels.stats.stattools import durbin_watson

dw = durbin_watson(results.resid)

• DW ≈ 2: no autocorrelation
• DW < 2: positive autocorrelation (values closer to 0 = stronger)
• DW > 2: negative autocorrelation (values closer to 4 = stronger)
• Also reported in results.summary() output
• Cannot be used when the model contains a lagged dependent variable

Breusch-Godfrey Test

More general than Durbin-Watson; works with lagged dependent variables:

from statsmodels.stats.diagnostic import acorr_breusch_godfrey

lm_stat, lm_pvalue, f_stat, f_pvalue = acorr_breusch_godfrey(
    results, nlags=4
)

• H0: no serial correlation up to lag nlags
• nlags should be set based on data frequency (e.g., 4 for quarterly, 12

for monthly)

• Returns LM statistic, LM p-value, F-statistic, F p-value

Ljung-Box Test

from statsmodels.stats.diagnostic import acorr_ljungbox

lb_result = acorr_ljungbox(results.resid, lags=10, return_df=True)
print(lb_result)
# DataFrame with columns: lb_stat, lb_pvalue
# One row per lag

• H0: residuals are independently distributed up to the specified lag
• return_df=True returns results as a DataFrame indexed by lag number
• Primarily used for time series residuals from ARIMA models, but applicable

to any regression residuals when serial correlation is a concern

---

Comprehensive Diagnostic Checklist

---

References and Further Reading

• statsmodels diagnostic tests: https://www.statsmodels.org/stable/diagnostic.html
• statsmodels influence measures: https://www.statsmodels.org/stable/generated/statsmodels.stats.outliers_influence.OLSInfluence.html
• statsmodels graphics: https://www.statsmodels.org/stable/graphics.html
• Greene, W.H. (2018). Econometric Analysis, 8th ed. — Chapters 4-5

(diagnostics)

• Belsley, Kuh, & Welsch (1980). Regression Diagnostics. Wiley.

statsmodels GLM and Discrete Models Reference

statsmodels v0.14.6 — syntax and library guidance only.

---

Contents

1. GLM Framework 2. Logit and Probit (Discrete Choice) 3. Multinomial Logit (MNLogit) 4. Poisson and Negative Binomial Count Models 5. Zero-Inflated Models 6. Separation Detection 7. Marginal Effects 8. Exposure and Offset for Rate Models 9. Ordered logit / probit 10. References and Further Reading

---

GLM Framework

The GLM class provides a unified interface for all generalized linear models. Specify the distribution of y via a family object and optionally override the default link function.

Basic Usage

import statsmodels.api as sm
import statsmodels.formula.api as smf

# Formula API
results = smf.glm("y ~ x1 + x2", data=df,
                   family=sm.families.Binomial()).fit()

# Array API
X = sm.add_constant(df[["x1", "x2"]])
results = sm.GLM(df["y"], X, family=sm.families.Binomial()).fit()

print(results.summary())

Family Classes

All families are under sm.families:

NegativeBinomial(alpha=1.0): alpha is the overdispersion parameter (also called the heterogeneity parameter). Larger alpha = more overdispersion. When alpha → 0, approaches Poisson.

Tweedie(var_power=p): variance = μ^p. p=1 is Poisson, p=2 is Gamma. Values between 1 and 2 produce a compound Poisson-Gamma (Tweedie) distribution with a point mass at zero.

Link Functions

All links are under sm.families.links:

Specifying a Non-Default Link

# Probit link with Binomial family (probit regression via GLM)
results = smf.glm("y ~ x1 + x2", data=df,
                   family=sm.families.Binomial(
                       link=sm.families.links.Probit()
                   )).fit()

# Log link with Gaussian family (log-linear model for continuous y)
results = smf.glm("y ~ x1 + x2", data=df,
                   family=sm.families.Gaussian(
                       link=sm.families.links.Log()
                   )).fit()

# Complementary log-log link (Binomial family — proportional hazards interpretation)
results = smf.glm("y ~ x1 + x2", data=df,
                   family=sm.families.Binomial(
                       link=sm.families.links.CLogLog()
                   )).fit()

GLM Result Attributes

results.params          # Coefficients on the link scale
results.pvalues         # P-values
results.conf_int()      # Confidence intervals (default alpha=0.05)
results.fittedvalues    # Fitted values on the response scale (μ̂)
results.resid_response  # Response residuals: y - μ̂
results.resid_pearson   # Pearson residuals: (y - μ̂) / sqrt(V(μ̂))
results.resid_deviance  # Deviance residuals (signed sqrt of deviance contributions)
results.deviance        # Residual deviance
results.null_deviance   # Null model deviance (intercept only)
results.df_resid        # Residual degrees of freedom
results.df_model        # Model degrees of freedom (excl. constant)
results.aic             # AIC = deviance + 2 * (df_model + 1)
results.bic             # BIC
results.llf             # Log-likelihood
results.nobs            # Number of observations

Dispersion parameter (phi):

results.scale           # Estimated dispersion (phi); for Poisson/Binomial fixed at 1

Prediction from GLM

# Predict on new data — returns values on the response scale (μ̂)
y_pred = results.predict(new_data)

# Predict with confidence intervals (linear predictor scale by default)
pred = results.get_prediction(new_data)
frame = pred.summary_frame(alpha=0.05)
# Columns: mean, mean_se, mean_ci_lower, mean_ci_upper

---

Logit and Probit (Discrete Choice)

Dedicated classes for binary dependent variables. These offer additional features beyond GLM(family=Binomial()): McFadden's R², pred_table(), and direct get_margeff().

Logit

# Formula API
results = smf.logit("y ~ x1 + x2 + x3", data=df).fit()

# Array API
results = sm.Logit(df["y"], sm.add_constant(df[["x1", "x2", "x3"]])).fit()

print(results.summary())

Logit coefficients are log-odds ratios. Exponentiate for odds ratios:

import numpy as np

or_df = pd.DataFrame({
    "OR": np.exp(results.params),
    "ci_lo": np.exp(results.conf_int()[0]),
    "ci_hi": np.exp(results.conf_int()[1])
})
print(or_df)

Probit

# Formula API
results = smf.probit("y ~ x1 + x2 + x3", data=df).fit()

# Array API
results = sm.Probit(df["y"], sm.add_constant(df[["x1", "x2", "x3"]])).fit()

Probit coefficients scale with the standard normal CDF. They are not directly interpretable as odds ratios; use marginal effects (see section below).

Key Result Attributes for Binary Models

results.llf             # Log-likelihood of the fitted model
results.llnull          # Log-likelihood of the null model (intercept only)
results.prsquared       # McFadden's pseudo R² = 1 - llf / llnull
                        # Values of 0.2-0.4 indicate good fit (lower than OLS R²)
results.pred_table()    # Confusion matrix as numpy array [[TN, FP], [FN, TP]]
                        # Default threshold: 0.5

Confusion matrix with labels:

table = results.pred_table()
print(f"True Negatives:  {table[0,0]}")
print(f"False Positives: {table[0,1]}")
print(f"False Negatives: {table[1,0]}")
print(f"True Positives:  {table[1,1]}")
accuracy = (table[0,0] + table[1,1]) / table.sum()
print(f"Accuracy: {accuracy:.3f}")

Predicted Probabilities

# In-sample fitted probabilities
prob_hat = results.predict()          # or results.fittedvalues

# Out-of-sample
prob_new = results.predict(new_data)

# Using get_prediction for confidence intervals
pred = results.get_prediction(new_data)
frame = pred.summary_frame()
# Columns: mean (probability), mean_se, mean_ci_lower, mean_ci_upper

---

Multinomial Logit (MNLogit)

For outcomes with 3 or more unordered categories. The base (reference) category is J=0.

Usage

# y must be integer-coded: 0, 1, 2, ...
# Category 0 is the reference

# Array API
results = sm.MNLogit(df["y"], sm.add_constant(df[["x1", "x2"]])).fit()

# Formula API
results = smf.mnlogit("y ~ x1 + x2", data=df).fit()

print(results.summary())

Coefficients and Shape

results.params          # Shape: (n_params, J-1)
                        # Columns are categories 1, 2, ..., J-1 (relative to category 0)
results.pvalues         # Same shape as params
results.conf_int()      # Shape: (n_params * (J-1), 2) — stacked by category

Extracting coefficients for a specific category:

# Coefficients for category 1 vs. reference (category 0)
coef_cat1 = results.params.iloc[:, 0]

# Coefficients for category 2 vs. reference
coef_cat2 = results.params.iloc[:, 1]

Predicted Probabilities

# Returns (n_obs, J) array — probability of each category for each observation
probs = results.predict()                    # in-sample
probs_new = results.predict(new_data_X)     # out-of-sample (array API: pass design matrix)

# For formula API new data:
probs_new = results.predict(new_data_df)

Changing the Reference Category

To use a different reference category, recode y before fitting:

# Make category 2 the reference by shifting values
df["y_recoded"] = df["y"].map({0: 1, 1: 2, 2: 0})  # category 2 → 0
results = smf.mnlogit("y_recoded ~ x1 + x2", data=df).fit()

---

Poisson and Negative Binomial Count Models

Poisson

Assumes equidispersion: E[y] = Var[y] = μ. Use when counts are genuinely Poisson-distributed.

# Formula API
results = smf.poisson("count ~ x1 + x2", data=df).fit()

# Array API
results = sm.Poisson(df["count"], sm.add_constant(df[["x1", "x2"]])).fit()

print(results.summary())

Checking for overdispersion (Var[y] >> E[y]):

mean_count = df["count"].mean()
var_count = df["count"].var()
dispersion_ratio = var_count / mean_count
print(f"Dispersion ratio: {dispersion_ratio:.2f}")
# >> 1 suggests overdispersion → use NegBin

Negative Binomial

Adds a random effect to the Poisson rate, producing Var[y] = μ + alpha*μ². Use when variance exceeds the mean (overdispersion).

# NegativeBinomial (NB2 parameterization — variance = μ + alpha*μ²)
results = sm.NegativeBinomial(
    df["count"], sm.add_constant(df[["x1", "x2"]])
).fit()

# Via GLM (fixes alpha instead of estimating it)
results = smf.glm("count ~ x1 + x2", data=df,
                   family=sm.families.NegativeBinomial(alpha=1.0)).fit()

# NegativeBinomialP (flexible parameterization: variance = μ + alpha*μ^p)
results = sm.NegativeBinomialP(
    df["count"], sm.add_constant(df[["x1", "x2"]])
).fit()

sm.NegativeBinomial estimates alpha from the data. sm.families.NegativeBinomial(alpha=...) fixes it.

Interpreting Count Model Coefficients

Poisson and NegBin coefficients are log rate ratios. Exponentiate for incidence rate ratios (IRR):

import numpy as np

irr = np.exp(results.params)
irr_ci = np.exp(results.conf_int())
irr_df = pd.concat([irr, irr_ci], axis=1)
irr_df.columns = ["IRR", "CI_lo", "CI_hi"]
print(irr_df)
# IRR > 1: positive association; IRR < 1: negative association

NB2 vs. NB-P

Accessing the estimated alpha (overdispersion) from NB2:

# alpha is the last element of params in NegativeBinomial
print(f"Estimated alpha: {results.params['alpha']:.4f}")

---

Zero-Inflated Models

For count data with more zeros than a Poisson or NegBin can explain. These model two latent processes: (1) whether an observation is a structural zero, and (2) the count conditional on being non-zero.

Zero-Inflated Poisson (ZIP)

# exog_infl: covariates for the zero-inflation (structural zero) equation
# Can differ from the count equation covariates
results = sm.ZeroInflatedPoisson(
    df["count"],
    sm.add_constant(df[["x1", "x2"]]),          # count model covariates
    exog_infl=sm.add_constant(df[["z1"]]),        # inflation model covariates
    inflation="logit"                             # logit model for P(structural zero)
).fit()

# If inflation model has same covariates as count model, omit exog_infl
results = sm.ZeroInflatedPoisson(
    df["count"],
    sm.add_constant(df[["x1", "x2"]]),
    inflation="logit"
).fit()

inflation options: "logit" (default) or "probit".

Zero-Inflated Negative Binomial (ZINB)

results = sm.ZeroInflatedNegativeBinomialP(
    df["count"],
    sm.add_constant(df[["x1", "x2"]]),
    exog_infl=sm.add_constant(df[["z1"]]),
    inflation="logit"
).fit()

Zero-Inflated Generalized Poisson (ZIGP)

results = sm.ZeroInflatedGeneralizedPoisson(
    df["count"],
    sm.add_constant(df[["x1", "x2"]]),
    exog_infl=sm.add_constant(df[["z1"]]),
    inflation="logit"
).fit()

Reading Zero-Inflated Results

The summary contains two blocks of coefficients:

print(results.summary())
# "inflate_*" parameters: coefficients for the zero-inflation (logit) equation
# Positive inflate coefficient → higher probability of structural zero
# Remaining parameters: count model coefficients (log scale)

Accessing coefficients by block:

all_params = results.params

# Count model coefficients (prefix varies — check results.params.index)
count_params = all_params[~all_params.index.str.startswith("inflate")]
inflate_params = all_params[all_params.index.str.startswith("inflate")]

Comparing ZIP vs. Poisson with Vuong Test

The Vuong test compares non-nested models (e.g., ZIP vs. Poisson):

from statsmodels.discrete.count_model import ZeroInflatedPoisson
from statsmodels.discrete.discrete_model import Poisson

poisson_res = sm.Poisson(df["count"], X).fit()
zip_res = sm.ZeroInflatedPoisson(df["count"], X, inflation="logit").fit()

# Vuong statistic > 1.96 favors ZIP; < -1.96 favors Poisson
# Not directly built into statsmodels — compute manually or use AIC comparison
print(f"Poisson AIC: {poisson_res.aic:.2f}")
print(f"ZIP AIC:     {zip_res.aic:.2f}")

---

Separation Detection

Perfect or quasi-perfect separation occurs when a predictor (or combination) perfectly predicts y=0 or y=1 in a binary model. The MLE does not exist; coefficients diverge to ±infinity.

Symptoms

results = smf.logit("y ~ x1 + x2 + C(state)", data=df).fit()

# Warning at convergence is the first signal:
# "Maximum Likelihood optimization failed to converge"
print(results.mle_retvals["converged"])   # False if not converged

# Very large coefficient magnitudes
print(results.params[results.params.abs() > 10])

# Very large standard errors
print(results.bse[results.bse > 100])

Diagnosing the Source

# For a binary predictor suspected of causing separation:
print(pd.crosstab(df["y"], df["x_binary"]))
# If any cell is 0, that predictor perfectly separates y

# For a continuous predictor:
print(df.groupby("y")["x_continuous"].agg(["min", "max"]))
# If ranges do not overlap, separation is present

# Check all categorical predictors
for col in df.select_dtypes("object").columns:
    ct = pd.crosstab(df["y"], df[col])
    zero_cells = (ct == 0).any().any()
    if zero_cells:
        print(f"Potential separation: {col}")
        print(ct)

Fixes

statsmodels does not natively implement Firth (penalized) logit. Options:

1. Reduce model complexity: remove the separating predictor, combine sparse categories 2. Aggregate rare categories: merge categories with few observations in either y=0 or y=1 3. External package: use the firthlogist package for Firth's penalized logit 4. Bayesian prior: use a weakly informative prior (not available in statsmodels directly)

# Check if dropping a variable resolves convergence
results_reduced = smf.logit("y ~ x1 + x2", data=df).fit()   # dropped x3
print(results_reduced.mle_retvals["converged"])

---

Marginal Effects

GLM and discrete model coefficients are on the link scale (log-odds, log rate). Marginal effects translate them to the response scale (probabilities, rates), enabling direct interpretation.

get_margeff() Method

Available on Logit, Probit, Poisson, MNLogit, and GLM results:

results = smf.logit("y ~ x1 + x2 + C(group)", data=df).fit()

# Average marginal effects (AME) — most common, averaged over observed X
mfx = results.get_margeff(at="overall", method="dydx")
print(mfx.summary())

at Parameter Options

AME (at="overall") is generally preferred over MEM because it does not evaluate the model at an artificial observation (the mean of a categorical variable is not a valid category value).

method Parameter Options

For count models, replace P with the expected count μ.

Dummy Variables in Marginal Effects

For binary/categorical predictors, use dummy=True to compute the discrete change in probability when the variable switches from 0 to 1:

mfx = results.get_margeff(at="overall", method="dydx", dummy=True)
print(mfx.summary())

Without dummy=True, marginal effects for binary variables are computed as a derivative (treating the variable as continuous), which is technically incorrect.

Accessing Marginal Effects Programmatically

mfx = results.get_margeff(at="overall", method="dydx", dummy=True)

mfx.margeff          # Marginal effect point estimates (array)
mfx.margeff_se       # Standard errors
mfx.pvalues          # P-values
mfx.conf_int()       # Confidence intervals — DataFrame with [lower, upper]
mfx.tvalues          # t-statistics

# As a DataFrame
mfx_df = pd.DataFrame({
    "margeff": mfx.margeff,
    "se": mfx.margeff_se,
    "pvalue": mfx.pvalues,
    "ci_lo": mfx.conf_int()[:, 0],
    "ci_hi": mfx.conf_int()[:, 1]
}, index=mfx.margeff_se.index if hasattr(mfx.margeff_se, "index") else None)

Marginal Effects from GLM

get_margeff() is available on GLM results, but the dedicated Logit/Probit/Poisson classes provide additional post-estimation features (e.g., pred_table(), McFadden's pseudo R-squared). Prefer the dedicated classes when available:

# Preferred for marginal effects
results = smf.logit("y ~ x1 + x2", data=df).fit()
mfx = results.get_margeff()

For count models:

results = smf.poisson("count ~ x1 + x2", data=df).fit()
mfx = results.get_margeff(at="overall", method="dydx")
print(mfx.summary())
# For Poisson, dydx marginal effect is ∂E[count]/∂x = β * E[count]

---

Exposure and Offset for Rate Models

When modeling event rates (events per unit of exposure), include exposure in the model so coefficients are interpretable as rate ratios.

exposure Parameter

Pass the raw exposure variable. statsmodels takes ln(exposure) internally:

# Poisson rate model: events per person-year
results = smf.poisson("events ~ x1 + x2", data=df,
                       exposure=df["person_years"]).fit()

# Array API
results = sm.Poisson(df["events"], X, exposure=df["person_years"]).fit()

offset Parameter

Pass the pre-computed log(exposure) manually:

import numpy as np

df["log_pop"] = np.log(df["population"])
results = smf.poisson("events ~ x1 + x2", data=df,
                       offset=df["log_pop"]).fit()

exposure and offset are equivalent: offset = log(exposure). Use offset when your exposure variable is already on the log scale, or when you need a custom offset that is not the log of a simple exposure.

Interpreting Rate Model Coefficients

# Coefficients represent log(rate ratio) per unit change in x
irr = np.exp(results.params)
print(irr)
# Exclude the constant when interpreting IRRs — it represents the baseline rate
# (at x=0 and offset=0), which may not be meaningful on its own

Negative Binomial Rate Model

# Overdispersed rate model
results = sm.NegativeBinomial(
    df["events"], X, exposure=df["person_years"]
).fit()

Rate Prediction

When predicting with an exposure/offset, include it in the prediction call:

# Predict expected counts (not rates) for new observations
new_X = sm.add_constant(new_df[["x1", "x2"]])
pred_counts = results.predict(new_X, exposure=new_df["person_years"])

# Predicted rates (events per unit exposure)
pred_rates = pred_counts / new_df["person_years"]

---

Ordered Logit and Probit (OrderedModel)

For ordinal outcomes (e.g., Likert scales, education levels, satisfaction ratings):

from statsmodels.miscmodels.ordinal_model import OrderedModel

# Ordered logit (proportional odds model)
results = OrderedModel(y, X, distr='logit').fit(method='bfgs')
print(results.summary())

# Ordered probit
results = OrderedModel(y, X, distr='probit').fit(method='bfgs')

• y must be an ordered categorical or integer-coded (0, 1, 2, ..., K-1)
• distr: 'logit' (proportional odds) or 'probit'
• Formula API is also available via OrderedModel.from_formula("y ~ x1 + x2", data=df, distr='logit')
• Returns threshold parameters (cutpoints) alongside regression coefficients
• Coefficients represent log-odds (logit) or z-score (probit) change per unit increase in x
• For interpretation, compute predicted probabilities for each category:

pred_probs = results.predict()  # N x K matrix of predicted probabilities

• Parallel regression assumption (proportional odds): the effect of each predictor is the same across all thresholds. Violations suggest an alternative like multinomial logit
• results.pred_table() for classification accuracy across ordered categories

---

References and Further Reading

• statsmodels GLM documentation: https://www.statsmodels.org/stable/glm.html
• statsmodels GLM families and links: https://www.statsmodels.org/stable/glm.html#families
• statsmodels discrete choice models: https://www.statsmodels.org/stable/discretemod.html
• statsmodels count model examples: https://www.statsmodels.org/stable/examples/notebooks/generated/count_model.html
• statsmodels marginal effects: https://www.statsmodels.org/stable/margeff.html
• statsmodels zero-inflated models: https://www.statsmodels.org/stable/generated/statsmodels.discrete.count_model.ZeroInflatedPoisson.html
• Cameron, A.C. & Trivedi, P.K. (2013). Regression Analysis of Count Data, 2nd ed. Cambridge University Press.
• Long, J.S. & Freese, J. (2014). Regression Models for Categorical Dependent Variables Using Stata, 3rd ed. Stata Press.
• Wooldridge, J.M. (2010). Econometric Analysis of Cross Section and Panel Data, 2nd ed. MIT Press — Chapters 15-16 (count models, binary and multinomial response).
• Hilbe, J.M. (2011). Negative Binomial Regression, 2nd ed. Cambridge University Press.
• McCullagh, P. & Nelder, J.A. (1989). Generalized Linear Models, 2nd ed. Chapman & Hall.

statsmodels Gotchas and Common Pitfalls

A reference for known sharp edges, silent failures, and library-specific behaviors in statsmodels v0.14.6. Each section describes a problem, why it happens, and how to fix it. No methodology content — syntax and library guidance only.

---

Contents

1. Constant Term Not Auto-Added in Array API 2. patsy vs formulaic Formula Parsing 3. Convergence Warnings in MLE Models 4. summary() vs summary2() 5. Large-Sample vs Exact Tests 6. results.predict() Requires exog in Array API 7. statsmodels vs pyfixest: When to Use Which 8. DataFrame Index Issues 9. Categorical Variable Gotchas 10. NaN/Missing Data Handling

---

1. Constant Term Not Auto-Added in Array API

This is the most common statsmodels error. The Array API (sm.OLS) never adds an intercept. If you forget sm.add_constant(), the model is fit through the origin and all coefficient estimates will be wrong — there is no warning.

import statsmodels.api as sm
import numpy as np

y = np.array([1, 2, 3, 4, 5])
X = np.array([[2, 1], [4, 3], [5, 2], [7, 5], [8, 4]])

# WRONG — missing intercept, coefficients will be incorrect
results = sm.OLS(y, X).fit()

# CORRECT — add constant column first
X_with_const = sm.add_constant(X)  # prepends a column of 1s
results = sm.OLS(y, X_with_const).fit()

sm.add_constant(X) prepends a column of ones named const. The column appears first in results.params.

The Formula API never has this problem — it adds an intercept automatically. To suppress it in the formula API, use "y ~ x - 1".

Edge case: If X already contains a column of constants, sm.add_constant() will detect this. Behavior is controlled by the has_constant parameter: 'skip' (default) silently skips, 'raise' raises an error, 'add' adds regardless. When in doubt:

# Check whether constant is present before adding
import numpy as np
if not np.any(np.all(X == X[0, :], axis=0)):
    X = sm.add_constant(X)

---

2. patsy vs formulaic Formula Parsing

statsmodels and pyfixest use different formula parsing libraries. If you are mixing both packages in the same analysis, the formula dialects are not interchangeable.

patsy-specific patterns (statsmodels only)

import statsmodels.formula.api as smf

# Numpy transformations work directly in patsy formulas
results = smf.ols("y ~ np.log(x1) + np.sqrt(x2)", data=df).fit()

# I() protects operators that patsy would otherwise interpret
results = smf.ols("y ~ x1 + I(x1**2)", data=df).fit()

# Explicit reference level for Treatment coding
results = smf.ols("y ~ C(region, Treatment(reference='West'))", data=df).fit()

# Sum coding (deviation from grand mean)
results = smf.ols("y ~ C(region, Sum)", data=df).fit()

# Helmert coding
results = smf.ols("y ~ C(region, Helmert)", data=df).fit()

What does NOT work in patsy (unlike formulaic)

• patsy does not support the | operator for fixed effects — "y ~ x | fe" will

error. Use pyfixest for fixed-effects absorption.

• Some R contrast functions (contr.treatment, contr.sum) are not directly

available; patsy has its own equivalents (Treatment, Sum, Helmert).

Rule of thumb: If transforming a variable before the formula (outside the string), the dialect difference doesn't matter. Conflicts arise only when using advanced formula features — check which library is active before assuming syntax will transfer.

---

3. Convergence Warnings in MLE Models

Models fit by maximum likelihood — Logit, Probit, MNLogit, GLM, ARIMA, and others — can fail to converge. statsmodels reports this as:

ConvergenceWarning: Maximum Likelihood optimization failed to converge.
Check mle_retvals

The results object is still returned but estimates are unreliable.

Common Causes and Fixes

Perfect or quasi-perfect separation (Logit/Probit)

A predictor perfectly or near-perfectly separates outcomes. The log-likelihood has no maximum — the coefficient diverges toward infinity.

Signs: very large coefficient values (e.g., 1e8), extremely small standard errors or extremely large ones, ConvergenceWarning.

# Check: look at coefficient magnitudes
print(results.params)

# Check: examine optimization return values
print(results.mle_retvals)
# 'converged': False indicates the problem

Fixes: remove the separating variable, use Firth's penalized logistic regression (not built into statsmodels — use logistf in R or a third-party package), or collect more data.

Scale issues

Predictors on vastly different scales make the likelihood surface poorly conditioned.

import pandas as pd
from sklearn.preprocessing import StandardScaler

# Standardize predictors before fitting
scaler = StandardScaler()
df_scaled = df.copy()
df_scaled[['x1', 'x2']] = scaler.fit_transform(df[['x1', 'x2']])
results = smf.logit("y ~ x1 + x2", data=df_scaled).fit()

Too many parameters relative to N

Reduce model complexity: drop variables, combine categories, or use regularization.

Optimizer settings

# Try a different optimizer
results = smf.logit("y ~ x1 + x2", data=df).fit(method='bfgs')
results = smf.logit("y ~ x1 + x2", data=df).fit(method='nm')       # Nelder-Mead
results = smf.logit("y ~ x1 + x2", data=df).fit(method='powell')

# Increase iteration limit
results = smf.logit("y ~ x1 + x2", data=df).fit(maxiter=1000)

# Provide starting values to help the optimizer
import numpy as np
start = np.zeros(len(results.params))  # or use OLS coefficients as starting point
results = smf.logit("y ~ x1 + x2", data=df).fit(start_params=start)

Available method values: 'newton' (default), 'nm', 'bfgs', 'lbfgs', 'powell', 'cg', 'ncg', 'basinhopping'.

GLM-specific: For GLM models, convergence failure can also indicate a poor choice of link function or variance structure for the data. Check that the family/link combination makes sense for the response variable.

---

4. summary() vs summary2()

Both methods describe the same model but have different layouts and programmatic accessibility.

summary()

Returns a Summary object. Calling print(results.summary()) renders three panels of formatted text. This is the standard output format widely shown in documentation.

results = smf.ols("y ~ x1 + x2", data=df).fit()
print(results.summary())       # Three-panel text layout

Extracting data from the text output is fragile. Prefer direct attribute access.

summary2()

Returns a Summary2 object. The layout differs — one table for model info, one for coefficients. The coefficient table is returned as a pandas DataFrame, making programmatic access easier.

s2 = results.summary2()
print(s2.tables[0])   # Model info as DataFrame
print(s2.tables[1])   # Coefficients, SEs, t-stats, p-values as DataFrame

# Access coefficient table directly
coef_table = s2.tables[1]
print(coef_table['Coef.'])
print(coef_table['P>|t|'])

Recommended: Direct Attribute Access

For any programmatic extraction of results, skip both summary methods and read attributes directly:

params    = results.params        # Coefficients (Series with variable names)
pvalues   = results.pvalues       # p-values
conf      = results.conf_int()    # CI as DataFrame, columns [0.025, 0.975]
bse       = results.bse           # Standard errors
rsq       = results.rsquared
rsq_adj   = results.rsquared_adj
aic       = results.aic
n         = results.nobs

This approach is robust to formatting changes across statsmodels versions and is not affected by display width or locale settings.

---

5. Large-Sample vs Exact Tests

statsmodels mostly uses asymptotic (large-sample) theory for inference. This matters in small samples where asymptotic approximations can be poor.

What is exact vs asymptotic in statsmodels OLS

Practical implication

With N < 50 and many parameters, Wald tests, LR tests, and robust-SE-based inference may diverge meaningfully from exact p-values. Consider:

• Bootstrapping for inference in small samples
• Exact permutation tests
• Being cautious about borderline p-values (e.g., 0.04 vs 0.06 is not

meaningfully different when N is small)

# Check N relative to parameters before interpreting results
print(f"N = {int(results.nobs)}, K = {len(results.params)}, ratio = {results.nobs/len(results.params):.1f}")
# Asymptotic tests are more reliable when N/K >> 10

---

6. results.predict() Requires exog in Array API

Formula API: pass a DataFrame, variable names handled automatically

results = smf.ols("y ~ x1 + x2", data=df).fit()

new_df = pd.DataFrame({"x1": [3, 6], "x2": [2, 4]})
predictions = results.predict(new_df)   # Works — patsy builds design matrix

Array API: you must replicate the exact design matrix structure

y = df["y"].values
X = df[["x1", "x2"]].values
X_with_const = sm.add_constant(X)
results = sm.OLS(y, X_with_const).fit()

# WRONG — missing constant
new_X = np.array([[3, 2], [6, 4]])
predictions = results.predict(new_X)   # Dimension mismatch error

# CORRECT — add constant to prediction data too
new_X = np.array([[3, 2], [6, 4]])
new_X_with_const = sm.add_constant(new_X)
predictions = results.predict(new_X_with_const)

get_prediction() for confidence and prediction intervals

results.predict() returns point predictions only. To get confidence intervals (interval around the conditional mean) or prediction intervals (interval for individual observations), use results.get_prediction():

# Formula API
pred = results.get_prediction(new_df)
frame = pred.summary_frame(alpha=0.05)  # Returns DataFrame

print(frame.columns)
# ['mean', 'mean_se', 'mean_ci_lower', 'mean_ci_upper',
#  'obs_ci_lower', 'obs_ci_upper']

# mean_ci_*: confidence interval for E[y | x]
# obs_ci_*: prediction interval for individual y

The obs_ci_* columns give the prediction interval, which is always wider than the confidence interval because it accounts for irreducible residual variance.

---

7. statsmodels vs pyfixest vs linearmodels: When to Use Which

Summary rule

• statsmodels for: GLMs, time series, diagnostic tests, mixed models,

discrete choice, models without fixed effects or panel structure

• pyfixest for: fixed effects absorption, IV with FE, DiD, fast panel

estimation, modern robust SE options (CRV3, wild bootstrap)

• linearmodels for: random effects, between/FD estimation, Fama-MacBeth,

IV without FE (LIML, GMM), system estimation (SUR, 3SLS), asset pricing, Driscoll-Kraay SEs

---

8. DataFrame Index Issues

statsmodels aligns inputs by DataFrame index. Non-default or non-contiguous indices can produce unexpected behavior.

Symptoms

• Predictions or residuals contain NaN unexpectedly
• Wrong number of observations after merging DataFrames
• IndexError or shape mismatch when passing new data to predict()

Fix: reset index before fitting

df = df.reset_index(drop=True)
results = smf.ols("y ~ x1 + x2", data=df).fit()

When this bites you most often

# Subsetting a DataFrame preserves original index
subset = df[df['year'] >= 2010]   # index is NOT 0, 1, 2, ...
results = smf.ols("y ~ x1", data=subset).fit()  # May cause issues

# Fix: reset after subsetting
subset = df[df['year'] >= 2010].reset_index(drop=True)
results = smf.ols("y ~ x1", data=subset).fit()

This is especially important when constructing prediction DataFrames that need to align with training data index.

---

9. Categorical Variable Gotchas

Default reference level

C(var) in a patsy formula drops the first level in sorted order as the reference (baseline). This is Treatment coding.

# If region has levels ['East', 'North', 'South', 'West'],
# 'East' is dropped as reference by default
results = smf.ols("y ~ C(region)", data=df).fit()
print(results.params)
# Intercept          ...
# C(region)[T.North] ...
# C(region)[T.South] ...
# C(region)[T.West]  ...

Setting an explicit reference level

# Set 'West' as the reference category
results = smf.ols("y ~ C(region, Treatment(reference='West'))", data=df).fit()

All levels with no intercept (cell-means parameterization)

# Drops intercept; each coefficient is the mean of y for that category
results = smf.ols("y ~ C(region) - 1", data=df).fit()

High-cardinality categoricals

Including a categorical with many levels as dummies is slow and memory-intensive. If the variable has more than ~50 levels and you want to control for it as a fixed effect, use pyfixest instead:

import pyfixest as pf
results = pf.feols("y ~ x1 | region", data=df)

Numeric variables accidentally treated as categorical (or vice versa)

# If 'year' is stored as int but you want it as categorical:
results = smf.ols("y ~ C(year)", data=df).fit()

# If 'code' is stored as string but you want it as numeric:
df['code'] = df['code'].astype(float)
results = smf.ols("y ~ code", data=df).fit()

---

10. NaN/Missing Data Handling

Formula API behavior

The formula API drops rows with NaN in any variable used in the formula (complete case analysis). A warning is issued if rows are dropped:

MissingDataError: No missing values in data, check formula

or, if values are dropped:

/path/to/statsmodels/...: n observations dropped due to missing values

Always verify:

print(f"Input rows: {len(df)}")
print(f"Observations used in model: {int(results.nobs)}")
# If these differ, rows were dropped due to NaN

Array API behavior

The array API does not handle NaN automatically. NaN values in y or X will propagate through linear algebra operations and produce NaN coefficients with no warning.

# WRONG — NaN in array silently corrupts results
results = sm.OLS(y, X).fit()
print(results.params)  # May show NaN for all coefficients

# CORRECT — clean data before fitting
mask = ~(np.isnan(y) | np.any(np.isnan(X), axis=1))
results = sm.OLS(y[mask], X[mask]).fit()

Diagnosing missing data before fitting

import pandas as pd

# Formula API: check which variables have NaN
print(df[['y', 'x1', 'x2']].isna().sum())

# Array API: check arrays directly
print(f"NaN in y: {np.isnan(y).sum()}")
print(f"NaN in X: {np.isnan(X).sum()}")

# After fitting with formula API: confirm row count
print(f"Rows used: {int(results.nobs)} of {len(df)}")

Imputation

statsmodels does not do imputation. If you need to impute before fitting, use sklearn.impute or pandas operations before passing data to statsmodels.

---

References and Further Reading

• statsmodels official documentation: https://www.statsmodels.org/stable/
• statsmodels FAQ: https://www.statsmodels.org/stable/faq.html
• patsy formula language reference: https://patsy.readthedocs.io/en/latest/formulas.html
• patsy categorical coding reference: https://patsy.readthedocs.io/en/latest/categorical-coding.html
• statsmodels regression diagnostics: https://www.statsmodels.org/stable/stats.html#regression-diagnostics
• statsmodels GLM reference: https://www.statsmodels.org/stable/glm.html

statsmodels Hypothesis Testing Reference

Hypothesis testing for regression coefficients, model comparisons, and multiple testing corrections using statsmodels v0.14.6. All examples assume a fitted OLS results object unless otherwise noted.

---

Contents

1. t-Tests on Individual Coefficients 2. F-Tests on Groups of Coefficients 3. Wald Tests 4. Likelihood Ratio Tests 5. Multiple Comparison Corrections 6. Chi-Squared Tests 7. Comparing Nested Models 8. Joint Significance Tests 9. References and Further Reading

---

t-Tests on Individual Coefficients

Standard coefficient tests are reported automatically in results.summary(). To access them programmatically or test non-zero hypotheses:

results = smf.ols("y ~ x1 + x2 + x3", data=df).fit()

# Individual coefficient t-statistics and p-values (all coefficients)
print(results.tvalues)   # Series indexed by variable name
print(results.pvalues)   # Two-sided p-values

# Custom hypothesis on a single coefficient
# H0: beta_x1 = 0
t_test = results.t_test("x1 = 0")
print(t_test)

# H0: beta_x1 = 1
t_test = results.t_test("x1 = 1")
print(t_test)

# Access specific attributes from the test result
print(t_test.tvalue)     # t-statistic
print(t_test.pvalue)     # two-sided p-value
print(t_test.conf_int()) # confidence interval

The string syntax in t_test() uses the variable names from the formula interface. For models fit with sm.OLS(y, X) (array interface), use the R-matrix form described in the Wald Tests section.

---

F-Tests on Groups of Coefficients

Joint Significance of Multiple Coefficients

# H0: beta_x2 = beta_x3 = 0
f_test = results.f_test("x2 = 0, x3 = 0")
print(f_test)            # Prints F-statistic, p-value, df numerator

# Access attributes directly
print(f_test.fvalue)     # F-statistic (scalar or array)
print(f_test.pvalue)     # p-value
print(f_test.df_num)     # numerator degrees of freedom
print(f_test.df_denom)   # denominator degrees of freedom

R-Matrix Notation (Array API)

For models fit without formula interface or when constructing hypotheses programmatically:

import numpy as np

# Coefficient order matches results.params index
# Intercept, x1, x2, x3 → indices 0, 1, 2, 3

R = np.array([[0, 0, 1, 0],   # coefficient on x2 = 0
              [0, 0, 0, 1]])   # coefficient on x3 = 0
f_test = results.f_test(R)

Testing Equality of Coefficients

# H0: beta_x1 = beta_x2
f_test = results.f_test("x1 = x2")
# Equivalent: "x1 - x2 = 0"
f_test = results.f_test("x1 - x2 = 0")

---

Wald Tests

General linear hypothesis testing using asymptotic chi-squared or F distribution:

import numpy as np

# Wald test (chi-squared by default)
wald = results.wald_test("x1 = 0, x2 = 0")
print(wald)

# Wald test with F-distribution (preferred for OLS in finite samples)
wald = results.wald_test("x1 = 0, x2 = 0", use_f=True)

# Complex linear restrictions
# H0: beta_x1 + beta_x2 = 1 AND beta_x3 = 0
wald = results.wald_test("x1 + x2 = 1, x3 = 0")

# R-matrix form with constraint vector q (Rβ = q)
R = np.array([[0, 1, 1, 0],   # x1 + x2 = 1
              [0, 0, 0, 1]])   # x3 = 0
q = np.array([1, 0])
wald = results.wald_test((R, q))

# Access attributes
print(wald.statistic)   # chi-squared or F statistic
print(wald.pvalue)
print(wald.df_num)      # degrees of freedom (numerator for F)

Key distinction between f_test and wald_test:

• f_test always returns an F-statistic
• wald_test returns chi-squared by default; use use_f=True for F
• For OLS, f_test and wald_test(..., use_f=True) give identical results
• Wald tests are asymptotic — use F-version for finite samples when available

---

Likelihood Ratio Tests

Compare nested models using the log-likelihood ratio. The LR statistic (-2 × log-likelihood ratio) follows a chi-squared distribution with degrees of freedom equal to the number of restrictions.

For OLS Models

LR tests are not directly available as a statsmodels function for OLS; compute manually:

from scipy.stats import chi2

# Fit restricted and unrestricted models
results_full = smf.ols("y ~ x1 + x2 + x3", data=df).fit()
results_restricted = smf.ols("y ~ x1", data=df).fit()

# Compute LR statistic
lr_stat = -2 * (results_restricted.llf - results_full.llf)
df_diff = results_full.df_model - results_restricted.df_model

lr_pvalue = chi2.sf(lr_stat, df_diff)

print(f"LR statistic: {lr_stat:.4f}")
print(f"p-value:      {lr_pvalue:.4f}")
print(f"df:           {df_diff}")

results.llf is the log-likelihood value (available on all statsmodels results objects that support MLE).

For MLE Models (Logit, Probit, GLM, ARIMA)

The same pattern applies — all MLE results objects expose .llf:

results_full = smf.logit("y ~ x1 + x2 + x3", data=df).fit()
results_restricted = smf.logit("y ~ x1", data=df).fit()

lr_stat = -2 * (results_restricted.llf - results_full.llf)
df_diff = (results_full.df_model - results_restricted.df_model)
lr_pvalue = chi2.sf(lr_stat, df_diff)

LR tests require models to be fit on exactly the same observations. Verify with results.nobs before computing.

---

Multiple Comparison Corrections

multipletests — Adjust P-Values for Multiple Testing

from statsmodels.stats.multitest import multipletests

pvals = [0.01, 0.04, 0.03, 0.15, 0.50]

# Returns: (reject array, corrected p-values, corrected alpha sidak, corrected alpha bonf)
reject, pvals_corrected, alphac_sidak, alphac_bonf = multipletests(
    pvals, alpha=0.05, method='fdr_bh'
)

print(reject)           # Boolean array: True where H0 is rejected after correction
print(pvals_corrected)  # Adjusted p-values

Available Methods

FWER (Family-Wise Error Rate): controls probability of any false positive. FDR (False Discovery Rate): controls expected proportion of false positives among rejections.

Choosing a Method

• Few tests (< 20): Holm-Bonferroni ('holm') is a safe default — step-down

procedure dominates Bonferroni without additional assumptions

• Many tests (> 20): Benjamini-Hochberg ('fdr_bh') is the standard choice —

valid under independence or positive dependence (PRDS condition)

• Unknown or arbitrary dependence structure: Benjamini-Yekutieli ('fdr_by')

is conservative but guaranteed valid

• Very conservative needed (confirmatory single-family): Bonferroni

('bonferroni')

Extracting Corrected P-Values Into a DataFrame

results_df = pd.DataFrame({
    'variable': ['x1', 'x2', 'x3', 'x4', 'x5'],
    'pvalue_raw': pvals
})

reject, pvals_corrected, _, _ = multipletests(
    results_df['pvalue_raw'], alpha=0.05, method='fdr_bh'
)

results_df['pvalue_corrected'] = pvals_corrected
results_df['reject_h0'] = reject
print(results_df)

---

Chi-Squared Tests

These tests come from scipy but are commonly used alongside statsmodels analyses.

Goodness of Fit

from scipy.stats import chisquare

# Compare observed frequencies to expected
stat, pvalue = chisquare(observed_freq, f_exp=expected_freq)
# H0: observed frequencies match expected distribution
# If f_exp is omitted, tests for uniform distribution

Independence (Contingency Table)

from scipy.stats import chi2_contingency

table = pd.crosstab(df['var1'], df['var2'])
chi2, p, dof, expected = chi2_contingency(table)

# expected: array of expected cell counts under independence
# Rule of thumb: chi2_contingency is unreliable if any expected cell < 5

For 2x2 tables with small expected counts, use Fisher's exact test instead:

from scipy.stats import fisher_exact

odds_ratio, pvalue = fisher_exact(table_2x2)

---

Comparing Nested Models

Using AIC and BIC

AIC and BIC are available on all statsmodels results objects. Lower values indicate better fit (penalizing for model complexity):

model1 = smf.ols("y ~ x1", data=df).fit()
model2 = smf.ols("y ~ x1 + x2", data=df).fit()
model3 = smf.ols("y ~ x1 + x2 + x3", data=df).fit()

comparison = pd.DataFrame({
    'Model': ['x1', 'x1+x2', 'x1+x2+x3'],
    'AIC': [m.aic for m in [model1, model2, model3]],
    'BIC': [m.bic for m in [model1, model2, model3]],
    'R2_adj': [m.rsquared_adj for m in [model1, model2, model3]],
    'nparams': [m.df_model + 1 for m in [model1, model2, model3]]
})
print(comparison)

AIC vs BIC: BIC penalizes complexity more heavily (penalty = log(n) per parameter vs 2 per parameter for AIC) — BIC is preferred when model selection for inference is the goal; AIC when predictive accuracy is the goal.

Using ANOVA (Nested OLS Models)

from statsmodels.stats.anova import anova_lm

# Models must be fit with formula API (smf.ols)
# Pass models in order from most restricted to least restricted
anova_table = anova_lm(model1, model2, model3)
print(anova_table)
# Columns: df_resid, ssr, df_diff, ss_diff, F, Pr(>F)

Each row tests the improvement from adding the variables in the current model relative to the previous model. Significant F indicates the added variables improve fit beyond chance.

Requirements for anova_lm:

• Models must be nested (each model is a restriction of the next)
• Models must be fit on the same observations (nobs must match)
• Models must be fit using the formula interface

---

Joint Significance Tests

The overall model F-test (H0: all slope coefficients = 0) is always reported in results.summary(). Access it programmatically:

# Overall F-test: H0 is that all slopes jointly equal zero
print(f"F-statistic: {results.fvalue:.4f}")
print(f"F p-value:   {results.f_pvalue:.6f}")
print(f"df model:    {results.df_model}")    # numerator df
print(f"df resid:    {results.df_resid}")    # denominator df

# Test a specific subset of variables jointly
# H0: coefficients on x2 AND x3 are jointly zero
f_test = results.f_test("x2 = 0, x3 = 0")
print(f"F: {f_test.fvalue:.4f}, p: {f_test.pvalue:.4f}")

For interaction models, testing the joint significance of an interaction term and its constituent main effects:

# Model with interaction
results_int = smf.ols("y ~ x1 + x2 + x1:x2", data=df).fit()

# H0: x2 and x1:x2 are jointly zero (testing whether x2 matters at all)
f_test = results_int.f_test("x2 = 0, x1:x2 = 0")
print(f_test)

Note: In formula syntax, interaction terms use : notation in f_test and wald_test strings, matching how they appear in results.params.

---

References and Further Reading

• statsmodels hypothesis tests: https://www.statsmodels.org/stable/stats.html
• statsmodels multitest: https://www.statsmodels.org/stable/generated/statsmodels.stats.multitest.multipletests.html
• statsmodels contrast and test documentation: https://www.statsmodels.org/stable/contrasts.html
• Greene, W.H. (2018). Econometric Analysis, 8th ed. — Chapters 5-6

(hypothesis testing)

• Wooldridge, J.M. (2019). Introductory Econometrics, 7th ed. — Chapter 4

(multiple regression: inference)

• Benjamini & Hochberg (1995). "Controlling the False Discovery Rate." JRSS-B.

statsmodels Linear Models Reference

statsmodels v0.14.6 — syntax and library guidance only.

---

Contents

1. OLS (Ordinary Least Squares) 2. WLS (Weighted Least Squares) 3. GLS (Generalized Least Squares) 4. Robust Regression (RLM) 5. Quantile Regression (QuantReg) 6. Polynomial and Spline Regression 7. Interaction Terms 8. Mixed Effects / Multilevel Models (MixedLM) 9. References and Further Reading

---

OLS (Ordinary Least Squares)

Formula API

import statsmodels.formula.api as smf

results = smf.ols("y ~ x1 + x2 + x3", data=df).fit()
print(results.summary())

Array API

import statsmodels.api as sm

X = sm.add_constant(df[["x1", "x2", "x3"]])
results = sm.OLS(df["y"], X).fit()
print(results.summary())

sm.add_constant() prepends a column of ones to the design matrix. Without it, the model has no intercept. The constant column is named "const" by default.

Accessing Coefficients

# Point estimates
results.params          # Series indexed by variable name

# P-values (two-tailed)
results.pvalues

# Confidence intervals (default alpha=0.05)
results.conf_int()              # 95% CI — returns DataFrame with columns [0.025, 0.975]
results.conf_int(alpha=0.01)    # 99% CI

Example — extract a single coefficient with its CI:

coef = results.params["x1"]
ci_lo, ci_hi = results.conf_int().loc["x1"]
print(f"x1: {coef:.4f} [{ci_lo:.4f}, {ci_hi:.4f}]")

Model Fit Statistics

results.rsquared        # R²
results.rsquared_adj    # Adjusted R²
results.aic             # Akaike Information Criterion
results.bic             # Bayesian Information Criterion
results.fvalue          # F-statistic (overall model significance)
results.f_pvalue        # p-value for F-test
results.nobs            # Number of observations
results.df_resid        # Residual degrees of freedom
results.df_model        # Model degrees of freedom (excl. constant)
results.mse_resid       # Mean squared error of residuals (= σ²)

Fitted Values and Residuals

results.fittedvalues    # ŷ — Series, same index as training data
results.resid           # ε = y - ŷ (raw residuals)
results.resid_pearson   # Pearson (standardized) residuals for OLS same as resid / RMSE

Studentized residuals (externally studentized, for outlier detection):

influence = results.get_influence()
student_resid = influence.resid_studentized_external  # values > |2| flag potential outliers

Prediction

In-sample fitted values:

results.fittedvalues

Out-of-sample point prediction:

new_data = pd.DataFrame({"x1": [1.5, 2.0], "x2": [3.0, 4.5], "x3": [0.1, 0.2]})
y_pred = results.predict(new_data)   # returns Series

Out-of-sample prediction with confidence and prediction intervals:

pred = results.get_prediction(new_data)
frame = pred.summary_frame(alpha=0.05)
# Columns: mean, mean_se, mean_ci_lower, mean_ci_upper, obs_ci_lower, obs_ci_upper
# mean_ci_*: confidence interval for the mean response
# obs_ci_*: prediction interval for an individual observation
print(frame)

Array API — pass an array with the same column structure as the training X:

X_new = sm.add_constant(pd.DataFrame({"x1": [1.5], "x2": [3.0], "x3": [0.1]}))
y_pred = results.predict(X_new)

Saving and Reloading Results

import pickle

# Save
with open("ols_results.pkl", "wb") as f:
    pickle.dump(results, f)

# Load
with open("ols_results.pkl", "rb") as f:
    results = pickle.load(f)

---

WLS (Weighted Least Squares)

Use when errors have non-constant variance (heteroskedasticity) but the variance structure is known or estimable. Observations with smaller error variance receive higher weight.

Formula API

# weight = 1 / variance: high variance → low weight → less influence
results = smf.wls("y ~ x1 + x2", data=df, weights=1/df["variance"]).fit()

Array API

X = sm.add_constant(df[["x1", "x2"]])
results = sm.WLS(df["y"], X, weights=1/df["variance"]).fit()

Common Weight Specifications

# Population-weighted regression (each obs represents pop_i people)
results = smf.wls("y ~ x1 + x2", data=df, weights=df["population"]).fit()

# Precision weights from measurement error variance
results = smf.wls("y ~ x1", data=df, weights=1/df["meas_var"]).fit()

# Inverse of squared fitted standard deviations (feasible WLS)
ols_resid = smf.ols("y ~ x1", data=df).fit().resid
df["abs_resid"] = ols_resid.abs()
aux = smf.ols("abs_resid ~ x1", data=df).fit()
fitted_sd = aux.fittedvalues
results = smf.wls("y ~ x1", data=df, weights=1/fitted_sd**2).fit()

The weights parameter in statsmodels WLS is the inverse of the error variance, not the square root. Larger weight = observation is more reliable.

Accessing Results

All standard OLS result attributes apply: .params, .pvalues, .conf_int(), .fittedvalues, .resid, .rsquared, etc.

The R² in WLS is weighted; interpret with care relative to OLS R².

---

GLS (Generalized Least Squares)

Use when errors are correlated or have non-constant variance and the full covariance structure of errors is known (or estimable).

Array API

import statsmodels.api as sm
import numpy as np

X = sm.add_constant(df[["x1", "x2"]])
y = df["y"]

# sigma: covariance matrix of errors (n x n)
results = sm.GLS(y, X, sigma=covariance_matrix).fit()

sigma can be:

• A full (n, n) covariance matrix
• A 1-d array of length n for diagonal (heteroskedastic only, equivalent to WLS)
• None (falls back to OLS)

GLSAR for AR(p) Error Structure

When errors follow an autoregressive process, GLSAR is more practical than specifying the full covariance matrix:

# AR(1) errors
results = sm.GLSAR(y, X, rho=1).fit()

# AR(2) errors
results = sm.GLSAR(y, X, rho=2).fit()

Iterative estimation (Cochrane-Orcutt / Prais-Winsten style):

from statsmodels.regression.linear_model import GLSAR

model = GLSAR(y, X, rho=1)
results = model.iterative_fit(maxiter=10)
print(f"Estimated rho: {model.rho}")
print(results.summary())

GLS vs. WLS vs. GLSAR

---

Robust Regression (RLM)

M-estimators that iteratively downweight influential observations (outliers). Use when OLS residual diagnostics show influential outliers but removing them is not justified.

Basic Usage

import statsmodels.api as sm

X = sm.add_constant(df[["x1", "x2"]])
y = df["y"]

# Default: Huber's T norm
results = sm.RLM(y, X).fit()
print(results.summary())

Formula API (via smf):

import statsmodels.formula.api as smf

results = smf.rlm("y ~ x1 + x2", data=df).fit()

M-Estimator Norms

All norms are in sm.robust.norms:

# Huber's T (default, tuning constant 1.345)
results = sm.RLM(y, X, M=sm.robust.norms.HuberT()).fit()

# Tukey's biweight (fully rejects extreme outliers)
results = sm.RLM(y, X, M=sm.robust.norms.TukeyBiweight()).fit()

# Andrew's wave (sinusoidal redescending)
results = sm.RLM(y, X, M=sm.robust.norms.AndrewWave()).fit()

# Hampel (three-part influence function)
results = sm.RLM(y, X, M=sm.robust.norms.Hampel()).fit()

# LeastSquares (equivalent to OLS — comparison baseline)
results = sm.RLM(y, X, M=sm.robust.norms.LeastSquares()).fit()

Accessing Weights and Diagnostics

results.weights         # Per-observation weights (0 = completely downweighted)
results.params          # Coefficient estimates
results.pvalues         # P-values
results.conf_int()      # Confidence intervals

Identifying downweighted observations:

import pandas as pd

weight_df = pd.DataFrame({
    "weight": results.weights,
    "resid": results.resid
}, index=df.index)

# Observations with weight below 0.5 are substantially downweighted
outliers = weight_df[weight_df["weight"] < 0.5]
print(f"{len(outliers)} observations downweighted below 0.5")
print(outliers.sort_values("weight"))

Comparing OLS and RLM

ols_res = sm.OLS(y, X).fit()
rlm_res = sm.RLM(y, X).fit()

comparison = pd.DataFrame({
    "OLS": ols_res.params,
    "RLM": rlm_res.params,
    "diff": rlm_res.params - ols_res.params
})
print(comparison)
# Large differences signal that outliers are influencing OLS estimates

---

Quantile Regression (QuantReg)

Models conditional quantiles of y rather than the conditional mean. Useful for understanding distributional effects, heterogeneous treatment effects, and when outliers make OLS mean-unstable.

Basic Usage

# Formula API — median regression
results = smf.quantreg("y ~ x1 + x2", data=df).fit(q=0.5)

# 25th percentile
results = smf.quantreg("y ~ x1 + x2", data=df).fit(q=0.25)

# 90th percentile
results = smf.quantreg("y ~ x1 + x2", data=df).fit(q=0.90)

# Array API
results = sm.QuantReg(y, X).fit(q=0.75)

The q parameter is the quantile to estimate (0 < q < 1). q=0.5 is least absolute deviations (LAD) regression, equivalent to median regression.

Standard Errors

By default, sparsity-based asymptotic standard errors are used. Bootstrap standard errors are often preferred in practice:

# Kernel-based robust standard errors (default)
results = smf.quantreg("y ~ x1 + x2", data=df).fit(
    q=0.5,
    vcov="robust",    # default; kernel-based (Powell sandwich)
    kernel="epa",     # Epanechnikov kernel for sparsity estimation
    bandwidth="hsheather"  # bandwidth selection method
)

Available kernels for bandwidth: "epa" (Epanechnikov), "cos", "gau" (Gaussian), "par" (Parzen), "tri" (triangular).

Fitting Multiple Quantiles

quantiles = [0.10, 0.25, 0.50, 0.75, 0.90]
qr_results = {}

for q in quantiles:
    qr_results[q] = smf.quantreg("y ~ x1 + x2", data=df).fit(q=q)

# Compare x1 coefficient across quantiles
x1_coefs = pd.Series(
    {q: qr_results[q].params["x1"] for q in quantiles},
    name="x1_coef"
)
print(x1_coefs)

Accessing Results

All standard result attributes apply:

results.params
results.pvalues
results.conf_int()
results.fittedvalues
results.resid
results.predict(new_data)

Pseudo R² (Koenker-Machado):

results.prsquared   # 1 - sum(check_loss) / sum(check_loss_null)

---

Polynomial and Spline Regression

Polynomial Terms via I()

The I() wrapper in patsy formulas passes the expression through as raw Python, bypassing formula operator interpretation:

# Quadratic (degree 2)
results = smf.ols("y ~ x1 + I(x1**2)", data=df).fit()

# Cubic (degree 3)
results = smf.ols("y ~ x1 + I(x1**2) + I(x1**3)", data=df).fit()

# Interaction with polynomial
results = smf.ols("y ~ x1 + I(x1**2) + x2 + x1:x2", data=df).fit()

In the array API, add polynomial columns to the DataFrame manually:

X = pd.DataFrame({
    "x1": df["x1"],
    "x1_sq": df["x1"] ** 2,
    "x1_cu": df["x1"] ** 3,
    "x2": df["x2"]
})
X = sm.add_constant(X)
results = sm.OLS(df["y"], X).fit()

B-Spline Basis via patsy bs()

from patsy import bs

# Cubic B-spline with 4 degrees of freedom (1 intercept + 3 basis functions)
results = smf.ols("y ~ bs(x1, df=4, degree=3)", data=df).fit()

# Natural spline with knots at specific values
results = smf.ols("y ~ bs(x1, knots=[25, 50, 75], degree=3)", data=df).fit()

# Spline without intercept (no separate constant in the spline)
results = smf.ols("y ~ bs(x1, df=5, degree=3, include_intercept=False)", data=df).fit()

df (degrees of freedom): number of basis functions, equal to degree + number_of_interior_knots. Knot locations default to quantiles of x1.

degree:

• 1 = piecewise linear
• 2 = piecewise quadratic
• 3 = cubic (most common, C2 continuity)

Generating the spline basis for prediction or inspection:

from patsy import dmatrix

# Fit on training data
basis_train = dmatrix("bs(x1, df=4, degree=3)", data=df, return_type="dataframe")
results = sm.OLS(df["y"], basis_train).fit()

# Apply same transformation to new data
basis_new = dmatrix("bs(x1, df=4, degree=3)", data=new_df, return_type="dataframe")
y_pred = results.predict(basis_new)

Using dmatrix ensures knot locations from training are reused at prediction time.

---

Interaction Terms

Formula API Operators

# Full interaction
results = smf.ols("y ~ x1 * x2", data=df).fit()

# Interaction term only (suppresses main effects)
results = smf.ols("y ~ x1:x2", data=df).fit()

# Categorical interaction — C() marks a column as categorical
results = smf.ols("y ~ x1 * C(group)", data=df).fit()

# Three-way interaction (all lower-order terms included)
results = smf.ols("y ~ x1 * x2 * x3", data=df).fit()

# All pairwise interactions among three variables
results = smf.ols("y ~ (x1 + x2 + x3)**2", data=df).fit()

Categorical Reference Level

By default, patsy uses the first category alphabetically as the reference level. Override:

# Set "control" as the reference level
results = smf.ols("y ~ x1 * C(group, Treatment('control'))", data=df).fit()

Array API — Manual Interaction Columns

X = df[["x1", "x2"]].copy()
X["x1_x2"] = X["x1"] * X["x2"]   # continuous × continuous

# Continuous × categorical (dummy)
X["treat"] = (df["group"] == "treatment").astype(int)
X["x1_treat"] = X["x1"] * X["treat"]

X = sm.add_constant(X)
results = sm.OLS(df["y"], X).fit()

Interpreting Interaction Coefficients

results = smf.ols("y ~ x1 * x2", data=df).fit()

# Coefficient on x1:x2 is the marginal effect of x1 on y
# per unit increase in x2 (and vice versa)
print(results.params)
print(results.conf_int())

# For categorical interactions: the coefficient on x1:C(group)[T.treated]
# is the difference in the x1 slope between treated and reference groups

---

Mixed Effects / Multilevel Models (MixedLM)

For hierarchical/nested data (e.g., students in schools, repeated measures):

import statsmodels.formula.api as smf

# Random intercept model (varying intercept by group)
results = smf.mixedlm(
    "score ~ ses + school_type",
    data=df,
    groups=df["school_id"]
).fit()
print(results.summary())

# Random intercept + random slope
results = smf.mixedlm(
    "score ~ ses + school_type",
    data=df,
    groups=df["school_id"],
    re_formula="~ses"           # ses has a random slope by school
).fit()

Array API:

import statsmodels.api as sm

model = sm.MixedLM(
    endog=y,
    exog=X,                     # fixed effects design matrix
    groups=groups,              # grouping variable
    exog_re=Z                   # random effects design matrix
)
results = model.fit()

Key parameters:

• groups: grouping variable defining the clusters (e.g., school ID, subject ID)
• re_formula: formula for random effects (default: random intercept only, "~1")
• vc_formula: variance components formula for crossed random effects

Results attributes:

• results.fe_params — fixed effects coefficients
• results.random_effects — dict of random effects by group
• results.cov_re — random effects covariance matrix
• results.bse_fe — standard errors for fixed effects
• results.bse_re — standard errors for random effects variance components

Prediction:

# Fixed effects prediction (population average)
pred = results.predict(new_data)

# Group-specific prediction (includes random effects)
pred = results.predict(new_data, exog_re=new_Z)

Common use cases in education research:

• Students nested in classrooms nested in schools (two-level or three-level)
• Repeated test scores for the same student over time (longitudinal)
• Schools nested in districts with district-level predictors

Convergence: MixedLM uses gradient-based optimization by default (BFGS, L-BFGS-B, CG sequence). For convergence issues, try:

results = model.fit(method=["powell"])     # Powell optimizer
results = model.fit(maxiter=200)           # Increase max iterations
results = model.fit(reml=False)            # ML instead of REML

---

References and Further Reading

• statsmodels OLS documentation: https://www.statsmodels.org/stable/regression.html
• statsmodels WLS documentation: https://www.statsmodels.org/stable/regression.html#wls
• statsmodels GLS documentation: https://www.statsmodels.org/stable/regression.html#gls
• statsmodels RLM documentation: https://www.statsmodels.org/stable/rlm.html
• statsmodels QuantReg documentation: https://www.statsmodels.org/stable/generated/statsmodels.regression.quantile_regression.QuantReg.html
• patsy formula language: https://patsy.readthedocs.io/en/latest/formulas.html
• patsy splines: https://patsy.readthedocs.io/en/latest/spline-regression.html
• Greene, W.H. (2018). Econometric Analysis, 8th ed. — Chapters 4-9 (OLS, GLS, WLS, robust)
• Wooldridge, J.M. (2019). Introductory Econometrics, 7th ed. — Chapters 4, 8 (OLS inference, heteroskedasticity)
• Koenker, R. (2005). Quantile Regression. Cambridge University Press.
• Maronna, R.A., Martin, R.D., Yohai, V.J., & Salibian-Barrera, M. (2019). Robust Statistics: Theory and Methods, 2nd ed.

statsmodels Quickstart Reference

A technical reference for using statsmodels v0.14.6 in Python. Covers installation, the two main APIs, model fitting, reading output, and key library patterns. No methodology content — syntax and library guidance only.

Contents

• Installation
• Formula API vs Array API
• Your First Model (Formula API)
• Your First Model (Array API)
• Reading the Summary Output
• Key Results Attributes
• Robust Standard Errors
• Formula Syntax Overview
• Quick Comparison: statsmodels vs pyfixest
• References and Further Reading

---

Installation

pip install statsmodels

Key dependencies (installed automatically):

Verify installation:

import statsmodels
print(statsmodels.__version__)  # Should print 0.14.x

---

Formula API vs Array API

statsmodels exposes two parallel interfaces. You will use both in different situations — understanding the distinction upfront prevents the most common errors.

Formula API

import statsmodels.formula.api as smf

• Accepts R-style formula strings: "y ~ x1 + x2"
• Uses patsy to parse formulas and build design matrices
• Automatically adds an intercept (constant term) — no manual step needed
• Categorical variables handled via C(var) in the formula
• Variable names come from a data= DataFrame argument

Array API

import statsmodels.api as sm

• Accepts numpy arrays or pandas DataFrames directly
• Does NOT add an intercept automatically — you must add it manually
• More verbose but gives full control over design matrices
• Useful when building matrices programmatically

Same Model, Both Ways

import statsmodels.formula.api as smf
import statsmodels.api as sm
import pandas as pd
import numpy as np

df = pd.DataFrame({
    "y":  [1, 2, 3, 4, 5],
    "x1": [2, 4, 5, 7, 8],
    "x2": [1, 3, 2, 5, 4],
})

# --- Formula API ---
results_f = smf.ols("y ~ x1 + x2", data=df).fit()

# --- Array API ---
y = df["y"].values
X = df[["x1", "x2"]].values
X = sm.add_constant(X)   # CRITICAL: adds a leading column of ones for the intercept
results_a = sm.OLS(y, X).fit()

# Both produce identical coefficient estimates
print(results_f.params)
print(results_a.params)

When to Prefer Each

---

Your First Model (Formula API)

import statsmodels.formula.api as smf
import pandas as pd

df = pd.DataFrame({
    "y":  [1, 2, 3, 4, 5],
    "x1": [2, 4, 5, 7, 8],
    "x2": [1, 3, 2, 5, 4],
})

results = smf.ols("y ~ x1 + x2", data=df).fit()
print(results.summary())

The .fit() call runs OLS and returns a RegressionResultsWrapper. The object holds all estimated parameters, test statistics, and diagnostics. The original model specification is accessible via results.model.

---

Your First Model (Array API)

import statsmodels.api as sm
import numpy as np

y = np.array([1, 2, 3, 4, 5])
X = np.array([[2, 1],
              [4, 3],
              [5, 2],
              [7, 5],
              [8, 4]])
X = sm.add_constant(X)  # CRITICAL: adds intercept column as the first column

results = sm.OLS(y, X).fit()
print(results.summary())

sm.add_constant(X) prepends a column of ones. Without it, the model has no intercept and coefficient estimates will be incorrect.

---

Reading the Summary Output

results.summary() prints a three-panel text table. Here is what each section contains.

Panel 1: Model Information

OLS Regression Results
==============================================================================
Dep. Variable:                      y   R-squared:                       0.986
Model:                            OLS   Adj. R-squared:                  0.972
Method:                 Least Squares   F-statistic:                     70.84
Date:                ...               Prob (F-statistic):              0.0137
Time:                ...               Log-Likelihood:                 -1.3485
No. Observations:                   5   AIC:                             8.697
Df Residuals:                       2   BIC:                             7.883
Df Model:                           2
Covariance Type:            nonrobust

Key fields:

Panel 2: Coefficient Table

==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     -0.3125      0.593     -0.527      0.653      -2.864       2.239
x1             0.5781      0.095      6.092      0.026       0.169       0.987
x2            -0.1563      0.165     -0.947      0.445      -0.867       0.554
==============================================================================

Panel 3: Diagnostic Statistics

==============================================================================
Omnibus:                        0.276   Durbin-Watson:                   2.499
Prob(Omnibus):                  0.871   Jarque-Bera (JB):                0.370
Skew:                          -0.266   Prob(JB):                        0.831
Kurtosis:                       1.812   Cond. No.                         33.2
==============================================================================

summary2()

results.summary2() returns a Summary2 object with a different layout. Access its tables programmatically:

s2 = results.summary2()
print(s2.tables[0])   # Model info table
print(s2.tables[1])   # Coefficient table as a DataFrame

For automated extraction of coefficients and p-values, prefer direct attribute access (results.params, results.pvalues) over parsing summary text.

---

Key Results Attributes

All attributes listed below are available on the RegressionResultsWrapper returned by .fit().

# Accessing attributes directly
print(results.params)
print(results.pvalues)
print(results.rsquared)
print(results.conf_int())           # 95% CI by default
print(results.conf_int(alpha=0.01)) # 99% CI

---

Robust Standard Errors

Pass cov_type to .fit() to override the default OLS standard errors.

import statsmodels.formula.api as smf

# HC1 — White's heteroskedasticity-consistent (good default for cross-sectional data)
results = smf.ols("y ~ x1 + x2", data=df).fit(cov_type='HC1')

# HC3 — bias-corrected; preferred for small samples
results = smf.ols("y ~ x1 + x2", data=df).fit(cov_type='HC3')

# Clustered standard errors
results = smf.ols("y ~ x1", data=df).fit(
    cov_type='cluster',
    cov_kwds={'groups': df['cluster_var']}
)

# HAC — Newey-West heteroskedasticity and autocorrelation consistent (time series)
results = smf.ols("y ~ x1", data=df).fit(
    cov_type='HAC',
    cov_kwds={'maxlags': 4}
)

Covariance types available:

The cov_type only changes standard errors and test statistics — point estimates (coefficients) are unaffected.

---

Formula Syntax Overview

statsmodels uses patsy to parse formula strings. Patsy syntax closely follows R's formula notation.

Common Patterns

Intercept Behavior

• Intercept is included by default: "y ~ x1" fits y = a + b*x1
• Suppress with -1 or +0: "y ~ x1 - 1" fits y = b*x1 through origin
• The array API never adds an intercept — use sm.add_constant(X) explicitly

Categorical Variables in Formulas

import statsmodels.formula.api as smf

# Default: first alphabetical level dropped as reference
results = smf.ols("y ~ C(region)", data=df).fit()

# Explicit reference level
results = smf.ols("y ~ C(region, Treatment(reference='West'))", data=df).fit()

# All levels, no intercept (for a cell-means parameterization)
results = smf.ols("y ~ C(region) - 1", data=df).fit()

In-Formula Transformations

import numpy as np
import statsmodels.formula.api as smf

# Log transformation — numpy functions work directly in patsy formulas
results = smf.ols("y ~ np.log(income) + age", data=df).fit()

# Polynomial — use I() to protect ** from patsy interpretation
results = smf.ols("y ~ age + I(age**2)", data=df).fit()

# Standardize inline (center and scale)
results = smf.ols("y ~ scale(income)", data=df).fit()

---

Quick Comparison: statsmodels vs pyfixest

When to Use Which

Use statsmodels when:

• Fitting GLMs (logit, probit, Poisson, negative binomial, gamma, etc.)
• Working with time series (ARIMA, VAR, state space models)
• Running diagnostic tests (heteroskedasticity, normality, autocorrelation)
• Estimating discrete choice models (MNLogit, ordered logit)
• You need marginal effects from nonlinear models
• The model does not involve multi-way fixed effects

Use pyfixest when:

• Multi-way fixed effects are required (absorbs FE via demeaning — much faster)
• Running IV regressions with fixed effects
• Implementing difference-in-differences or event study designs
• You need CRV3 or wild bootstrap clustered standard errors
• Estimating many models in a loop over panel data

Both libraries support heteroskedasticity-robust and clustered standard errors for OLS. For plain OLS without FE, either works — statsmodels has more post-estimation diagnostics.

---

References and Further Reading

• statsmodels official documentation: https://www.statsmodels.org/stable/
• statsmodels GitHub repository: https://github.com/statsmodels/statsmodels
• Seabold, S. & Perktold, J. (2010). "statsmodels: Econometric and Statistical

Modeling with Python." Proceedings of the 9th Python in Science Conference.

• patsy formula documentation: https://patsy.readthedocs.io/
• patsy formula language reference: https://patsy.readthedocs.io/en/latest/formulas.html
• statsmodels regression example gallery: https://www.statsmodels.org/stable/examples/index.html#regression