statsmodels-statistical-modeling

statsmodels

概述

Statsmodels 为 Python 提供经典统计建模和严谨的推断。它涵盖线性模型、广义线性模型、离散选择、时间序列和全面的诊断工具。与以预测为重点的 scikit-learn 不同，statsmodels 强调系数解读、p 值、置信区间和模型诊断。

适用场景

• 拟合线性回归（OLS、WLS、GLS），并查看详细的系数表和诊断结果
• 进行逻辑回归，获取临床或流行病学研究所需的优势比和边际效应
• 使用 Poisson 或负二项回归分析计数数据
• 使用 ARIMA、SARIMAX 或指数平滑进行时间序列预测
• 执行 ANOVA、t 检验或非参数检验，并进行适当的校正
• 检验模型假设（异方差、自相关、残差正态性）
• 使用 AIC/BIC 或似然比检验比较模型
• 使用类似 R 的公式接口（y ~ x1 + x2 + C(group)）直观地指定模型
• 若要进行以预测为重点的机器学习，包括交叉验证和超参数调优，请改用 scikit-learn
• 若要进行基于贝叶斯的建模和后验推断，请改用 pymc

前提条件

• Python 包：statsmodels、numpy、pandas、scipy
• 可选：matplotlib（用于诊断图）、patsy（用于公式 API，随 statsmodels 一同提供）
• 数据：pandas DataFrame 或 NumPy 数组形式的表格数据

pip install statsmodels numpy pandas matplotlib

快速入门

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

# Generate sample data
np.random.seed(42)
n = 100
df = pd.DataFrame({
    "x1": np.random.randn(n),
    "x2": np.random.randn(n),
    "group": np.random.choice(["A", "B"], n)
})
df["y"] = 2 + 3 * df["x1"] - 1.5 * df["x2"] + np.random.randn(n)

# OLS with formula API (R-style)
results = smf.ols("y ~ x1 + x2 + C(group)", data=df).fit()
print(results.summary())
print(f"R²: {results.rsquared:.3f}, AIC: {results.aic:.1f}")

核心 API

模块 1：线性回归（OLS、WLS、GLS）

带有全面诊断功能的标准线性模型。

import statsmodels.api as sm
import numpy as np

# Generate data
np.random.seed(42)
X = np.random.randn(200, 3)
y = 1 + 2*X[:, 0] - 0.5*X[:, 1] + np.random.randn(200)

# ALWAYS add constant for intercept
X_const = sm.add_constant(X)
results = sm.OLS(y, X_const).fit()

print(results.summary())
print(f"\nCoefficients: {results.params}")
print(f"P-values: {results.pvalues}")
print(f"R²: {results.rsquared:.4f}")

# Predictions with confidence intervals
pred = results.get_prediction(X_const[:5])
print(pred.summary_frame())

# Robust standard errors (heteroskedasticity-consistent)
results_robust = sm.OLS(y, X_const).fit(cov_type="HC3")
print("Robust SEs:", results_robust.bse)

# Weighted Least Squares
weights = 1 / np.abs(results.resid + 0.1)  # Example weights
results_wls = sm.WLS(y, X_const, weights=weights).fit()
print(f"WLS R²: {results_wls.rsquared:.4f}")

模块 2：广义线性模型（GLM）

将回归扩展到非正态结果（二元、计数、连续正值）。

import statsmodels.api as sm
import numpy as np

# Poisson regression for count data
np.random.seed(42)
X = np.random.randn(200, 2)
X_const = sm.add_constant(X)
y_counts = np.random.poisson(np.exp(0.5 + 0.3*X[:, 0]))

model = sm.GLM(y_counts, X_const, family=sm.families.Poisson())
results = model.fit()
print(results.summary())

# Rate ratios
rate_ratios = np.exp(results.params)
print(f"Rate ratios: {rate_ratios}")

# Check overdispersion
overdispersion = results.pearson_chi2 / results.df_resid
print(f"Overdispersion ratio: {overdispersion:.2f}")
if overdispersion > 1.5:
    print("→ Consider Negative Binomial model")

模块 3：离散选择模型（Logit、Probit、计数）

二元、多项和计数结果模型。

import statsmodels.api as sm
import numpy as np

# Logistic regression
np.random.seed(42)
X = np.random.randn(300, 2)
X_const = sm.add_constant(X)
prob = 1 / (1 + np.exp(-(0.5 + X[:, 0] - 0.5*X[:, 1])))
y_binary = np.random.binomial(1, prob)

