statsmodels

Statsmodels：统计建模与计量经济学

概述

Statsmodels 是 Python 中首屈一指的统计建模库，提供涵盖广泛统计方法的估计、推断和诊断工具。使用此技能进行严谨的统计分析，从简单的线性回归到复杂的时间序列模型和计量经济学分析。

何时使用此技能

此技能应在以下情况下使用：

• 拟合回归模型（OLS、WLS、GLS、分位数回归）
• 执行广义线性建模（逻辑回归、泊松、Gamma 等）
• 分析离散结果（二分类、多分类、计数、有序）
• 进行时间序列分析（ARIMA、SARIMAX、VAR、预测）
• 运行统计检验和诊断
• 检验模型假设（异方差性、自相关、正态性）
• 检测异常值和有影响的观测
• 比较模型（AIC/BIC、似然比检验）
• 估计因果效应
• 生成可发表级别的统计表格和推断

快速入门指南

线性回归（OLS）

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

# Prepare data - ALWAYS add constant for intercept
X = sm.add_constant(X_data)

# Fit OLS model
model = sm.OLS(y, X)
results = model.fit()

# View comprehensive results
print(results.summary())

# Key results
print(f"R-squared: {results.rsquared:.4f}")
print(f"Coefficients:\\n{results.params}")
print(f"P-values:\\n{results.pvalues}")

# Predictions with confidence intervals
predictions = results.get_prediction(X_new)
pred_summary = predictions.summary_frame()
print(pred_summary)  # includes mean, CI, prediction intervals

# Diagnostics
from statsmodels.stats.diagnostic import het_breuschpagan
bp_test = het_breuschpagan(results.resid, X)
print(f"Breusch-Pagan p-value: {bp_test[1]:.4f}")

# Visualize residuals
import matplotlib.pyplot as plt
plt.scatter(results.fittedvalues, results.resid)
plt.axhline(y=0, color='r', linestyle='--')
plt.xlabel('Fitted values')
plt.ylabel('Residuals')
plt.show()

逻辑回归（二分类结果）

from statsmodels.discrete.discrete_model import Logit

# Add constant
X = sm.add_constant(X_data)

# Fit logit model
model = Logit(y_binary, X)
results = model.fit()

print(results.summary())

# Odds ratios
odds_ratios = np.exp(results.params)
print("Odds ratios:\\n", odds_ratios)

# Predicted probabilities
probs = results.predict(X)

# Binary predictions (0.5 threshold)
predictions = (probs > 0.5).astype(int)

# Model evaluation
from sklearn.metrics import classification_report, roc_auc_score

print(classification_report(y_binary, predictions))
print(f"AUC: {roc_auc_score(y_binary, probs):.4f}")

# Marginal effects
marginal = results.get_margeff()
print(marginal.summary())

时间序列（ARIMA）

from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# Check stationarity
from statsmodels.tsa.stattools import adfuller

adf_result = adfuller(y_series)
print(f"ADF p-value: {adf_result[1]:.4f}")

if adf_result[1] > 0.05:
    # Series is non-stationary, difference it
    y_diff = y_series.diff().dropna()

# Plot ACF/PACF to identify p, q
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
plot_acf(y_diff, lags=40, ax=ax1)
plot_pacf(y_diff, lags=40, ax=ax2)
plt.show()

# Fit ARIMA(p,d,q)
model = ARIMA(y_series, order=(1, 1, 1))
results = model.fit()

print(results.summary())

# Forecast
forecast = results.forecast(steps=10)
forecast_obj = results.get_forecast(steps=10)
forecast_df = forecast_obj.summary_frame()

print(forecast_df)  # includes mean and confidence intervals

# Residual diagnostics
results.plot_diagnostics(figsize=(12, 8))
plt.show()

广义线性模型（GLM）

import statsmodels.api as sm

# Poisson regression for count data
X = sm.add_constant(X_data)
model = sm.GLM(y_counts, X, family=sm.families.Poisson())
results = model.fit()

print(results.summary())

# Rate ratios (for Poisson with log link)
rate_ratios = np.exp(results.params)
print("Rate ratios:\\n", rate_ratios)

# Check overdispersion
overdispersion = results.pearson_chi2 / results.df_resid
print(f"Overdispersion: {overdispersion:.2f}")

if overdispersion > 1.5:
    # Use Negative Binomial instead
    from statsmodels.discrete.count_model import NegativeBinomial
    nb_model = NegativeBinomial(y_counts, X)
    nb_results = nb_model.fit()
    print(nb_results.summary())

核心统计建模能力

1. 线性回归模型

