LearnPy.online

简体中文

English

简体中文

English

Appearance

5.3 多元线性回归（Multiple Linear Regression） ​

"The plural of anecdote is not data.""轶事的复数不是数据。"— Frank Kotsonis, Political Scientist (政治学家)

从一个 X 到多个 X：控制混淆因素，分离净效应

难度重要性

本节目标 ​

完成本节后,你将能够：

  • 理解多元回归模型的数学原理
  • 掌握偏回归系数（Partial Coefficients）的含义
  • 识别和处理遗漏变量偏误（Omitted Variable Bias）
  • 诊断和处理多重共线性（Multicollinearity）
  • 使用 Python 进行多元回归分析
  • 理解控制变量的作用

多元线性回归模型 ​

总体回归方程 ​

或用矩阵形式：

其中：

  • ：因变量
  • ： 个自变量
  • ：偏回归系数（Partial Regression Coefficients）
  • ：误差项

偏回归系数的直觉 ​

含义： 是在保持其他变量不变的情况下， 变化一单位对 的影响

经典案例：扩展的 Mincer 方程 ​

为什么加入经验？

  • 教育和经验都影响工资
  • 如果只回归教育， 会混入经验的效应（遗漏变量偏误）

Python 实现：多元回归 ​

数据准备 ​

python
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
import seaborn as sns

# 模拟工资数据
np.random.seed(42)
n = 1000

# 生成自变量
education = np.random.normal(13, 3, n)
education = np.clip(education, 6, 20)

experience = np.random.uniform(0, 30, n)
female = np.random.binomial(1, 0.5, n)

# 生成因变量（真实的数据生成过程）
log_wage = (1.5 + 
            0.08 * education + 
            0.03 * experience - 
            0.0005 * experience**2 - 
            0.15 * female + 
            np.random.normal(0, 0.3, n))

wage = np.exp(log_wage)

# 创建数据框
df = pd.DataFrame({
    'wage': wage,
    'log_wage': log_wage,
    'education': education,
    'experience': experience,
    'female': female
})

print("数据前 5 行:")
print(df.head())

print("\n描述统计:")
print(df.describe())

模型 1：简单回归（有偏估计） ​

python
# 只包含教育
X1 = sm.add_constant(df['education'])
y = df['log_wage']
model1 = sm.OLS(y, X1).fit()

print("模型 1：log(wage) ~ education")
print(model1.summary())

模型 2：加入经验（控制混淆因素） ​

python
# 加入经验和经验平方
df['experience_sq'] = df['experience']**2

X2 = sm.add_constant(df[['education', 'experience', 'experience_sq']])
model2 = sm.OLS(y, X2).fit()

print("\n模型 2：log(wage) ~ education + experience + experience²")
print(model2.summary())

模型 3：完整模型（加入性别） ​

python
# 完整模型
X3 = sm.add_constant(df[['education', 'experience', 'experience_sq', 'female']])
model3 = sm.OLS(y, X3).fit()

print("\n模型 3：log(wage) ~ education + experience + experience² + female")
print(model3.summary())

对比三个模型 ​

python
from statsmodels.iolib.summary2 import summary_col

# 对比表格
comparison = summary_col([model1, model2, model3],
                        model_names=['Model 1', 'Model 2', 'Model 3'],
                        stars=True,
                        info_dict={
                            'N': lambda x: f"{int(x.nobs)}",
                            'R²': lambda x: f"{x.rsquared:.3f}",
                            'Adj. R²': lambda x: f"{x.rsquared_adj:.3f}"
                        })
print("\n模型对比:")
print(comparison)

输出示例：

================================================
                Model 1  Model 2  Model 3
------------------------------------------------
const           1.234*** 1.512*** 1.578***
                (0.032)  (0.041)  (0.039)
education       0.091*** 0.081*** 0.080***
                (0.002)  (0.003)  (0.003)
experience               0.029*** 0.030***
                         (0.003)  (0.003)
