使用statsmodels进行普通最小二乘(OLS)回归分析实战指南

前言

statsmodels是Python中一个强大的统计分析库，提供了多种统计模型和检验方法。本文将重点介绍如何使用statsmodels进行普通最小二乘(OLS)回归分析，通过实际案例演示从基础到进阶的应用。

基础OLS回归分析

1. 数据准备与模型构建

首先我们创建一个简单的线性回归模型。假设y与x存在二次关系：

import numpy as np
import statsmodels.api as sm

# 设置随机种子保证结果可复现
np.random.seed(9876789)

# 生成100个样本点
nsample = 100
x = np.linspace(0, 10, 100)
X = np.column_stack((x, x**2))  # 包含x和x²
beta = np.array([1, 0.1, 10])   # 真实系数
e = np.random.normal(size=nsample)  # 误差项

# 添加常数项(截距)
X = sm.add_constant(X)
y = np.dot(X, beta) + e  # 生成y值

2. 模型拟合与结果解读

model = sm.OLS(y, X)  # 创建OLS模型
results = model.fit()  # 拟合模型
print(results.summary())  # 输出完整结果

模型结果摘要包含丰富信息：

• 系数估计值及其显著性检验
• R²和调整R²
• F检验结果
• 模型诊断信息

3. 关键指标提取

print("参数估计:", results.params)
print("R平方值:", results.rsquared)
print("调整R平方:", results.rsquared_adj)

非线性关系的线性模型

1. 构建非线性关系数据

nsample = 50
sig = 0.5
x = np.linspace(0, 20, nsample)
X = np.column_stack((x, np.sin(x), (x-5)**2, np.ones(nsample)))
beta = [0.5, 0.5, -0.02, 5.0]

y_true = np.dot(X, beta)
y = y_true + sig * np.random.normal(size=nsample)

2. 模型拟合与可视化

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

# 可视化结果
pred_ols = res.get_prediction()
iv_l = pred_ols.summary_frame()["obs_ci_lower"]
iv_u = pred_ols.summary_frame()["obs_ci_upper"]

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(x, y, "o", label="数据点")
ax.plot(x, y_true, "b-", label="真实关系")
ax.plot(x, res.fittedvalues, "r--.", label="OLS拟合")
ax.plot(x, iv_u, "r--")
ax.plot(x, iv_l, "r--")
ax.legend(loc="best")
plt.show()

分类变量的处理

1. 创建包含分类变量的数据

nsample = 50
groups = np.zeros(nsample, int)
groups[20:40] = 1
groups[40:] = 2

dummy = pd.get_dummies(groups).values
x = np.linspace(0, 20, nsample)
X = np.column_stack((x, dummy[:,1:]))  # 删除参考类别
X = sm.add_constant(X, prepend=False)

beta = [1.0, 3, -3, 10]
y_true = np.dot(X, beta)
y = y_true + np.random.normal(size=nsample)

2. 模型拟合与假设检验

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

# 联合假设检验
R = [[0, 1, 0, 0], [0, 0, 1, 0]]
print(res2.f_test(R))  # 检验两个分类变量系数是否同时为0

多重共线性问题

1. 经典案例：Longley数据集

from statsmodels.datasets.longley import load_pandas
data = load_pandas()
y = data.endog
X = data.exog
X = sm.add_constant(X)

# 拟合模型
ols_model = sm.OLS(y, X)
ols_results = ols_model.fit()
print(ols_results.summary())

2. 多重共线性诊断

# 计算条件数
norm_x = X.values
for i, name in enumerate(X):
    if name == "const":
        continue
    norm_x[:,i] = X[name]/np.linalg.norm(X[name])
norm_xtx = np.dot(norm_x.T, norm_x)

eigs = np.linalg.eigvals(norm_xtx)
condition_number = np.sqrt(eigs.max()/eigs.min())
print("条件数:", condition_number)  # 大于20表示严重多重共线性

3. 影响分析

# 删除一个观测值的影响
ols_results2 = sm.OLS(y.iloc[:14], X.iloc[:14]).fit()
print("参数变化百分比:", (ols_results2.params - ols_results.params)/ols_results.params*100)

# DFBETAS分析
infl = ols_results.get_influence()
print(infl.summary_frame().filter(regex="dfb"))  # 标准化系数变化

结语

本文通过statsmodels演示了OLS回归的完整流程，包括：

1. 基础线性回归实现
2. 非线性关系的线性建模
3. 分类变量处理与假设检验
4. 多重共线性诊断与影响分析

statsmodels提供了丰富的统计工具和诊断方法，是Python数据分析生态中不可或缺的一部分。掌握这些技术将大大提升您的数据分析能力。