针对连续结果、具有各种误差结构的全面线性模型套件。

可用模型：

• OLS：误差为独立同分布的标准线性回归
• WLS：用于异方差误差的加权最小二乘法
• GLS：用于任意协方差结构的广义最小二乘法
• GLSAR：用于时间序列的具有自回归误差的 GLS
• 分位数回归：条件分位数（对异常值稳健）
• 混合效应：具有随机效应的分层/多水平模型
• 递归/滚动：时变参数估计

关键特性：

• 全面的诊断检验
• 稳健标准误（HC、HAC、聚类稳健）
• 影响统计量（Cook 距离、杠杆值、DFFITS）
• 假设检验（F 检验、Wald 检验）
• 模型比较（AIC、BIC、似然比检验）
• 带置信区间和预测区间的预测

何时使用： 连续结果变量、需要系数推断、需要诊断

参考： 关于模型选择、诊断和最佳实践的详细指导，请参见 references/linear_models.md。

2. 广义线性模型（GLM）

将线性模型扩展到非正态分布的灵活框架。

分布族：

• 二项分布：二分类结果或比例（逻辑回归）
• 泊松分布：计数数据
• 负二项分布：过度离散的计数
• Gamma 分布：正值连续、右偏数据
• 逆高斯分布：具有特定方差结构的正值连续数据
• 高斯分布：等价于 OLS
• Tweedie 分布：用于半连续数据的灵活分布族

连接函数：

• Logit、Probit、Log、Identity、Inverse、Sqrt、CLogLog、Power
• 根据解释需求和模型拟合进行选择

关键特性：

• 通过 IRLS 进行最大似然估计
• 偏差残差和 Pearson 残差
• 拟合优度统计量
• 伪 R 方度量
• 稳健标准误

何时使用： 非正态结果、需要灵活的方差和连接函数设定

参考： 关于分布族选择、连接函数、解释和诊断，请参见 references/glm.md。

3. 离散选择模型

用于分类和计数结果的模型。

二值模型：

• Logit：逻辑回归（优势比）
• Probit：Probit 回归（正态分布）

多类别模型：

• MNLogit：无序类别（3 个及以上水平）
• 条件 Logit：含备选项特定变量的选择模型
• 有序模型：有序结果（有序类别）

计数模型：

• Poisson：标准计数模型
• 负二项：过度分散计数
• 零膨胀：过多零值（ZIP、ZINB）
• 栅栏模型：针对零值过多数据的两阶段模型

主要特性：

• 最大似然估计
• 均值处边际效应或平均边际效应
• 通过 AIC/BIC 进行模型比较
• 预测概率与分类
• 拟合优度检验

适用场景： 二值、分类或计数结果

参考： 模型选择、解释和评估参见 references/discrete_choice.md。

4. 时间序列分析

全面的时间序列建模与预测能力。

单变量模型：

• AutoReg (AR)：自回归模型
• ARIMA：自回归积分滑动平均
• SARIMAX：带外生变量的季节性 ARIMA
• 指数平滑：简单、Holt、Holt-Winters
• ETS： innovations 状态空间模型

多变量模型：

• VAR：向量自回归
• VARMAX：带 MA 和外生变量的 VAR
• 动态因子模型：提取共同因子
• VECM：向量误差修正模型（协整）

高级模型：

• 状态空间：卡尔曼滤波、自定义设定
• 区制转换：马尔可夫转换模型
• ARDL：自回归分布滞后

主要特性：

• 用于模型识别的 ACF/PACF 分析
• 平稳性检验（ADF、KPSS）
• 带预测区间的预测
• 残差诊断（Ljung-Box、异方差）
• Granger 因果检验
• 脉冲响应函数（IRF）
• 预测误差方差分解（FEVD）

适用场景： 按时间排序的数据、预测、理解时间动态

参考： 模型选择、诊断和预测方法参见 references/time_series.md。

5. 统计检验与诊断

用于模型验证的广泛检验与诊断能力。

残差诊断：

• 自相关检验（Ljung-Box、Durbin-Watson、Breusch-Godfrey）
• 异方差检验（Breusch-Pagan、White、ARCH）
• 正态性检验（Jarque-Bera、Omnibus、Anderson-Darling、Lilliefors）
• 设定检验（RESET、Harvey-Collier）

影响与离群值：

• 杠杆值（hat 值）
• Cook 距离
• DFFITS 和 DFBETAs
• 学生化残差
• 影响图

假设检验：

• t 检验（单样本、双样本、配对）
• 比例检验
• 卡方检验
• 非参数检验（Mann-Whitney、Wilcoxon、Kruskal-Wallis）
• 方差分析（单因素、双因素、重复测量）

