statsmodels

Statistical modeling with Python: linear and generalized linear regression, ANOVA, hypothesis tests, model diagnostics, and classic time-series models — with the full summary() tables statisticians expect.

Table of Contents

• Installation / Setup
• Core Syntax
• Common Operations
• Practical Examples
• Common Patterns
• Best Practices
• Common Mistakes
• Troubleshooting
• Related Cheat Sheets
• Official Documentation

Installation / Setup

pip install statsmodels      # or: uv add statsmodels

statsmodels builds on pandas and NumPy. Two APIs:

import statsmodels.api as sm        # array/formula-free API
import statsmodels.formula.api as smf   # formula API: "y ~ x1 + x2"

The formula API is most convenient with a DataFrame. Formulas use patsy syntax.

Core Syntax

Linear regression (OLS)

# Array API: add the intercept column yourself
X = sm.add_constant(df[["x1", "x2"]])     # add_const_constant -> constant column
model = sm.OLS(df["y"], X).fit()

# Formula API: intercept handled automatically
model = smf.ols("y ~ x1 + x2", data=df).fit()

model.summary()        # the classic stats table
model.params           # coefficients (pandas Series)
model.pvalues          # per-coefficient p-values
model.rsquared         # R²
model.conf_int()       # 95% CI per coefficient
model.resid            # residuals
model.predict(df)      # predictions (formula) / fitted
model.fittedvalues

Model families

Model	API	Use
OLS	sm.OLS(y, X).fit() / smf.ols("y ~ x", data)	Continuous outcome
Logit	sm.Logit(y, X).fit() / smf.logit("y ~ x", data)	Binary outcome
Probit	sm.Probit(...) / smf.probit(...)	Binary outcome (normal CDF)
Poisson	sm.Poisson(...) / smf.poisson(...)	Counts
GLM	sm.GLM(y, X, family=sm.families.Binomial()).fit()	Generalized linear models

logit = smf.logit("churned ~ age + plan", data=df).fit()
logit.summary()
odds = np.exp(logit.params)       # odds ratios
proba = logit.predict(df)         # predicted probabilities

Formula syntax highlights

Formula	Meaning
y ~ x1 + x2	Main effects
y ~ x1 * x2	Both mains + interaction x1:x2
y ~ x1 + x2 - 1	No intercept
y ~ np.log(x)	Transform inline
y ~ C(group)	Treat as categorical (dummy-coded)
y ~ C(group, Treatment(reference='b'))	Choose reference level
y ~ x1 + I(x1**2)	Raw polynomial term

Common Operations

ANOVA

import statsmodels.api as sm
from statsmodels.formula.api import ols

m = ols("score ~ C(group)", data=df).fit()
table = sm.stats.anova_lm(m, typ=2)

Model diagnostics

sm.stats.durbin_watson(model.resid)         # autocorrelation (≈2 is good)
sm.stats.jarque_bera(model.resid)           # normality of residuals
sm.stats.het_breuschpagan(model.resid, model.model.exog)   # heteroscedasticity
sm.qqplot(model.resid, line="s")            # QQ plot (use plt.show())

Statistical helpers

sm.stats.ttest_ind(a, b)                    # t, p, df
sm.stats.proportions_ztest(count, nobs)     # proportion test
sm.stats.multipletests(pvals, method="fdr_bh")   # multiple-comparison correction
sm.stats.DescrStatsW(x).tconfint_mean()     # CI of the mean

Time series

import statsmodels.api as sm

# Stationarity: ADF test
adf, p, *_ = sm.tsa.adfuller(series)
# p < 0.05 -> stationary (no unit root)

# ACF / PACF for order selection
sm.graphics.tsa.plot_acf(series, lags=40)
sm.graphics.tsa.plot_pacf(series, lags=40)

# Seasonal decomposition
decomp = sm.tsa.seasonal_decompose(series, model="additive", period=12)
decomp.trend / decomp.seasonal / decomp.resid

# ARIMA / SARIMA
model = sm.tsa.ARIMA(series, order=(1, 1, 1)).fit()
forecast = model.forecast(steps=12)

# Exponential smoothing (good baseline)
from statsmodels.tsa.holtwinters import ExponentialSmoothing
fit = ExponentialSmoothing(series, trend="add", seasonal="add", seasonal_periods=12).fit()
pred = fit.forecast(12)