logit_results = sm.Logit(y_binary, X_const).fit()
print(logit_results.summary())

# Odds ratios
odds_ratios = np.exp(logit_results.params)
print(f"Odds ratios: {odds_ratios}")

# Marginal effects (at means)
margeff = logit_results.get_margeff()
print(margeff.summary())

# Predicted probabilities
probs = logit_results.predict(X_const[:5])
print(f"Predicted P(Y=1): {probs}")

模块 4：时间序列（ARIMA、SARIMAX）

单变量和多变量时间序列建模与预测。

import statsmodels.api as sm
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.stattools import adfuller
import numpy as np
import pandas as pd

# Generate time series
np.random.seed(42)
dates = pd.date_range("2020-01-01", periods=200, freq="D")
y = np.cumsum(np.random.randn(200)) + 50
ts = pd.Series(y, index=dates)

# Stationarity test
adf_result = adfuller(ts)
print(f"ADF statistic: {adf_result[0]:.4f}, p-value: {adf_result[1]:.4f}")
print("Stationary" if adf_result[1] < 0.05 else "Non-stationary → difference")

# Fit ARIMA
model = ARIMA(ts, order=(1, 1, 1))
results = model.fit()
print(results.summary())

# Forecast with confidence intervals
forecast = results.get_forecast(steps=30)
forecast_df = forecast.summary_frame()
print(f"30-day forecast:\n{forecast_df.head()}")

# Seasonal ARIMA (SARIMAX)
from statsmodels.tsa.statespace.sarimax import SARIMAX

# Monthly data with yearly seasonality
model_sarima = SARIMAX(ts, order=(1, 1, 1), seasonal_order=(1, 1, 1, 12))
results_sarima = model_sarima.fit(disp=False)
print(f"AIC: {results_sarima.aic:.1f}")

# Diagnostic plots
results_sarima.plot_diagnostics(figsize=(12, 8))

模块 5：统计检验与诊断

假设检验、假设验证和模型验证。

import statsmodels.api as sm
from statsmodels.stats.diagnostic import het_breuschpagan, acorr_ljungbox
from statsmodels.stats.stattools import jarque_bera
import numpy as np

# Fit a model first
np.random.seed(42)
X = sm.add_constant(np.random.randn(200, 2))
y = 1 + 2*X[:, 1] + np.random.randn(200) * X[:, 1]  # Heteroskedastic
results = sm.OLS(y, X).fit()

# Heteroskedasticity test (Breusch-Pagan)
bp_stat, bp_p, _, _ = het_breuschpagan(results.resid, X)
print(f"Breusch-Pagan p-value: {bp_p:.4f} {'→ heteroskedastic' if bp_p < 0.05 else '→ OK'}")

# Normality test (Jarque-Bera)
jb_stat, jb_p, _, _ = jarque_bera(results.resid)
print(f"Jarque-Bera p-value: {jb_p:.4f} {'→ non-normal' if jb_p < 0.05 else '→ OK'}")

# Autocorrelation test (Ljung-Box)
lb_result = acorr_ljungbox(results.resid, lags=[10], return_df=True)
print(f"Ljung-Box p-value (lag 10): {lb_result['lb_pvalue'].values[0]:.4f}")

# Variance Inflation Factor (multicollinearity)
from statsmodels.stats.outliers_influence import variance_inflation_factor

vif_data = pd.DataFrame({
    "Variable": [f"x{i}" for i in range(X.shape[1])],
    "VIF": [variance_inflation_factor(X, i) for i in range(X.shape[1])]
})
print(vif_data)  # VIF > 10 suggests multicollinearity

模块 6：公式 API（类似 R）

使用公式直观地指定模型，并自动进行虚拟变量编码。

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

np.random.seed(42)
df = pd.DataFrame({
    "y": np.random.randn(100),
    "x1": np.random.randn(100),
    "x2": np.random.randn(100),
    "group": np.random.choice(["A", "B", "C"], 100),
})

# Formula with categoricals (auto dummy-coded)
res = smf.ols("y ~ x1 + x2 + C(group)", data=df).fit()
print(res.summary())

# Interactions
res2 = smf.ols("y ~ x1 * x2", data=df).fit()  # x1 + x2 + x1:x2
print(f"Interaction term p-value: {res2.pvalues['x1:x2']:.4f}")