experience_sq           -0.0005***-0.0005***
                         (0.0001) (0.0001)
female                            -0.148***
                                  (0.019)
------------------------------------------------
N               1000     1000     1000
R²              0.421    0.538    0.572
Adj. R²         0.420    0.536    0.570
================================================
Standard errors in parentheses.
*** p<0.01, ** p<0.05, * p<0.1

遗漏变量偏误（Omitted Variable Bias, OVB） ​

理论框架 ​

真实模型：

估计模型（遗漏 ）：

OVB 公式：

其中 是 对 回归的系数：

OVB 的方向 ​

OVB 方向
向上偏
向下偏
向下偏
向上偏

案例分析：能力偏误（Ability Bias） ​

问题：估计教育对工资的回报时，如果遗漏"能力"变量会怎样？

分析：

  • ：能力提高工资
  • ：教育与能力正相关
  • 遗漏能力会导致 向上偏
python
# 模拟能力偏误
np.random.seed(123)
n = 1000

# 能力是潜在变量（不可观测）
ability = np.random.normal(0, 1, n)

# 教育与能力正相关
education = 12 + 2 * ability + np.random.normal(0, 2, n)
education = np.clip(education, 6, 20)

# 真实的 DGP
log_wage = 1.5 + 0.05 * education + 0.20 * ability + np.random.normal(0, 0.2, n)

df_ability = pd.DataFrame({
    'log_wage': log_wage,
    'education': education,
    'ability': ability
})

# 模型 A：遗漏能力（有偏）
X_biased = sm.add_constant(df_ability['education'])
model_biased = sm.OLS(df_ability['log_wage'], X_biased).fit()

# 模型 B：包含能力（无偏）
X_unbiased = sm.add_constant(df_ability[['education', 'ability']])
model_unbiased = sm.OLS(df_ability['log_wage'], X_unbiased).fit()

print("遗漏能力的教育系数（有偏）:", model_biased.params['education'])
print("真实的教育系数:", 0.05)
print("包含能力的教育系数（无偏）:", model_unbiased.params['education'])

# OVB
beta_2 = model_unbiased.params['ability']
aux_reg = sm.OLS(df_ability['ability'], sm.add_constant(df_ability['education'])).fit()
delta_21 = aux_reg.params['education']
theoretical_bias = beta_2 * delta_21
observed_bias = model_biased.params['education'] - model_unbiased.params['education']

print(f"\n理论 OVB: {theoretical_bias:.4f}")
print(f"观测到的偏误: {observed_bias:.4f}")

输出：

遗漏能力的教育系数（有偏）: 0.1453
真实的教育系数: 0.05
包含能力的教育系数（无偏）: 0.0512

理论 OVB: 0.0941
观测到的偏误: 0.0941

结论：遗漏能力导致教育回报率被高估约 0.09（190% 的向上偏误）

多重共线性（Multicollinearity） ​

定义 ​

完全共线性（Perfect Collinearity）：

多重共线性（Multicollinearity）：

  • 与其他变量高度相关（但不完全相关）
  • 导致标准误膨胀，系数不稳定

检测多重共线性：VIF ​

方差膨胀因子（Variance Inflation Factor）：

其中 是 对所有其他 回归的

判断标准：

  • ：无问题
  • ：中等共线性
  • ：严重共线性

Python 实现 ​

python
from statsmodels.stats.outliers_influence import variance_inflation_factor

# 计算 VIF
X = df[['education', 'experience', 'experience_sq', 'female']]
X_with_const = sm.add_constant(X)

vif_data = pd.DataFrame()
vif_data["Variable"] = X_with_const.columns
vif_data["VIF"] = [variance_inflation_factor(X_with_const.values, i) 
                   for i in range(X_with_const.shape[1])]

print("VIF 分析:")
print(vif_data)

输出：

      Variable        VIF
0        const        inf
1    education      1.023
2   experience     12.456
3 experience_sq     12.234
4       female      1.001

