基本多元回归
本章使用的数据
本章使用的数据主要放在工作目录下的 data 文件夹中。读取时用相对路径，例如 data/wage1.csv。
data/wage1.csv：伍德里奇《计量经济学导论》中的工资数据。一行是一名劳动者，包含工资、对数工资、教育年数、工作经验、在当前工作的年限、性别等变量，用于演示多元线性回归。下载data/jtrain2.csv：Lalonde 就业培训实验数据。一行是一名参加 National Supported Work Demonstration 项目的劳动者，包含是否被分配到就业培训、培训前收入、培训后收入、年龄、教育年数等变量，用于本章新版作业。train 表示是否被分配到培训，re78 表示 1978 年真实收入，单位为千美元。该数据的培训分配接近随机实验，因此比一般观察数据更适合讨论因果解释。下载，变量说明data/gpa1.csv：学生成绩数据。一行是一名学生，包含大学第一年 GPA、高中 GPA、ACT 成绩等变量。该数据保留为旧版作业示例，当前已标记为废弃。下载grunfeld：statsmodels 内置的企业投资面板数据，通过 statsmodels.datasets.grunfeld 读取，不需要单独下载，用于选学固定效应示例。
- 本章默认本地数据文件放在当前章节 notebook 同级目录下的
data 文件夹中。
- 出现
FileNotFoundError，通常是工作目录不对，或者数据文件没有放在 data 文件夹中。
- 固定效应部分属于选学内容，重点是看懂面板数据模型的基本形式和代码写法。
分析目标和数据
这个案例来自伍德里奇的《计量经济学导论》。
本案关心的是一个因果问题：教育水平是否会影响工资？即：每多读一年书，工资是否会提高？
要注意，提出因果问题并不等于已经完成因果论证。仅仅做一个横截面 OLS 回归，还不能自动证明教育对工资的因果效应。本章先用这个例子演示回归模型的建立、OLS 的运行和回归系数的解释；后续更严格的因果识别问题再继续展开。
数据采用《导论》中的数据Wage1。
# 导入包并设置绘图主题、色板、字体和输出格式
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
okabe_ito = ["#0072B2", "#E69F00", "#009E73", "#D55E00", "#CC79A7", "#56B4E9", "#999999"]
sns.set_theme(
style="ticks", # 使用带刻度、少装饰的主题
context="notebook", # 使用适合notebook阅读的字号
)
sns.set_palette(okabe_ito) # 使用色盲友好色板
plt.rcParams["font.sans-serif"] = ["Hiragino Sans GB", "Microsoft YaHei", "WenQuanYi Micro Hei"] # 支持中文
plt.rcParams["axes.unicode_minus"] = False # 正常显示负号
%config InlineBackend.figure_formats = ["svg"] # 使用SVG格式
# 加载并查看wage1.csv数据集
wage1_data = pd.read_csv("data/wage1.csv")
wage1_data.iloc[:5, :10] # 显示前5行和前10列
| 0 |
3.10 |
11 |
2 |
0 |
0 |
1 |
0 |
2 |
1 |
0 |
| 1 |
3.24 |
12 |
22 |
2 |
0 |
1 |
1 |
3 |
1 |
0 |
| 2 |
3.00 |
11 |
2 |
0 |
0 |
0 |
0 |
2 |
0 |
0 |
| 3 |
6.00 |
8 |
44 |
28 |
0 |
0 |
1 |
0 |
1 |
0 |
| 4 |
5.30 |
12 |
7 |
2 |
0 |
0 |
1 |
1 |
0 |
0 |
wage1_data.columns # 显示所有列名
Index(['wage', 'educ', 'exper', 'tenure', 'nonwhite', 'female', 'married',
'numdep', 'smsa', 'northcen', 'south', 'west', 'construc', 'ndurman',
'trcommpu', 'trade', 'services', 'profserv', 'profocc', 'clerocc',
'servocc', 'lwage', 'expersq', 'tenursq'],
dtype='object')
先看两个变量的关系
在进入多元回归之前，可以先用散点图观察教育年数和对数工资之间的简单关系。散点图的作用是帮助我们形成直观认识：教育年数较高的人，对数工资是否整体上也更高。
这里使用 sns.regplot()，而不是 sns.scatterplot()，是因为 regplot() 可以通过参数加入一条简单回归线，便于观察整体趋势。
但这一步还不能完成因果论证。教育年数不同的人，可能在性别、工作经验、在当前工作的年限等方面也不同。下面的散点图只是起点，后面的多元回归才会把这些变量加入模型，让教育系数具有“其他条件保持不变”的含义。
这里使用对数工资 lwage，是为了和后面的正式回归模型保持一致。工资通常右偏，取对数后更适合线性回归；同时，系数也可以近似解释为工资的百分比变化。
# 绘制教育年数和对数工资的散点图，并加入简单回归线
fig, ax = plt.subplots(figsize=(8, 5))
sns.regplot(
data=wage1_data,
x="educ",
y="lwage",
fit_reg=True, # 加入简单回归线，便于观察整体趋势
scatter_kws={"alpha": 0.35}, # alpha设置散点透明度，减轻点重叠
line_kws={"color": okabe_ito[3], "linewidth": 2}, # 设置回归线颜色和粗细
ax=ax,
)
ax.set_title("教育年数和对数工资的简单关系")
ax.set_xlabel("教育年数")
ax.set_ylabel("对数工资")
plt.show()
从图中可以看到，教育年数和对数工资基本呈现正相关：教育年数较高的样本，对数工资整体上也更高。图中的回归线向右上方倾斜，说明两个变量之间存在正向的简单关系。
不过，散点分布也比较分散，同一教育年数下的工资差异仍然很大。这说明工资还受到其他因素影响，例如工作经验、在当前工作的年限、性别、行业和地区等。因此，散点图只能作为初步观察，不能直接完成因果解释。下面要进入多元回归，把一些重要因素作为控制变量加入模型。
注意：为什么到这里还不行？
我们说的“因果关系”，其实是在说：\(x\) 的变化，引起了 \(y\) 的变化。放在这里，就是教育年数的增加，会引起工资的增加。
但很多因素，比如“智力”，可能同时影响教育年数和工资。会不会是这样：智力更高的人，更容易读更多年书，同时工作能力也更强，所以工资更高？如果是这样，我们观察到的正相关，就不一定是“教育影响工资”。甚至有可能教育本身并不影响工资，只是智力同时影响了教育和工资，导致数据里“看起来”像是教育提高了工资。
这就是遗漏变量问题：一个重要因素没有进入模型，但它同时影响解释变量和被解释变量。此时，简单相关关系不能直接解释为因果关系。
当然，这只是其中一种问题。现实研究里还可能有选择偏误、反向因果、测量误差等其他问题。这里先抓住最关键的一点：正相关不等于因果关系。
为什么要加入控制变量
那怎么处理这种问题？一个最直观的方法是：只在“智力相同”的人之间比较。
如果在智力相同的人里面，还是发现读书更多的人，工资也更高，那就可以排除“智力”这个因素的影响。至少在这个比较里，工资差异就不能再简单归因于智力差异。
这就是控制变量。这里的“控制”是一个动词，意思是把某个因素固定住，尽量只比较其他条件相同的人。
对于教育和工资这个问题，如果我们把所有可能影响工资的因素都控制起来，可以理解为全部都在“同一个水平内比较”：同智力、同地区、同年代、同家庭背景、同工作经验，等等。假设我们能把这些因素造成的影响都排除掉，剩下的差异主要就是教育年数不同，那么我们就可以说：在其他因素保持不变的情况下，教育水平会影响收入。这就是这里说的因果关系。
当然，实际研究没有这么理想。有些信息永远拿不到，比如智力、家庭环境、个人选择；这就是遗漏变量问题。很多时候，\(x\) 和 \(y\) 还会互相影响，这叫反向因果。除此之外还有其他问题。所以这个问题还没完，但已经不是本节要展开的内容。
本节先做最基本的一步：把数据里能直接观察到的一些变量放进模型。下面的例子会控制性别、工作经验和在当前工作的年限。
回归公式
基于这个思路，可以建立以下多元回归模型。这里使用对数工资 lwage 作为因变量：
\[
\text{lwage} = \beta_0 + \delta_0 \cdot \text{female} + \beta_1 \cdot \text{educ} + \beta_2 \cdot \text{exper} + \beta_3 \cdot \text{tenure} + u
\]
其中：
- \(\text{lwage}\) 是工资的自然对数。对数工资为负并不表示工资为负；只要原始工资是正数，就可以取对数。当原始工资小于 1 时，对数值会小于 0。
- \(\text{female}\) 是一个虚拟变量，女性为1，男性为0。
- \(\text{educ}\), \(\text{exper}\), 和 \(\text{tenure}\) 分别代表教育年数、工作经验和在当前工作的年限。
- \(u\) 是误差项。
从理论上看，这类线性回归模型通常可以写成一般形式： \[
y = \beta_0 + \beta_1 x_1 + \cdots + \beta_k x_k + u,
\] 其中 \(E(u\mid x_1,\dots,x_k)=0\) 是一个核心假设，它表示在给定解释变量之后，误差项的条件期望为零。普通最小二乘法（OLS）就是通过让残差平方和尽量小来估计这些系数。
教育水平对工资的影响体现在系数\(\beta_1\)中，这就是我们关心的目标。由于因变量是对数工资，\(\beta_1\) 可以近似解释为教育年数增加一年时，工资变化的百分比。
进行回归分析
在进行回归分析之前，我们需要准备数据和定义模型。这个过程包括设置自变量（解释变量）和因变量（被解释变量），以及确保数据格式适合进行回归分析。
定义因变量和自变量
因变量 \(y\): 因变量是我们想要预测或解释的变量。在本例中，因变量是 lwage，即对数工资。我们希望了解其他变量如何影响工资的百分比变化。
自变量 \(X\): 自变量是可能影响因变量的变量。在这个模型中，我们选择了四个自变量：female, educ, exper, 和 tenure。这些变量分别表示：
female: 表示性别，是一个二分类变量，女性为1，男性为0。educ: 表示受教育年数，是一个连续变量。exper: 表示总工作经验年数，是一个连续变量。tenure: 表示在当前工作的年限，也是一个连续变量。
构造自变量 \(X\) 和因变量 \(y\)
在开始构建模型之前，我们需要从数据集中提取这些变量。使用Pandas库，我们可以从DataFrame中选择需要的列。同时，为了进行回归分析，我们需要向自变量矩阵 \(X\) 中添加一个常数项，以便模型包含截距 \(\beta_0\)。
import statsmodels.api as sm
# 选择自变量和因变量，对数工资lwage作为因变量
X = wage1_data[["female", "educ", "exper", "tenure"]]
y = wage1_data["lwage"]
# 向自变量矩阵X添加常数项，以便包含截距
X = sm.add_constant(X)
X.head() # 注意，多了一列const
| 0 |
1.0 |
1 |
11 |
2 |
0 |
| 1 |
1.0 |
1 |
12 |
22 |
2 |
| 2 |
1.0 |
0 |
11 |
2 |
0 |
| 3 |
1.0 |
0 |
8 |
44 |
28 |
| 4 |
1.0 |
0 |
12 |
7 |
2 |
拟合回归模型
使用statsmodels库的OLS类，我们可以定义并拟合一个普通最小二乘回归模型。OLS类的第一个参数是因变量 \(y\)，第二个参数是自变量 \(X\)。在前面的模型假设下，OLS 给出的系数是通过最小化残差平方和得到的估计值，可以用来构造样本回归函数（或预测方程）。
# 创建OLS回归模型的实例并拟合数据
model = sm.OLS(y, X)
results = model.fit()
查看回归结果
模型拟合后，我们可以通过调用 summary() 方法来查看详细的回归分析结果，包括每个变量的系数、标准误、t 值、p 值等统计指标。每个系数的 t 值和 p 值都对应于一个假设检验，通常是检验原假设 \(H_0: \beta_j = 0\)，用来判断该解释变量在统计上是否对因变量有影响。
# 输出回归结果摘要
print(results.summary())
OLS Regression Results
==============================================================================
Dep. Variable: lwage R-squared: 0.392
Model: OLS Adj. R-squared: 0.388
Method: Least Squares F-statistic: 84.07
Date: Mon, 08 Jun 2026 Prob (F-statistic): 4.68e-55
Time: 10:06:25 Log-Likelihood: -282.46
No. Observations: 526 AIC: 574.9
Df Residuals: 521 BIC: 596.2
Df Model: 4
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 0.5013 0.102 4.920 0.000 0.301 0.702
female -0.3011 0.037 -8.085 0.000 -0.374 -0.228
educ 0.0875 0.007 12.605 0.000 0.074 0.101
exper 0.0046 0.002 2.845 0.005 0.001 0.008
tenure 0.0174 0.003 5.835 0.000 0.012 0.023
==============================================================================
Omnibus: 12.037 Durbin-Watson: 1.775
Prob(Omnibus): 0.002 Jarque-Bera (JB): 22.360
Skew: 0.012 Prob(JB): 1.40e-05
Kurtosis: 4.010 Cond. No. 141.
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
稳健标准误
上面的回归表默认使用常规标准误。如果截面数据存在异方差，常规标准误可能不够可靠。一个常见做法是报告异方差稳健标准误。这里先不展开理论，只演示如何在 statsmodels 中取得 HC1 稳健标准误，并和常规标准误作比较。
# 比较常规标准误和HC1稳健标准误
robust_results = results.get_robustcov_results(cov_type="HC1")
robust_table = pd.DataFrame(
{
"coef": results.params,
"常规标准误": results.bse,
"HC1稳健标准误": robust_results.bse,
"常规p值": results.pvalues,
"HC1稳健p值": robust_results.pvalues,
},
index=results.params.index,
)
robust_table.round(3)
| const |
0.501 |
0.102 |
0.115 |
0.000 |
0.000 |
| female |
-0.301 |
0.037 |
0.038 |
0.000 |
0.000 |
| educ |
0.087 |
0.007 |
0.008 |
0.000 |
0.000 |
| exper |
0.005 |
0.002 |
0.002 |
0.005 |
0.004 |
| tenure |
0.017 |
0.003 |
0.004 |
0.000 |
0.000 |
样本回归函数
可以得到样本回归函数（或者预测方程）：
\[
\widehat{lwage} = 0.501 - 0.301 \cdot female + 0.087 \cdot educ + 0.005 \cdot exper + 0.017 \cdot tenure
\]
这里给出的就是样本回归函数，用估计得到的系数来给出在不同特征组合下的对数工资预测值。与前面的“包含误差项 \(u\) 的模型方程”不同，这里已经不再带误差项。
因为这是多元回归，严格说它不是二维图形中的一条“回归线”，而是包含多个解释变量的预测方程。
解释回归结果
从回归结果中，我们可以提取以下关键信息：
- 模型拟合情况: \(R^2\) 值为 0.392，说明模型解释了 39.2% 的因变量（
lwage）的方差。这表明模型有一定的解释能力，但在线性回归中，我们通常不会仅凭 \(R^2\) 来判断模型优劣，更关心的是系数的含义及其显著性。
- 系数解释:
- const（截距）: \(0.501\)。截距项在常规标准误和 HC1 稳健标准误下都显著。它在形式上表示当所有解释变量为0时预测的对数工资。在本例中，“教育年数=0、工作经验=0、在现单位年限=0、且为男性”的情形本身并不常见，因此截距更多是起到调整预测方程位置的技术作用，通常不作强经济解释。
- female: \(-0.301\)，显著为负。由于因变量是对数工资，可以先近似理解为：在控制教育、总工作经验和在现单位年限之后，女性工资低于男性。这个系数的绝对值不算很小，若要换算成更精确的百分比，可以计算 \(e^{-0.301}-1\approx -26.0\%\)。这不是简单的男女平均工资差，而是在这些控制变量相同条件下的差异。
- educ（教育年数）: \(0.087\)，显著为正。它是本例最关心的系数，因为它对应我们提出的因果问题：教育年数增加是否会提高工资。在这个简单回归中，如果暂时接受控制变量已经处理了主要混杂因素，那么可以按因果问题的语言理解为：在性别、总工作经验和当前工作年限相同的条件下，教育年数每增加一年，工资平均提高约 8.7%。
- exper（工作经验）: \(0.005\)，在 5% 水平下显著，表示在其他变量保持不变时，每增加一年总工作经验，工资平均提高约 0.5%。
- tenure（在当前工作的年限）: \(0.017\)，显著为正，表示在其他变量保持不变时，每增加一年在当前工作的年限，工资平均提高约 1.7%。
- 模型诊断:
- Durbin-Watson 统计量为 1.775，接近2，未显示明显的一阶自相关迹象。不过本例是截面数据，Durbin-Watson 不是最核心的诊断指标。
- Omnibus 和 Jarque-Bera 测试的 \(p\) 值较小，表明残差不完全符合正态分布。在截面回归中，也常结合稳健标准误来处理异方差带来的标准误问题。
回归结果的经济学解释
在经济学研究中，分析者往往特别关注模型中的系数（\(\beta_i\)），因为这些系数对应我们关心的因果问题。教育系数 \(0.087\) 的含义，不是简单地说“教育年数多的人工资更高”，而是：在性别、总工作经验和当前工作年限相同的条件下，教育年数每增加一年，工资平均提高约 8.7%。这种“其他条件保持不变”的解释，正是加入控制变量后的斜率含义。
但要注意，这个因果解释要成立，需要更强的论证：例如教育年数不能和遗漏在误差项中的能力、家庭背景、地区、职业选择等因素系统相关。本章先学习回归的基本操作和系数解读，后续再讨论如何更严格地支持这种因果解释。
对系数的解释
大小：系数的大小告诉我们，当自变量增加一个单位时，因变量平均会增加或减少多少。如果因变量是对数形式，系数通常可以近似解释为因变量的百分比变化。
方向：系数的符号（正或负）表明变量作用方向。例如，如果某个变量的系数为正，则表明该变量增加时，被解释变量平均上升；如果系数为负，则平均下降。
显著性：统计显著性（通常通过p值来衡量）告诉我们系数是否在统计上显著，即我们观察到的关系是否有可能仅仅是由于抽样误差。p值较小（通常小于0.05）意味着我们有足够的证据认为该系数在总体中不为零。
比较多个回归模型（进阶）
在经济学论文和研究报告中，经常把多个回归模型放在同一张表中比较。常见做法是：先放一个简单模型，再逐步加入控制变量，观察核心解释变量的系数、显著性和模型拟合情况是否发生变化。
本节使用 statsmodels.iolib.summary2 中的 summary_col()。它可以把多个回归结果整理成一张表，并用星号表示统计显著性。
本节用到的方法和函数：
summary_col(results, model_names, stars, float_format, info_dict)：把多个回归结果合并成一张表。results：放入多个已经拟合好的回归结果。model_names：设置每一列的模型名称。stars=True：显示显著性星号。float_format：控制小数位数。info_dict：在表格底部添加样本量等信息。
下面比较三个模型：
- 模型1：只放入教育年数
educ。
- 模型2：加入性别变量
female。
- 模型3：继续加入工作经验
exper 和在当前工作的年限 tenure。
# 比较逐步加入控制变量的三个回归模型
from statsmodels.iolib.summary2 import summary_col
X_m1 = sm.add_constant(wage1_data[["educ"]]) # 模型1：只加入教育年数
m1 = sm.OLS(wage1_data["lwage"], X_m1).fit()
X_m2 = sm.add_constant(wage1_data[["educ", "female"]]) # 模型2：加入性别控制变量
m2 = sm.OLS(wage1_data["lwage"], X_m2).fit()
X_m3 = sm.add_constant(wage1_data[["educ", "female", "exper", "tenure"]]) # 模型3：加入更多控制变量
m3 = sm.OLS(wage1_data["lwage"], X_m3).fit()
model_table = summary_col(
results=[m1, m2, m3], # 放入多个回归结果
model_names=["模型1", "模型2", "模型3"], # 设置表格列名
stars=True, # 显示显著性星号
float_format="%.3f", # 系数和标准误保留3位小数
info_dict={"N": lambda x: f"{int(x.nobs)}"}, # 在表格底部显示样本量
)
print(model_table)
===========================================
模型1 模型2 模型3
-------------------------------------------
const 0.584*** 0.826*** 0.501***
(0.097) (0.094) (0.102)
educ 0.083*** 0.077*** 0.087***
(0.008) (0.007) (0.007)
female -0.361*** -0.301***
(0.039) (0.037)
exper 0.005***
(0.002)
tenure 0.017***
(0.003)
R-squared 0.186 0.300 0.392
R-squared Adj. 0.184 0.298 0.388
N 526 526 526
===========================================
Standard errors in parentheses.
* p<.1, ** p<.05, ***p<.01
简单读表：三个模型中，educ 的系数都为正，并且都有显著性星号，说明教育年数和对数工资之间存在稳定的正向关系。模型1中，educ 的系数约为 0.083；加入 female 后，系数变为约 0.077；继续加入 exper 和 tenure 后，系数变为约 0.087。
这说明加入控制变量会影响系数的大小，也会改变我们对系数的解释。模型3中的 educ 系数含义更接近本章关心的问题：在性别、总工作经验和当前工作年限保持不变的条件下，教育年数每增加一年，工资平均提高约 8.7%。同时，\(R^2\) 从模型1的 0.186 提高到模型3的 0.392，说明加入控制变量后，模型解释的样本差异更多。
作业：就业培训是否提高收入
这道作业使用 data/jtrain2.csv。数据来自 Lalonde 的就业培训研究，研究对象是参加 National Supported Work Demonstration 项目的劳动者。这个项目中，是否被分配到就业培训具有实验性质，因此本题比普通的成绩预测练习更适合讨论因果关系。
研究问题：
就业培训是否会提高劳动者之后的收入？
主要变量：
train：是否被分配到就业培训，1 表示被分配到培训组，0 表示对照组。re78：1978 年真实收入，单位为千美元，是本题的因变量。age：年龄。educ：受教育年数。black：是否为黑人。hisp：是否为西班牙裔。married：是否已婚。nodegree：是否没有高中学历。re74、re75：1974 年和 1975 年真实收入，单位为千美元，表示培训前收入。
本题用到的方法和函数：
pd.read_csv()：读取本地 CSV 文件。DataFrame.describe()：查看变量的描述性统计。sm.add_constant()：为自变量矩阵添加截距项。sm.OLS(y, X).fit()：估计普通最小二乘回归模型。get_robustcov_results(cov_type="HC1")：得到 HC1 异方差稳健标准误。
任务：
请仿照前面的工资案例完成，不需要写很长的文字，但每一步都要说明你在做什么、结果应该怎么看。
- 读取数据，查看核心变量的描述性统计。
- 估计不加控制变量的模型：
\[
re78_i = \beta_0 + \beta_1 train_i + u_i
\]
- 估计加入控制变量的模型：
\[
\begin{aligned}
re78_i ={}& \beta_0 + \beta_1 train_i + \beta_2 age_i + \beta_3 educ_i \\
&+ \beta_4 black_i + \beta_5 hisp_i + \beta_6 married_i + \beta_7 nodegree_i \\
&+ \beta_8 re74_i + \beta_9 re75_i + u_i
\end{aligned}
\]
- 比较两个模型中
train 的系数方向、大小和显著性。
- 解释
train 的系数含义：在什么条件下，它可以被理解为就业培训对收入的因果影响？
- 简要说明加入控制变量以后，
train 的系数是否发生明显变化。
多元回归练习题（废弃）
这一版作业保留作为旧示例，不再作为本章推荐作业。它更适合练习条件相关关系和预测，不够适合作为因果关系训练。新版作业请使用上面的就业培训数据。
在一项研究中，研究人员收集了一批学生的高中成绩（hsGPA）、大学入学考试成绩（ACT）以及他们在大学第一年的成绩（colGPA）。该研究目的是探索高中GPA和ACT成绩对大学第一年GPA的影响。
数据集包含以下变量：
hsGPA：高中GPAACT：大学入学考试成绩colGPA：大学第一年GPA
研究问题： 高中GPA和ACT成绩如何影响大学第一年的GPA？
回归模型： 使用以下回归模型分析数据： \[
\text{colGPA} = \beta_0 + \beta_1 \times \text{hsGPA} + \beta_2 \times \text{ACT} + u
\]
其中：
- \(\beta_0\) 是截距
- \(\beta_1\) 是高中GPA对大学GPA的影响。
- \(\beta_2\) 是ACT成绩对大学GPA的影响。
- \(u\) 是误差项。
在这个模型中，可以把系数理解为在控制了另一个成绩之后，各自对大学第一年GPA的边际影响。需要注意的是，如果还有其他重要但未观测的因素同时影响 hsGPA、ACT 和 colGPA，这些遗漏变量可能会使系数偏离理想的因果效应，因此在解释时要结合实际背景进行判断。
任务：
- 利用提供的数据，估计上述回归模型。
- 解释高中GPA和ACT成绩的系数，包括意义和显著性。
- 思考题：讨论模型的适用性和限制。
固定效应（可选）
更进一步，我们可以使用 linearmodels来完成基准回归，并且加入固定效应。本节属于面板数据回归的预告，重点是看懂基本形式和代码写法。
在面板数据中，同一企业在不同年份反复出现，每个企业往往有一些不随时间变化但又影响因变量的特质，例如管理水平、企业文化等。如果这些特质与解释变量有关，简单的回归可能会有偏。固定效应模型通过为每个个体（以及每个年份）引入一个单独的截距，把这种“时间不变但个体不同”的因素吸收掉。
固定效应可以控制企业不随时间变化的特征，以及年份共同冲击，因此比完全不控制这些因素的模型更接近我们想要的因果比较。但它仍然不能自动解决所有问题：如果存在随时间变化、同时影响投资和企业价值的遗漏因素，或者存在反向因果，系数仍然不能简单视为严格因果效应。
如果导入 linearmodels 时出错，可以在Jupyter Notebook的python单元格中执行如下代码，注意开头有个感叹号：
!pip install -i https://pypi.tuna.tsinghua.edu.cn/simple linearmodels
数据准备
这里使用来自statsmodels的数据集 grunfeld ：“The Determinants of Corporate Investment”
这里我们试图使用企业的价值和资本来解释企业的投资行为。
from statsmodels.datasets import grunfeld
data = grunfeld.load_pandas().data
data.head()
| 0 |
317.6 |
3078.5 |
2.8 |
General Motors |
1935.0 |
| 1 |
391.8 |
4661.7 |
52.6 |
General Motors |
1936.0 |
| 2 |
410.6 |
5387.1 |
156.9 |
General Motors |
1937.0 |
| 3 |
257.7 |
2792.2 |
209.2 |
General Motors |
1938.0 |
| 4 |
330.8 |
4313.2 |
203.4 |
General Motors |
1939.0 |
变量的定义如下：
- invest - Gross investment in 1947 dollars
- value - Market value as of Dec. 31 in 1947 dollars
- capital - Stock of plant and equipment in 1947 dollars
- firm - General Motors, US Steel, General Electric, Chrysler, Atlantic Refining, IBM, Union Oil, Westinghouse, Goodyear, Diamond Match, American Steel
- year - 1935 - 1954
首先把表示个体和时间的列（准备要控制的固定效应）都转为index
data = data.set_index(["firm", "year"])
data.head()
| firm |
year |
|
|
|
| General Motors |
1935.0 |
317.6 |
3078.5 |
2.8 |
| 1936.0 |
391.8 |
4661.7 |
52.6 |
| 1937.0 |
410.6 |
5387.1 |
156.9 |
| 1938.0 |
257.7 |
2792.2 |
209.2 |
| 1939.0 |
330.8 |
4313.2 |
203.4 |
进行回归
考虑企业投资的决定，我们设置2个回归模型：
不带有固定效应: \[
invest_{it} = \beta_0 + \beta_1 value_{it} + \beta_2 capital_{it} + \varepsilon_{it}
\]
带有个体和时间固定效应，其中\(\nu_i\)、\(\lambda_t\)分别代表个体和时间固定效应： \[
invest_{it} = \beta_0 + \beta_1 value_{it} + \beta_2 capital_{it} + \nu_i + \lambda_t + \varepsilon_{it}
\]
这里 \(\nu_i\) 可以理解为企业 \(i\) 的“固定特征”，在所有年份都相同；\(\lambda_t\) 是年份 \(t\) 的共同冲击，对所有企业都相同。例如，宏观经济环境的变化就可以通过年份固定效应来近似刻画。
要控制固定效应，只要在公式中添加 EntityEffects或TimeEffects 即可。拟合两个模型：
from linearmodels import PanelOLS
# 无固定效应
m1 = PanelOLS.from_formula(
"invest ~ 1 + value + capital", data=data
).fit()
# 有固定效应
m2 = PanelOLS.from_formula(
"invest ~ 1 + value + capital + EntityEffects + TimeEffects", data=data
).fit()
# 有固定效应，聚类到个体
m3 = PanelOLS.from_formula(
"invest ~ 1 + value + capital + EntityEffects + TimeEffects", data=data
).fit(cov_type="clustered", cluster_entity=True)
比较结果：
from linearmodels.panel.results import compare
print(compare(results = {'OLS':m1,"FE":m2, 'FE-cluster':m3}))
Model Comparison
===================================================================
OLS FE FE-cluster
-------------------------------------------------------------------
Dep. Variable invest invest invest
Estimator PanelOLS PanelOLS PanelOLS
No. Observations 220 220 220
Cov. Est. Unadjusted Unadjusted Clustered
R-squared 0.8179 0.7253 0.7253
R-Squared (Within) 0.7357 0.7566 0.7566
R-Squared (Between) 0.8426 0.7944 0.7944
R-Squared (Overall) 0.8179 0.7856 0.7856
F-statistic 487.28 248.15 248.15
P-value (F-stat) 0.0000 0.0000 0.0000
===================== ============ ============ ===========
Intercept -38.410 -72.394 -72.394
(-4.5654) (-5.6861) (-3.5000)
value 0.1145 0.1167 0.1167
(20.753) (9.0219) (10.368)
capital 0.2275 0.3514 0.3514
(9.3904) (16.696) (7.4836)
======================= ============== ============== =============
Effects Entity Entity
Time Time
-------------------------------------------------------------------
T-stats reported in parentheses
注意到
- 2个“FE”列的最下方，显示这2个模型已经控制了个体和时间固定效应。
- Cov. Est.行， 显示FE-Cluster模型使用了聚类标准误。