# Logit via formula
df["binary"] = (df["y"] > 0).astype(int)
logit_res = smf.logit("binary ~ x1 + x2 + C(group)", data=df).fit()
print(f"Logit AIC: {logit_res.aic:.1f}")

常见工作流

工作流 1：完整的回归分析

目标：拟合 OLS、验证假设，并在需要时使用稳健标准误（SE）。

import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.stats.outliers_influence import variance_inflation_factor
import numpy as np
import pandas as pd

# 1. Fit initial model
np.random.seed(42)
df = pd.DataFrame({"y": np.random.randn(200), "x1": np.random.randn(200), "x2": np.random.randn(200)})
df["y"] = 2 + 3*df["x1"] - df["x2"] + np.random.randn(200)

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

# 2. Check heteroskedasticity
bp_stat, bp_p, _, _ = het_breuschpagan(results.resid, results.model.exog)
print(f"Breusch-Pagan p: {bp_p:.4f}")

# 3. If heteroskedastic, use robust SEs
if bp_p < 0.05:
    results = smf.ols("y ~ x1 + x2", data=df).fit(cov_type="HC3")
    print("Using HC3 robust standard errors")

# 4. Check multicollinearity
X = results.model.exog
for i in range(1, X.shape[1]):  # skip constant
    print(f"VIF x{i}: {variance_inflation_factor(X, i):.2f}")

# 5. Final results
print(results.summary())
print(f"\nAIC: {results.aic:.1f}, BIC: {results.bic:.1f}")

工作流 2：模型比较

目标：使用适当的标准比较嵌套和非嵌套模型。

import statsmodels.formula.api as smf
from scipy import stats
import pandas as pd
import numpy as np

np.random.seed(42)
df = pd.DataFrame({"y": np.random.randn(200), "x1": np.random.randn(200),
                    "x2": np.random.randn(200), "x3": np.random.randn(200)})
df["y"] = 1 + 2*df["x1"] - df["x2"] + 0.1*df["x3"] + np.random.randn(200)

# Fit nested models
m1 = smf.ols("y ~ x1", data=df).fit()
m2 = smf.ols("y ~ x1 + x2", data=df).fit()
m3 = smf.ols("y ~ x1 + x2 + x3", data=df).fit()

# Compare via AIC/BIC (lower = better)
comparison = pd.DataFrame({
    "R²": [m.rsquared for m in [m1, m2, m3]],
    "AIC": [m.aic for m in [m1, m2, m3]],
    "BIC": [m.bic for m in [m1, m2, m3]],
}, index=["y~x1", "y~x1+x2", "y~x1+x2+x3"])
print(comparison)

# Likelihood ratio test (nested: m2 vs m3)
lr_stat = 2 * (m3.llf - m2.llf)
p_val = 1 - stats.chi2.cdf(lr_stat, df=m3.df_model - m2.df_model)
print(f"\nLR test (m3 vs m2): stat={lr_stat:.2f}, p={p_val:.4f}")

工作流 3：时间序列预测流程

目标：检验平稳性、确定模型阶数并进行预测。

from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.stattools import adfuller
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# Generate data
np.random.seed(42)
ts = pd.Series(np.cumsum(np.random.randn(200)) + 100,
               index=pd.date_range("2020-01-01", periods=200, freq="D"))

# 1. Test stationarity
adf_p = adfuller(ts)[1]
print(f"ADF p-value: {adf_p:.4f} → {'stationary' if adf_p < 0.05 else 'non-stationary'}")

# 2. Identify order from ACF/PACF (on differenced series)
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 6))
plot_acf(ts.diff().dropna(), lags=20, ax=ax1)
plot_pacf(ts.diff().dropna(), lags=20, ax=ax2)
plt.savefig("acf_pacf.png", dpi=150, bbox_inches="tight")

# 3. Fit and forecast
model = ARIMA(ts[:180], order=(1, 1, 1))
results = model.fit()
forecast = results.get_forecast(steps=20)
fc_df = forecast.summary_frame()
print(f"ARIMA AIC: {results.aic:.1f}")
print(f"Forecast (first 5 days):\n{fc_df.head()}")

关键参数