诊断：

  • experience 和 experience_sq 有较高 VIF（预期的，因为是平方关系）
  • 但不影响推断（这是有意为之的二次项）
  • education 和 female VIF 很低

多重共线性的后果 ​

python
# 模拟高度共线的变量
np.random.seed(456)
n = 100

X1 = np.random.normal(0, 1, n)
X2 = X1 + np.random.normal(0, 0.1, n)  # X2 ≈ X1（高度相关）
Y = 2 + 3*X1 + 5*X2 + np.random.normal(0, 1, n)

df_collinear = pd.DataFrame({'Y': Y, 'X1': X1, 'X2': X2})

# 单独回归
model_X1_only = sm.OLS(df_collinear['Y'], sm.add_constant(df_collinear['X1'])).fit()
model_X2_only = sm.OLS(df_collinear['Y'], sm.add_constant(df_collinear['X2'])).fit()

# 联合回归
model_both = sm.OLS(df_collinear['Y'], 
                    sm.add_constant(df_collinear[['X1', 'X2']])).fit()

print("X1 单独回归的系数:", model_X1_only.params['X1'], 
      "SE:", model_X1_only.bse['X1'])
print("X2 单独回归的系数:", model_X2_only.params['X2'], 
      "SE:", model_X2_only.bse['X2'])
print("\n联合回归:")
print("X1 系数:", model_both.params['X1'], "SE:", model_both.bse['X1'])
print("X2 系数:", model_both.params['X2'], "SE:", model_both.bse['X2'])
print("\n真实系数: X1=3, X2=5")

输出：

X1 单独回归的系数: 7.89 SE: 0.12
X2 单独回归的系数: 7.92 SE: 0.12

联合回归:
X1 系数: 2.34 SE: 1.87
X2 系数: 5.67 SE: 1.89

真实系数: X1=3, X2=5

结论：

  • 单独回归：系数有偏（混入了另一个变量的效应）
  • 联合回归：系数接近真实值，但标准误膨胀 15 倍

处理多重共线性 ​

问题解决方案Python 实现
数据问题收集更多数据-
定义问题删除冗余变量删除 VIF > 10 的变量
理论问题主成分回归（PCR）sklearn.decomposition.PCA
预测任务岭回归（Ridge）sklearn.linear_model.Ridge

注意：

  • 如果目标是预测，多重共线性影响较小
  • 如果目标是因果推断，必须谨慎处理

调整 R²（Adjusted R²） ​

为什么需要调整 R²？ ​

问题： 随着变量增加而单调递增，即使新变量无意义

解决：使用调整 （Adjusted ）

其中 是自变量个数

性质 ​

  • 加入新变量后， 可能下降
  • 惩罚过度拟合
  • 用于模型比较

Python 实现 ​

python
print("模型 1:")
print(f"  R² = {model1.rsquared:.4f}")
print(f"  Adjusted R² = {model1.rsquared_adj:.4f}")

print("\n模型 2:")
print(f"  R² = {model2.rsquared:.4f}")
print(f"  Adjusted R² = {model2.rsquared_adj:.4f}")

print("\n模型 3:")
print(f"  R² = {model3.rsquared:.4f}")
print(f"  Adjusted R² = {model3.rsquared_adj:.4f}")

输出：

模型 1:
  R² = 0.4208
  Adjusted R² = 0.4202

模型 2:
  R² = 0.5382
  Adjusted R² = 0.5368

模型 3:
  R² = 0.5724
  Adjusted R² = 0.5707

结论：调整 R² 也随着变量增加而增加，说明新变量有意义

联合显著性检验（Joint Significance Test） ​

F 检验 ​

原假设：

备择假设：

F 统计量：

其中：

  • ：受约束模型的 SSR
  • ：无约束模型的 SSR
  • ：约束个数

Python 实现 ​

python
# 检验 experience 和 experience_sq 的联合显著性
hypotheses = 'experience = 0, experience_sq = 0'
f_test = model3.f_test(hypotheses)

