statsmodels

Statsmodels - 统计建模与推断

Statsmodels 是 Python 与 R 风格严谨统计分析之间的桥梁。它允许用户使用公式（通过 patsy）估计模型、执行广泛的诊断检验，并生成学术出版领域的标准详细摘要表。

适用场景

• 估计具有详细诊断的线性回归模型（OLS、WLS）。
• 广义线性模型（GLM）：Logistic、Poisson、Gamma 回归。
• 时间序列分析（ARIMA、SARIMAX、VAR、状态空间模型）。
• 统计假设检验（t 检验、ANOVA、正态性、异方差性）。
• 生存分析（Kaplan-Meier、Cox 比例风险）。
• 估计处理效应并进行因果推断。
• 非参数统计（核密度估计）。

参考文档

官方文档：https://www.statsmodels.org/stable/
公式 API：https://www.statsmodels.org/stable/example_formulas.html
搜索模式：sm.OLS、smf.ols、sm.tsa、results.summary()、statsmodels.api

核心原则

Statsmodels 与 scikit-learn 的区别

特性	scikit-learn	Statsmodels
目标	预测（准确性）	推断（解释／p 值）
接口	fit / predict	fit / summary
预处理	Pipeline 对象	公式（Patsy）或设计矩阵
诊断	交叉验证	残差分析、p 值、CI

两种 API

• API（statsmodels.api）：需要显式添加常数项（截距），并使用类似 NumPy 的数组。
• 公式 API（statsmodels.formula.api）：使用 R 风格公式（y ~ x1 + x2），并可直接处理 Pandas DataFrame。（推荐大多数用户使用。）

快速参考

安装

pip install statsmodels patsy

标准导入

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

基本模式 - 普通最小二乘法（OLS）

import statsmodels.formula.api as smf

# 1. Define model with R-style formula
# 'y ~ x1 + x2' means: y = beta0 + beta1*x1 + beta2*x2
model = smf.ols('tip ~ total_bill + size', data=df_tips)

# 2. Fit the model
results = model.fit()

# 3. Print the comprehensive results table
print(results.summary())

# 4. Access specific values
p_values = results.pvalues
params = results.params # beta coefficients

关键规则

✅ 应该做

• 检查残差 - 始终绘制并检验残差（results.resid）的正态性和同方差性。
• 添加常数项 - 如果使用 sm.api（而非公式），请记得执行 X = sm.add_constant(X)，否则模型会经过原点（beta0 = 0）。
• 使用分类变量 - 在公式中使用 C() 运算符（例如 y ~ C(region)）自动创建虚拟变量。
• 指定协方差类型 - 如果怀疑存在非恒定方差（异方差性），请使用 cov_type='HC3' 或 'cluster'。
• 谨慎解读 R-squared - 如果残差呈现模式，即使 R-squared 很高，也不代表模型良好。
• 检查多重共线性 - 使用 VIF（方差膨胀因子）确保预测变量之间没有高度相关。

❌ 不要做

• 将预测等同于推断 - 模型的 R-squared 很高，并不意味着其系数代表现实世界中的因果效应。
• 忽略截距 - 大多数物理和社会过程都需要常数项。
• 使用过多预测变量导致过拟合 - 使用 AIC/BIC 指标对复杂模型进行惩罚。
• 超出数据范围进行外推 - 统计模型仅在训练数据的适用范围内有效。

反模式（绝对不要）

import statsmodels.api as sm

# ❌ BAD: Forgetting the intercept in Array API
model = sm.OLS(y, X) # This forces the line through (0,0)
results = model.fit()

# ✅ GOOD: Add the constant
X = sm.add_constant(X)
model = sm.OLS(y, X).fit()

# ❌ BAD: Ignoring categorical data strings
# smf.ols('price ~ color', data=df) # Might fail if 'color' isn't numeric

# ✅ GOOD: Explicitly tell Statsmodels it's categorical
smf.ols('price ~ C(color)', data=df).fit()

# ❌ BAD: Only looking at R-squared
# print(results.rsquared) # Tells only part of the story

# ✅ GOOD: Inspecting the whole summary
print(results.summary()) # Check p-values, F-stat, Jarque-Bera

回归分析

线性模型（OLS、WLS）

# Multiple Linear Regression with interactions
# 'x1 * x2' includes x1, x2, and the interaction term x1:x2
model = smf.ols('y ~ x1 * x2 + np.log(x3)', data=df).fit()

# Weighted Least Squares (for heteroscedastic data)
wls_model = sm.WLS(y, X, weights=1.0/variance_estimates).fit()

广义线性模型（GLM）

# Logistic Regression (Binary outcome)
logit_model = smf.logit('admit ~ gre + gpa + C(rank)', data=df).fit()

