1.3 OLS 回归详解：从单变量到多变量

"Correlation is not causation but it sure is a hint.""相关不是因果，但它确实是一个提示。"— Edward Tufte, Statistician (统计学家)

掌握 Python 中最常用的回归方法

本节目标

使用 statsmodels 进行 OLS 回归
理解回归输出的每一项指标
与 Stata/R 进行对比
提取和使用回归结果

案例：研究教育回报率

我们使用一个经典的劳动经济学问题：教育对工资的影响有多大？

python

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

# 生成模拟数据：1000 个观测
np.random.seed(2024)
n = 1000

data = pd.DataFrame({
    'wage': np.random.normal(5000, 2000, n),
    'education': np.random.randint(9, 22, n),
    'experience': np.random.randint(0, 30, n),
    'female': np.random.choice([0, 1], n),
    'urban': np.random.choice([0, 1], n)
})

# 让数据更符合现实：教育和经验正向影响工资
data['wage'] = (2000 +
                300 * data['education'] +
                50 * data['experience'] -
                800 * data['female'] +
                500 * data['urban'] +
                np.random.normal(0, 500, n))

print(data.head())
print(f"\n数据维度: {data.shape}")

模型 1：简单线性回归

先从最简单的模型开始：只用教育年限解释工资。

Python 代码

python

# 模型 1: wage = β0 + β1*education + ε
X1 = sm.add_constant(data['education'])
y = data['wage']

model1 = sm.OLS(y, X1).fit()
print(model1.summary())

关键输出：

                 coef    std err          t      P>|t|
------------------------------------------------------
const       2134.567     89.234     23.918      0.000
education    312.456     5.823      53.667      0.000
------------------------------------------------------
R-squared:                       0.743

Stata 对比

stata

* Stata 代码
regress wage education

R 对比

# R 代码
model1 <- lm(wage ~ education, data = data)
summary(model1)

解读：

教育年限每增加 1 年，工资增加约 312 元
R² = 0.743，教育解释了 74.3% 的工资变异
p < 0.001，教育的影响高度显著

模型 2：多元线性回归

加入更多控制变量：工作经验、性别、城乡。

python

# 模型 2: wage = β0 + β1*education + β2*experience + β3*female + β4*urban + ε
X2 = sm.add_constant(data[['education', 'experience', 'female', 'urban']])
y = data['wage']

model2 = sm.OLS(y, X2).fit()
print(model2.summary())

完整输出：

                            OLS Regression Results
==============================================================================
Dep. Variable:                   wage   R-squared:                       0.952
Model:                            OLS   Adj. R-squared:                  0.952
Method:                 Least Squares   F-statistic:                     4923.
Date:                ...              Prob (F-statistic):               0.00
No. Observations:                1000   AIC:                         1.234e+04
Df Residuals:                     995   BIC:                         1.236e+04
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const       2045.234     45.678     44.778      0.000    1955.570    2134.898
education    298.765      2.345    127.418      0.000     294.165     303.365
experience    49.876      1.234     40.418      0.000      47.456      52.296
female      -789.234     23.456    -33.642      0.000    -835.267    -743.201
urban        487.654     23.456     20.789      0.000     441.621     533.687
==============================================================================

解读：

教育：每多 1 年，工资增加 299 元（控制其他变量后）
经验：每多 1 年，工资增加 50 元
性别：女性工资比男性低 789 元（性别工资差距）
城乡：城市居民工资比农村高 488 元
R² 提升到 0.952，模型拟合度显著提高

模型对比：三种语言并列

Python (statsmodels)

python

import statsmodels.api as sm

# 准备数据
X = sm.add_constant(data[['education', 'experience', 'female', 'urban']])
y = data['wage']

# 拟合模型
model = sm.OLS(y, X).fit()

# 查看结果
print(model.summary())

# 提取系数
print(model.params)

# 提取 R²
print(f"R² = {model.rsquared:.4f}")

# 预测
predictions = model.predict(X)

Stata

stata

* 读取数据
use wage_data.dta, clear

* 拟合模型
regress wage education experience female urban

* 查看系数
display _b[education]

* 预测
predict wage_hat

R

# 读取数据
data <- read.csv("wage_data.csv")

# 拟合模型
model <- lm(wage ~ education + experience + female + urban, data = data)

# 查看结果
summary(model)

# 提取系数
coef(model)

# 提取 R²
summary(model)$r.squared

# 预测
predictions <- predict(model)

提取回归结果

获取系数

python

# 所有系数
print(model2.params)
'''
const        2045.234
education     298.765
experience     49.876
female       -789.234
urban         487.654
'''

# 单个系数
edu_coef = model2.params['education']
print(f"教育系数: {edu_coef:.2f}")

获取 p 值

python

# 所有 p 值
print(model2.pvalues)

# 判断显著性
if model2.pvalues['education'] < 0.05:
    print("教育系数在 5% 水平上显著")

获取置信区间

python