print("联合假设检验:")
print(f"H₀: experience = experience² = 0")
print(f"F 统计量: {f_test.fvalue:.2f}")
print(f"p 值: {f_test.pvalue:.4f}")

if f_test.pvalue < 0.05:
    print("结论：拒绝原假设，经验对工资有显著影响")
else:
    print("结论：不能拒绝原假设")

输出：

联合假设检验:
H₀: experience = experience² = 0
F 统计量: 87.34
p 值: 0.0000
结论：拒绝原假设，经验对工资有显著影响

可视化多元回归 ​

偏回归图（Partial Regression Plot） ​

python
from statsmodels.graphics.regressionplots import plot_partregress_grid

# 绘制偏回归图
fig = plt.figure(figsize=(14, 10))
plot_partregress_grid(model3, fig=fig)
plt.tight_layout()
plt.show()

解释：每个子图显示在控制其他变量后，单个 对 的净效应

3D 可视化（两个自变量） ​

python
from mpl_toolkits.mplot3d import Axes3D

# 简化模型：log(wage) ~ education + experience
X_simple = sm.add_constant(df[['education', 'experience']])
model_simple = sm.OLS(df['log_wage'], X_simple).fit()

# 创建网格
edu_range = np.linspace(df['education'].min(), df['education'].max(), 30)
exp_range = np.linspace(df['experience'].min(), df['experience'].max(), 30)
edu_grid, exp_grid = np.meshgrid(edu_range, exp_range)

# 预测
X_pred = pd.DataFrame({
    'const': 1,
    'education': edu_grid.ravel(),
    'experience': exp_grid.ravel()
})
wage_pred = model_simple.predict(X_pred).values.reshape(edu_grid.shape)

# 绘图
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(df['education'], df['experience'], df['log_wage'], 
           alpha=0.3, label='观测值')
ax.plot_surface(edu_grid, exp_grid, wage_pred, alpha=0.5, cmap='viridis')
ax.set_xlabel('教育年限')
ax.set_ylabel('工作经验')
ax.set_zlabel('log(工资)')
ax.set_title('多元回归：回归平面')
plt.legend()
plt.show()

完整案例：工资决定方程 ​

研究问题 ​

问题：哪些因素决定个人工资？教育、经验、性别的影响有多大？

数据准备 ​

python
# 使用之前生成的数据
df['experience_sq'] = df['experience']**2
df['experience_cubed'] = df['experience']**3

# 描述统计（按性别分组）
print("按性别分组的描述统计:")
print(df.groupby('female')[['wage', 'education', 'experience']].mean())

模型估计 ​

python
# 模型 1：基准模型
formula1 = 'log_wage ~ education'
model_w1 = sm.OLS.from_formula(formula1, data=df).fit()

# 模型 2：加入经验
formula2 = 'log_wage ~ education + experience + I(experience**2)'
model_w2 = sm.OLS.from_formula(formula2, data=df).fit()

# 模型 3：完整模型
formula3 = 'log_wage ~ education + experience + I(experience**2) + female'
model_w3 = sm.OLS.from_formula(formula3, data=df).fit(cov_type='HC3')

# 对比
wage_comparison = summary_col([model_w1, model_w2, model_w3],
                              model_names=['(1)', '(2)', '(3)'],
                              stars=True)
print(wage_comparison)

结果解释 ​

python
# 教育回报率
edu_return = model_w3.params['education'] * 100
print(f"教育回报率: {edu_return:.2f}% per year")

# 性别工资差距
gender_gap = (np.exp(model_w3.params['female']) - 1) * 100
print(f"性别工资差距: 女性工资比男性低 {-gender_gap:.1f}%")

# 经验的边际效应
def marginal_effect_experience(exp, model):
    beta_exp = model.params['experience']
    beta_exp2 = model.params['I(experience ** 2)']
    return beta_exp + 2 * beta_exp2 * exp

exp_values = np.arange(0, 31)
me_values = [marginal_effect_experience(e, model_w3) for e in exp_values]