多重比较：

• Tukey HSD
• Bonferroni 校正
• 错误发现率（FDR）

效应量与检验力：

• Cohen's d、eta 平方
• t 检验、比例检验的检验力分析
• 样本量计算

稳健推断：

• 异方差一致标准误（HC0-HC3）
• HAC 标准误（Newey-West）
• 聚类稳健标准误

适用场景： 验证假设、发现问题、确保稳健推断

参考： 全面的检验与诊断流程参见 references/stats_diagnostics.md。

公式 API（R 风格）

Statsmodels 支持 R 风格公式，以便直观地指定模型：

import statsmodels.formula.api as smf

# OLS with formula
results = smf.ols('y ~ x1 + x2 + x1:x2', data=df).fit()

# Categorical variables (automatic dummy coding)
results = smf.ols('y ~ x1 + C(category)', data=df).fit()

# Interactions
results = smf.ols('y ~ x1 * x2', data=df).fit()  # x1 + x2 + x1:x2

# Polynomial terms
results = smf.ols('y ~ x + I(x**2)', data=df).fit()

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

# Poisson
results = smf.poisson('count ~ x1 + x2', data=df).fit()

# ARIMA (not available via formula, use regular API)

模型选择与比较

信息准则

# Compare models using AIC/BIC
models = {
    'Model 1': model1_results,
    'Model 2': model2_results,
    'Model 3': model3_results
}

comparison = pd.DataFrame({
    'AIC': {name: res.aic for name, res in models.items()},
    'BIC': {name: res.bic for name, res in models.items()},
    'Log-Likelihood': {name: res.llf for name, res in models.items()}
})

print(comparison.sort_values('AIC'))
# Lower AIC/BIC indicates better model

似然比检验（嵌套模型）

# For nested models (one is subset of the other)
from scipy import stats

lr_stat = 2 * (full_model.llf - reduced_model.llf)
df = full_model.df_model - reduced_model.df_model
p_value = 1 - stats.chi2.cdf(lr_stat, df)

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

if p_value < 0.05:
    print("Full model significantly better")
else:
    print("Reduced model preferred (parsimony)")

交叉验证

from sklearn.model_selection import KFold
from sklearn.metrics import mean_squared_error

kf = KFold(n_splits=5, shuffle=True, random_state=42)
cv_scores = []

for train_idx, val_idx in kf.split(X):
    X_train, X_val = X.iloc[train_idx], X.iloc[val_idx]
    y_train, y_val = y.iloc[train_idx], y.iloc[val_idx]

    # Fit model
    model = sm.OLS(y_train, X_train).fit()

    # Predict
    y_pred = model.predict(X_val)

    # Score
    rmse = np.sqrt(mean_squared_error(y_val, y_pred))
    cv_scores.append(rmse)

print(f"CV RMSE: {np.mean(cv_scores):.4f} ± {np.std(cv_scores):.4f}")

最佳实践

数据准备

1. 始终添加常数项：除非排除截距，否则使用 sm.add_constant()
2. 检查缺失值：在拟合前处理或插补
3. 按需缩放：有助于收敛和解释（但树模型不需要）
4. 编码分类变量：使用公式 API 或手动哑变量编码

模型构建

1. 从简单开始：先使用基础模型，按需增加复杂度
2. 检查假设：检验残差、异方差性、自相关
3. 使用合适的模型：使模型匹配结果类型（二值→Logit，计数→Poisson）
4. 考虑替代方案：若假设被违反，使用稳健方法或不同模型

推断

1. 报告效应量：而不仅仅是 p 值
2. 使用稳健标准误：存在异方差性或聚类时
3. 多重比较：检验多个假设时进行校正
4. 置信区间：始终与点估计一起报告

模型评估

1. 检查残差：绘制残差与拟合值图、Q-Q 图
2. 影响力诊断：识别并调查有影响力的观测
3. 样本外验证：在留出集上测试或交叉验证
4. 比较模型：非嵌套模型使用 AIC/BIC，嵌套模型使用 LR 检验

报告

1. 全面总结：使用 .summary() 获取详细输出
2. 记录决策：说明变换、排除的观测
3. 谨慎解释：考虑连接函数（例如对数连接使用 exp(β)）
4. 可视化：绘制预测值、置信区间、诊断图

常见工作流

工作流 1：线性回归分析

1. 探索数据（图表、描述性统计）
2. 拟合初始 OLS 模型
3. 检查残差诊断
4. 检验异方差性、自相关
5. 检查多重共线性（VIF）
6. 识别有影响力的观测
7. 如需要，使用稳健标准误重新拟合
8. 解释系数和推断
9. 在留出集上或通过 CV 验证