# Poisson Regression (Count data)
poisson_model = smf.poisson('num_awards ~ math + C(prog)', data=df).fit()

# Negative Binomial (For over-dispersed counts)
nb_model = smf.glm('y ~ x1', data=df, family=sm.families.NegativeBinomial()).fit()

时间序列分析（tsa）

平稳性与建模

from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.statespace.sarimax import SARIMAX

# 1. Test for Stationarity (Augmented Dickey-Fuller)
adf_result = adfuller(df['sales'])
print(f"ADF P-value: {adf_result[1]}") # p < 0.05 means stationary

# 2. SARIMAX (Seasonal ARIMA with eXogenous variables)
model = SARIMAX(df['sales'], 
                order=(1, 1, 1), 
                seasonal_order=(1, 1, 0, 12),
                exog=df['advertising'])
results = model.fit()

# 3. Forecasting
forecast = results.get_forecast(steps=12)
conf_int = forecast.conf_int()

ANOVA 与假设检验

from statsmodels.stats.anova import anova_lm

# Perform ANOVA on OLS model
model = smf.ols('yield ~ C(fertilizer) + C(soil)', data=df).fit()
anova_table = anova_lm(model, typ=2)

# Post-hoc tests (Tukey's HSD)
from statsmodels.stats.multicomp import pairwise_tukeyhsd
tukey = pairwise_tukeyhsd(endog=df['yield'], groups=df['fertilizer'], alpha=0.05)
print(tukey)

模型诊断

残差分析

import statsmodels.stats.api as sms

# 1. Test for Heteroscedasticity (Breusch-Pagan)
name = ['Lagrange multiplier statistic', 'p-value', 'f-value', 'f p-value']
test = sms.het_breuschpagan(results.resid, results.model.exog)
print(dict(zip(name, test)))

# 2. Test for Normality (Omnibus)
# Included in results.summary() by default. 

# 3. Check for Outliers (Influence)
influence = results.get_influence()
cooks_d = influence.cook_distance[0]

实用工作流

1. 稳健的科学报告流程

def analyze_experiment(df):
    """Rigorous analysis of an experimental dataset."""
    # 1. Fit model
    model = smf.ols('outcome ~ treatment + age + gender', data=df).fit()
    
    # 2. Diagnostic Plots (requires matplotlib)
    import matplotlib.pyplot as plt
    sm.graphics.plot_regress_exog(model, 'treatment')
    
    # 3. Check for Multicollinearity
    from statsmodels.stats.outliers_influence import variance_inflation_factor
    # (calculation for VIF...)
    
    return model.summary()

2. 市场组合建模（归因）

def estimate_attribution(df):
    # Log-log model to calculate elasticities
    # Coefficients will be % change in sales for 1% change in spend
    model = smf.ols('np.log(sales) ~ np.log(tv_spend) + np.log(digital_spend)', data=df).fit()
    return model.params

3. 生存分析

from statsmodels.duration.hazard_regression import PHReg

# Cox Proportional Hazards Model
model = PHReg.from_formula('time ~ age + C(treatment)', data=df, status=df['event'])
results = model.fit()
print(results.summary())

性能优化

使用 numba 计算似然

Statsmodels 主要使用 Python 和 Cython 编写，但一些较新的时间序列模块会使用经过优化的数值后端，以更快地拟合状态空间模型。

公式与设计矩阵

对于非常大的数据集（100 万行以上），使用 patsy 创建设计矩阵可能会大量消耗内存。在这种情况下，请手动构造 X 矩阵，并使用 sm.OLS(y, X)。

常见问题与解决方案

奇异矩阵错误

“LinAlgError: Singular matrix” 表示预测变量之间存在完全相关性（例如，同时包含 temp_celsius 和 temp_fahrenheit）。

# ✅ Solution: Remove redundant columns
df = df.drop('redundant_col', axis=1)

分类变量泄漏（虚拟变量陷阱）

同时包含截距项和所有类别的虚拟变量会导致完全多重共线性。

# ❌ Problem: y ~ dummy_cat1 + dummy_cat2 + dummy_cat3 + intercept 
# ✅ Solution: Statsmodels/Patsy automatically drops one category (ref category)
# to avoid this. Don't try to force all dummies in!

非平稳时间序列

预测非平稳序列会导致“伪回归”。

# ✅ Solution: Difference your data first
df['diff_y'] = df['y'].diff()
# Or use ARIMA with integrated term (d=1)

Statsmodels 是 Python 生态系统中统计有效性的黄金标准。它超越了黑盒预测，提供了高风险科学和经济决策所需的透明度与数学严谨性。