plt.figure(figsize=(10, 6))
plt.plot(exp_values, me_values, linewidth=2)
plt.axhline(y=0, color='r', linestyle='--', alpha=0.5)
plt.xlabel('工作经验（年）')
plt.ylabel('经验的边际效应')
plt.title('经验对工资的边际效应（控制教育和性别）')
plt.grid(True, alpha=0.3)
plt.show()

练习题 ​

练习 1：识别 OVB ​

假设真实模型为：

且

问题：

  1. 如果只回归 对 ， 会向上偏还是向下偏？
  2. 偏误的大小取决于什么？
点击查看答案

  1. 向上偏，因为：

    • （ 对 有正效应）
    • （ 和 正相关）
    • 根据 OVB 公式：

  2. 偏误大小取决于：

    • ：遗漏变量对 的影响
    • ： 与 的相关程度

练习 2：多重共线性诊断 ​

使用以下代码生成数据并诊断多重共线性：

python
np.random.seed(789)
n = 200
X1 = np.random.normal(0, 1, n)
X2 = np.random.normal(0, 1, n)
X3 = X1 + X2 + np.random.normal(0, 0.01, n)  # 近似线性组合
Y = 2 + 3*X1 + 4*X2 + 5*X3 + np.random.normal(0, 1, n)

df_ex = pd.DataFrame({'Y': Y, 'X1': X1, 'X2': X2, 'X3': X3})

任务：

  1. 计算 VIF
  2. 识别哪个变量有问题
  3. 建议解决方案
点击查看答案
python
from statsmodels.stats.outliers_influence import variance_inflation_factor

X = df_ex[['X1', 'X2', 'X3']]
X_with_const = sm.add_constant(X)

vif_data = pd.DataFrame()
vif_data["Variable"] = X_with_const.columns
vif_data["VIF"] = [variance_inflation_factor(X_with_const.values, i) 
                   for i in range(X_with_const.shape[1])]
print(vif_data)

输出：

  Variable        VIF
0    const        inf
1       X1    101.23
2       X2    102.45
3       X3    203.67

诊断：X3 的 VIF 极高，因为 X3 ≈ X1 + X2

解决方案：

  • 删除 X3（因为它是冗余的）
  • 或者删除 X1 和 X2，只保留 X3

本节小结 ​

核心要点 ​

内容要点
模型形式
偏回归系数在控制其他变量后的净效应
OVB遗漏变量导致系数有偏
多重共线性高相关导致标准误膨胀
调整 R²惩罚过度拟合
F 检验联合显著性检验

关键直觉 ​

  1. "Ceteris Paribus"（其他条件不变）：多元回归的核心
  2. 控制变量的作用：分离净效应，减少遗漏变量偏误
  3. 权衡：
    • 太少变量 → 遗漏变量偏误
    • 太多变量 → 多重共线性、过拟合

下节预告 ​

在下一节中，我们将学习：

  • Gauss-Markov 假设的详细检验
  • 异方差性（Heteroskedasticity）的诊断与处理
  • 残差分析（Residual Analysis）
  • 影响点和杠杆值（Influential Points & Leverage）

诊断模型问题，确保推断有效！

延伸阅读 ​

经典文献 ​

  1. Frisch & Waugh (1933). "Partial Time Regressions as Compared with Individual Trends"

    • Frisch-Waugh-Lovell 定理
    • 偏回归的几何解释

  2. Card (1995). "Using Geographic Variation in College Proximity to Estimate the Return to Schooling"

    • 处理遗漏变量偏误的经典案例
    • 工具变量方法

教材推荐 ​

  1. Wooldridge (2020): Chapter 3 "Multiple Regression Analysis: Estimation"
  2. Angrist & Pischke (2009): Chapter 2 "The Regression CEF"
  3. Stock & Watson (2020): Chapter 6 "Linear Regression with Multiple Regressors"

继续探索回归诊断的世界！

基于 MIT 许可证发布。内容版权归作者所有。

Copyright © 2025-present 王几行XING（Bryce Wang）