工作流 2：二分类

1. 拟合逻辑回归（Logit）
2. 检查收敛问题
3. 解释优势比
4. 计算边际效应
5. 评估分类性能（AUC、混淆矩阵）
6. 检查有影响力的观测
7. 与替代模型（Probit）比较
8. 在测试集上验证预测

工作流 3：计数数据分析

1. 拟合 Poisson 回归
2. 检查过度离散
3. 若过度离散，拟合负二项模型
4. 检查过多零值（考虑 ZIP/ZINB）
5. 解释率比
6. 评估拟合优度
7. 通过 AIC 比较模型
8. 验证预测

工作流 4：时间序列预测

1. 绘制序列图，检查趋势/季节性
2. 检验平稳性（ADF、KPSS）
3. 若非平稳则差分
4. 从 ACF/PACF 识别 p、q
5. 拟合 ARIMA 或 SARIMAX
6. 检查残差诊断（Ljung-Box）
7. 生成带置信区间的预测
8. 在测试集上评估预测精度

参考文档

此技能包含全面的参考文件以提供详细指导：

references/linear_models.md

详细涵盖线性回归模型，包括：

• OLS、WLS、GLS、GLSAR、分位数回归
• 混合效应模型
• 递归回归和滚动回归
• 全面诊断（异方差性、自相关、多重共线性）
• 影响力统计量和异常值检测
• 稳健标准误（HC、HAC、聚类）
• 假设检验和模型比较

references/glm.md

广义线性模型完整指南：

• 所有分布族（Binomial、Poisson、Gamma 等）
• 连接函数及各自适用场景
• 模型拟合与解释
• 伪 R 方和拟合优度
• 诊断和残差分析
• 应用（逻辑回归、Poisson 回归、Gamma 回归）

references/discrete_choice.md

离散结果模型的综合指南：

• 二值模型（Logit、Probit）
• 多项模型（MNLogit、Conditional Logit）
• 计数模型（Poisson、Negative Binomial、Zero-Inflated、Hurdle）
• 有序模型
• 边际效应与解释
• 模型诊断与比较

references/time_series.md

深入的时间序列分析指导：

• 单变量模型（AR、ARIMA、SARIMAX、Exponential Smoothing）
• 多变量模型（VAR、VARMAX、Dynamic Factor）
• 状态空间模型
• 平稳性检验与诊断
• 预测方法与评估
• Granger 因果检验、IRF、FEVD

references/stats_diagnostics.md

全面的统计检验与诊断：

• 残差诊断（自相关、异方差、正态性）
• 影响点与离群点检测
• 假设检验（参数与非参数）
• ANOVA 与事后检验
• 多重比较校正
• 稳健协方差矩阵
• 功效分析与效应量

何时参考：

• 需要详细的参数解释
• 在相似模型之间做选择
• 排查收敛或诊断问题
• 理解具体的检验统计量
• 寻找高级功能的代码示例

搜索模式：

# Find information about specific models
grep -r "Quantile Regression" references/

# Find diagnostic tests
grep -r "Breusch-Pagan" references/stats_diagnostics.md

# Find time series guidance
grep -r "SARIMAX" references/time_series.md

常见陷阱与规避方法

1. 忘记常数项：始终使用 sm.add_constant()，除非不需要截距
2. 忽略假设：检查残差、异方差、自相关
3. 结果类型对应错误的模型：二值→Logit/Probit，计数→Poisson/NB，而非 OLS
4. 未检查收敛：留意优化警告
5. 误解系数：记住链接函数（log、logit 等）
6. 在存在过度离散时使用 Poisson：检查离散度，必要时使用 Negative Binomial
7. 未使用稳健标准误：当存在异方差或聚类时
8. 过拟合：参数相对于样本量过多
9. 数据泄漏：在测试数据上拟合或使用未来信息
10. 未验证预测：始终检查样本外表现
11. 比较非嵌套模型：使用 AIC/BIC，而非 LR 检验
12. 忽略有影响的观测：检查 Cook's distance 和杠杆值
13. 多重检验：在检验大量假设时校正 p 值
14. 未对时间序列做差分：在非平稳数据上拟合 ARIMA
15. 混淆预测区间与置信区间：预测区间更宽

获取帮助

详细文档与示例：

• 官方文档：https://www.statsmodels.org/stable/
• 用户指南：https://www.statsmodels.org/stable/user-guide.html
• 示例：https://www.statsmodels.org/stable/examples/index.html
• API 参考：https://www.statsmodels.org/stable/api.html