# 95% 置信区间
ci = model2.conf_int(alpha=0.05)
print(ci)
'''
                   0           1
const       1955.570    2134.898
education    294.165     303.365
experience    47.456      52.296
female      -835.267    -743.201
urban        441.621     533.687
'''

获取拟合优度

python

print(f"R² = {model2.rsquared:.4f}")
print(f"Adjusted R² = {model2.rsquared_adj:.4f}")
print(f"F-statistic = {model2.fvalue:.2f}")
print(f"Prob(F-statistic) = {model2.f_pvalue:.4e}")

预测新观测

python

# 创建新数据：一个 18 年教育、5 年经验的城市男性
new_data = pd.DataFrame({
    'const': [1],
    'education': [18],
    'experience': [5],
    'female': [0],
    'urban': [1]
})

# 预测工资
predicted_wage = model2.predict(new_data)
print(f"预测工资: {predicted_wage.values[0]:.2f} 元")

# 预测区间
predictions = model2.get_prediction(new_data)
pred_summary = predictions.summary_frame(alpha=0.05)
print(pred_summary)
'''
       mean    mean_se  mean_ci_lower  mean_ci_upper  obs_ci_lower  obs_ci_upper
0  8234.56      23.45        8188.59        8280.53       7245.67       9223.45
'''

解读：

点预测：8234.56 元
95% 置信区间：[8188.59, 8280.53]（均值的不确定性）
95% 预测区间：[7245.67, 9223.45]（个体的不确定性，更宽）

三语言对比总结

任务	Python	Stata	R
拟合模型	`sm.OLS(y, X).fit()`	`regress y x1 x2`	`lm(y ~ x1 + x2)`
查看结果	`model.summary()`	自动显示	`summary(model)`
获取系数	`model.params`	`_b[x1]`	`coef(model)`
获取 p 值	`model.pvalues`	`test x1`	`summary(model)$coefficients`
获取 R²	`model.rsquared`	`e(r2)`	`summary(model)$r.squared`
预测	`model.predict(X)`	`predict yhat`	`predict(model)`
置信区间	`model.conf_int()`	`regress, level(95)`	`confint(model)`

️ 常见错误

错误 1：忘记添加常数项

python

#  错误
X = data[['education', 'experience']]
model = sm.OLS(y, X).fit()  # 没有截距项！

#  正确
X = sm.add_constant(data[['education', 'experience']])
model = sm.OLS(y, X).fit()

错误 2：变量名拼写错误

python

#  错误
model.params['educaton']  # 拼写错误

#  正确
model.params['education']

错误 3：混淆 X 和 y 的顺序

python

#  错误（会报错）
model = sm.OLS(X, y).fit()  # 顺序反了

#  正确
model = sm.OLS(y, X).fit()  # y 在前，X 在后

实践练习

运行以下完整代码：

python

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

# 生成数据
np.random.seed(42)
n = 500
data = pd.DataFrame({
    'wage': np.random.normal(4000, 1500, n),
    'education': np.random.randint(9, 22, n),
    'age': np.random.randint(22, 60, n),
    'female': np.random.choice([0, 1], n)
})

# 构造真实关系
data['wage'] = (1000 + 250 * data['education'] +
                30 * data['age'] - 500 * data['female'] +
                np.random.normal(0, 400, n))

# 拟合模型
X = sm.add_constant(data[['education', 'age', 'female']])
y = data['wage']
model = sm.OLS(y, X).fit()

# 查看结果
print(model.summary())

# 提取关键信息
print(f"\n教育回报率: {model.params['education']:.2f} 元/年")
print(f"性别工资差距: {model.params['female']:.2f} 元")
print(f"R² = {model.rsquared:.4f}")

下一步

文章 03：Logit 回归 - 处理二元因变量
文章 04：summary_col() - 优雅地展示多个模型

** 你已经掌握了 Python 中的 OLS 回归！**

1.3 OLS 回归详解：从单变量到多变量 ​

本节目标 ​

案例：研究教育回报率 ​

模型 1：简单线性回归 ​

Python 代码 ​

Stata 对比 ​

R 对比 ​

模型 2：多元线性回归 ​

模型对比：三种语言并列 ​

Python (statsmodels) ​

Stata ​

R ​

提取回归结果 ​

获取系数 ​

获取 p 值 ​

获取置信区间 ​

获取拟合优度 ​

预测新观测 ​

三语言对比总结 ​

️ 常见错误 ​

错误 1：忘记添加常数项 ​

错误 2：变量名拼写错误 ​

错误 3：混淆 X 和 y 的顺序 ​

实践练习 ​

下一步 ​

1.3 OLS 回归详解：从单变量到多变量

本节目标

案例：研究教育回报率

模型 1：简单线性回归

Python 代码

Stata 对比

R 对比

模型 2：多元线性回归

模型对比：三种语言并列

Python (statsmodels)

Stata

R

提取回归结果

获取系数

获取 p 值

获取置信区间

获取拟合优度

预测新观测

三语言对比总结

️ 常见错误

错误 1：忘记添加常数项

错误 2：变量名拼写错误

错误 3：混淆 X 和 y 的顺序

实践练习

下一步