ARIMA (state-space) is the modern estimator; the older SARIMAX supports exogenous regressors: sm.tsa.SARIMAX(y, exog=X, order=(p,d,q)).

Practical Examples

1. Regression with interpretation

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

df = sm.datasets.get_rdataset("mtcars").data

m = smf.ols("mpg ~ wt + hp + C(cyl)", data=df).fit()
print(m.summary())

# Interpretation
coef = m.params["wt"]                 # change in mpg per +1 unit weight
ci = m.conf_int().loc["wt"]
print(f"wt effect: {coef:.2f} mpg per unit (95% CI {ci[0]:.2f}, {ci[1]:.2f})")

2. Logistic regression for churn

df = sm.datasets.get_rdataset("iris").data
df["is_setosa"] = (df["Species"] == "setosa").astype(int)

m = smf.logit("is_setosa ~ sepal_length + petal_length", data=df).fit()
print(m.summary())
df["p_setosa"] = m.predict(df)

3. ARIMA forecast of monthly data

import pandas as pd
from statsmodels.tsa.arima.model import ARIMA

idx = pd.date_range("2020-01", periods=60, freq="ME")
series = pd.Series(np.sin(np.arange(60)/4)*5 + 100, index=idx)

fit = ARIMA(series, order=(1, 1, 1)).fit()
fc = fit.forecast(steps=12)
print(fit.summary())                    # AIC, coefficients, Ljung-Box test

Common Patterns

# Train/test honest model comparison
train_mask = df["date"] < "2024-01-01"
m = smf.ols("y ~ x1 + x2", data=df[train_mask]).fit()
mse = ((df[~train_mask]["y"] - m.predict(df[~train_mask])) ** 2).mean()

# Clustered standard errors (panel / repeated measures)
m.get_robustcov_results(cov_type="cluster", clusters=df["group"])

# Compare nested models with an F-test
small = smf.ols("y ~ x1", data=df).fit()
large = smf.ols("y ~ x1 + x2 + x3", data=df).fit()
sm.stats.anova_lm(small, large)          # significant? -> larger model wins

Best Practices

• Check model.summary() fields: coefficient signs, p-values, R², F-statistic, and the residual diagnostics at the bottom.
• sm.add_constant(X) is required in the array API — forgetting it silently fits through the origin.
• Formula API is safer with DataFrames — it handles intercepts, dummies, and interactions for you.
• Test assumptions (normality/heteroscedasticity of residuals) before trusting p-values.
• Split data before fitting when comparing models — in-sample R² overstates fit.
• Logit coefficients are log-odds — exponentiate for odds ratios.
• For forecasting, fit on train, evaluate on a holdout tail, and compare against a naive/seasonal baseline.

Common Mistakes

Mistake	Consequence	Correct approach
Forgetting add_constant (array API)	Intercept forced through origin	sm.add_constant(X)
Reading summary() p-values on non-normal residuals	Invalid inference	Run jarque_bera/het_breuschpagan first
y ~ x1 string when using sm.OLS	TypeError (expects arrays)	Use smf.ols for formulas
Reporting in-sample R² as model quality	Overconfidence	Holdout / CV evaluation
Fitting ARIMA without checking stationarity	Garbage forecasts	adfuller; difference once (d=1)
Forgetting the target for prediction on new data	KeyError/length mismatch	Ensure same column names & order as training

Troubleshooting

Symptom	Cause	Fix
ValueError: endog has evaluated to an array with multiple columns	y given as 2-D / DataFrame	Pass a Series: df["y"]
PatsyError	Formula typo or unknown name	Check column names and operators
PerfectSeparationError (Logit)	Classes perfectly separable	Regularize or remove a feature
NaN in output	Missing values in inputs	df.dropna(subset=[...]) first
LinAlgWarning: matrix is singular	Collinear/duplicate columns	Drop correlated features
Forecast returns NaN	Non-stationary or too-short series	Difference, more data, simpler model

Related Cheat Sheets

• SciPy — scipy.stats hypothesis tests feed model diagnostics
• pandas.md — prepare and clean data before modeling
• NumPy — arrays under the hood
• Machine Learning — predictive modeling complements statistical inference
• Mathematics & Statistics — theory behind regression
• Home · Central Navigation Index

Official Documentation

• statsmodels documentation
• statsmodels: regression examples
• patsy formula language
• Time series analysis

Last Updated

2026-09-06 — Initial version.