参数	模块	默认值	范围 / 选项	作用
cov_type	OLS/WLS/GLM	"nonrobust"	"HC0"-"HC3"、"HAC"、"cluster"	稳健协方差估计量
family	GLM	必填	Poisson()、Binomial()、Gamma() 等	分布族
order	ARIMA	必填	(p, d, q) 元组	AR 阶数、差分阶数、MA 阶数
seasonal_order	SARIMAX	(0,0,0,0)	(P, D, Q, s) 元组	季节性 ARIMA 参数
alpha	summary()、conf_int()	0.05	0.01-0.10	置信区间的显著性水平
maxiter	所有 .fit()	35-100	50-1000	最大优化迭代次数
method	.fit()	取决于模型	"newton"、"bfgs"、"lbfgs"、"powell"	优化算法
lags	ACF/PACF	None	10-50	要显示的滞后阶数

最佳实践

1. OLS/GLM 始终添加常数项：使用 sm.add_constant(X)，或使用公式 API（会自动添加截距）

2. 根据结果类型选择模型：二元结果 → Logit/Probit，计数结果 → Poisson/NegBin，连续结果 → OLS/WLS，时间序列 → ARIMA

3. 解读结果前先检查诊断指标：对残差进行 Breusch-Pagan（异方差）、Jarque-Bera（正态性）和 Ljung-Box（自相关）检验

4. 假设不成立时使用稳健标准误：对于异方差稳健推断，使用 results = model.fit(cov_type="HC3")

5. 报告效应量，而不只是 p 值：除显著性检验外，还应包含系数、置信区间和 R²

6. 探索性分析优先使用公式 API：smf.ols("y ~ x1 * x2 + C(group)", data=df) 可读性更好，且会自动处理分类变量

7. 时间序列建模前先检验平稳性：使用 ADF 检验；如果序列非平稳，则进行差分

常用方法

方法：ANOVA 与事后检验

适用场景：比较 3 个或更多组的均值。

import statsmodels.formula.api as smf
from statsmodels.stats.multicomp import pairwise_tukeyhsd
import pandas as pd
import numpy as np

np.random.seed(42)
df = pd.DataFrame({"value": np.concatenate([np.random.normal(m, 1, 30) for m in [5, 6, 7]]),
                    "group": np.repeat(["A", "B", "C"], 30)})

# One-way ANOVA
anova = smf.ols("value ~ C(group)", data=df).fit()
print(sm.stats.anova_lm(anova))

# Post-hoc Tukey HSD
tukey = pairwise_tukeyhsd(df["value"], df["group"], alpha=0.05)
print(tukey)

方法：样本量的功效分析

适用场景：在研究开始前确定所需样本量。

from statsmodels.stats.power import TTestIndPower

analysis = TTestIndPower()
# What sample size for medium effect (d=0.5), 80% power, alpha=0.05?
n = analysis.solve_power(effect_size=0.5, alpha=0.05, power=0.8)
print(f"Required n per group: {n:.0f}")

# Power for given sample size
power = analysis.solve_power(effect_size=0.5, alpha=0.05, nobs1=50)
print(f"Power with n=50: {power:.3f}")

方法：混合效应模型

适用场景：分层/聚类数据（医院中的患者、学校中的学生）。

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

np.random.seed(42)
df = pd.DataFrame({
    "y": np.random.randn(100), "x": np.random.randn(100),
    "group": np.repeat(range(10), 10)
})

# Random intercept model
model = smf.mixedlm("y ~ x", data=df, groups=df["group"])
results = model.fit()
print(results.summary())

故障排除

问题	原因	解决方案
MissingDataError	数据中存在 NaN 值	删除缺失值：df.dropna()，或在拟合前进行插补
结果中没有截距	忘记使用 sm.add_constant()	始终添加常数项，或使用公式 API smf.ols()
ConvergenceWarning	优化失败	增大 maxiter，尝试不同的 method，或缩放变量
Poisson 模型过度离散	方差 > 均值	切换到 NegativeBinomial，或使用 GLM(family=NegativeBinomial())
时间序列非平稳	存在趋势或单位根	对序列进行差分（ts.diff()），或增大 ARIMA 中的 d
矩阵奇异错误	完全多重共线性	移除冗余变量；检查 VIF 是否 > 10
与 R 的结果不同	默认设置不同	检查：常数项、链接函数、优化器、标准误类型
Logit 中出现 PerfectSeparationError	预测变量可完美区分类别	使用正则化逻辑回归（惩罚 MLE）或 Firth 方法

参考资料

• statsmodels 文档 — 官方文档
• statsmodels 用户指南 — 教程
• statsmodels API 参考 — 完整 API
• Seabold, S. & Perktold, J. (2010). Statsmodels：使用 Python 进行计量经济学与统计建模。第 9 届 Python 科学大会论文集。