快速提示：
- 经济学中的回归，基本目标是解释2个变量之间的关系：x如何影响y？
- 这里更重视解释变量的系数 \(\beta_i\) 及其显著性，而不是太关心 \(R^2\) 等统计量
- 因此，回归模型的设计，以及对系数的进行统计检验以及解释，是经济学中回归的重点。
基本多元回归
分析目标和数据
这个案例来自伍德里奇的《计量经济学导论》。
本案希望研究的问题是：教育水平对工资有什么影响？即：每多读一年书，工资会有什么变化？
数据采用《导论》中的数据Wage1。
# 加载并查看wage1.csv数据集
import pandas as pd
# 加载数据
wage1_data = pd.read_csv("data/wage1.csv")
# 查看数据集的前几行和列名
wage1_data.head(), wage1_data.columns
( wage educ exper tenure nonwhite female married numdep smsa \
0 3.10 11 2 0 0 1 0 2 1
1 3.24 12 22 2 0 1 1 3 1
2 3.00 11 2 0 0 0 0 2 0
3 6.00 8 44 28 0 0 1 0 1
4 5.30 12 7 2 0 0 1 1 0
northcen ... trcommpu trade services profserv profocc clerocc \
0 0 ... 0 0 0 0 0 0
1 0 ... 0 0 1 0 0 0
2 0 ... 0 1 0 0 0 0
3 0 ... 0 0 0 0 0 1
4 0 ... 0 0 0 0 0 0
servocc lwage expersq tenursq
0 0 1.131402 4 0
1 1 1.175573 484 4
2 0 1.098612 4 0
3 0 1.791759 1936 784
4 0 1.667707 49 4
[5 rows x 24 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'))
回归公式
我们建立的回归模型可以表达为：
\[
\text{wage} = \beta_0 + \delta_0 \cdot \text{female} + \beta_1 \cdot \text{educ} + \beta_2 \cdot \text{exper} + \beta_3 \cdot \text{tenure} + u
\]
其中：
- \(\text{wage}\) 是工资。
- \(\text{female}\) 是一个虚拟变量，女性为1，男性为0。
- \(\text{educ}\), \(\text{exper}\), 和 \(\text{tenure}\) 分别代表教育年数、工作经验和在当前工作的年限。
- \(u\) 是误差项。
教育水平对工资的影响体现在系数\(\beta_1\)中，这就是我们关心的目标。
进行回归分析
在进行回归分析之前，我们需要准备数据和定义模型。这个过程包括设置自变量（解释变量）和因变量（被解释变量），以及确保数据格式适合进行回归分析。
定义因变量和自变量
因变量 \(y\): 因变量是我们想要预测或解释的变量。在本例中，因变量是 wage，即工资，我们希望了解其他变量如何影响工资。
自变量 \(X\): 自变量是可能影响因变量的变量。在这个模型中，我们选择了四个自变量：female, educ, exper, 和 tenure。这些变量分别表示：
female: 表示性别，是一个二分类变量，女性为1，男性为0。educ: 表示受教育年数，是一个连续变量。exper: 表示总工作经验年数，是一个连续变量。tenure: 表示在当前工作的年限，也是一个连续变量。
构造自变量 \(X\) 和因变量 \(y\)
在开始构建模型之前，我们需要从数据集中提取这些变量。使用Pandas库，我们可以轻松地从DataFrame中选择需要的列。同时，为了进行回归分析，我们需要向自变量矩阵 \(X\) 中添加一个常数项，以便模型包含截距 \(\beta_0\)。
import statsmodels.api as sm
# 选择自变量和因变量
X = wage1_data[["female", "educ", "exper", "tenure"]]
y = wage1_data["wage"]
# 向自变量矩阵X添加常数项，以便包含截距
X = sm.add_constant(X)
拟合回归模型
使用statsmodels库的OLS类，我们可以定义并拟合一个普通最小二乘回归模型。OLS类的第一个参数是因变量 \(y\)，第二个参数是自变量 \(X\)。
# 创建OLS回归模型的实例并拟合数据
model = sm.OLS(y, X)
results = model.fit()
查看回归结果
模型拟合后，我们可以通过调用 summary() 方法来查看详细的回归分析结果，包括每个变量的系数、标准误、t 值、p 值等统计指标。
# 输出回归结果摘要
print(results.summary())
OLS Regression Results
==============================================================================
Dep. Variable: wage R-squared: 0.364
Model: OLS Adj. R-squared: 0.359
Method: Least Squares F-statistic: 74.40
Date: Thu, 12 Dec 2024 Prob (F-statistic): 7.30e-50
Time: 03:11:57 Log-Likelihood: -1314.2
No. Observations: 526 AIC: 2638.
Df Residuals: 521 BIC: 2660.
Df Model: 4
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const -1.5679 0.725 -2.164 0.031 -2.991 -0.145
female -1.8109 0.265 -6.838 0.000 -2.331 -1.291
educ 0.5715 0.049 11.584 0.000 0.475 0.668
exper 0.0254 0.012 2.195 0.029 0.003 0.048
tenure 0.1410 0.021 6.663 0.000 0.099 0.183
==============================================================================
Omnibus: 185.864 Durbin-Watson: 1.794
Prob(Omnibus): 0.000 Jarque-Bera (JB): 715.580
Skew: 1.589 Prob(JB): 4.11e-156
Kurtosis: 7.749 Cond. No. 141.
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
回归线
可以得到回归线（或者样本回归函数）：
\[
wage = -1.568 + (-1.811) \cdot female + 0.572 \cdot educ + 0.025 \cdot exper + 0.141 \cdot tenure + u
\]
解释回归结果
从回归结果中，我们可以提取以下关键信息：
- 模型拟合情况: \(R^2\) 值为 0.364，说明模型解释了 36.4% 的因变量（wage）的方差。这表明模型有一定的解释能力，但仍有大部分变异未被模型捕捉。
- 系数解释:
- const（截距）: \(-1.568\), 在统计上显著（\(p < 0.05\)），表示当所有解释变量为0时，预测的工资值。
- female: \(-1.811\), 显著地负相关（\(p < 0.001\)），表明女性的工资平均比男性低约 $1.81，条件是其他变量保持不变。
- educ（教育年数）: \(0.572\), 显著正相关（\(p < 0.001\)），意味着每增加一年的教育，工资平均增加约 $0.57，其他条件不变。
- exper（工作经验）: \(0.025\), 在边际上显著（\(p = 0.029\)），表示每增加一年工作经验，工资增加约 $0.03，其他条件不变。
- tenure（在当前工作的年限）: \(0.141\), 显著正相关（\(p < 0.001\)），每增加一年在当前工作的年限，工资增加约 $0.14，其他条件不变。
- 模型诊断:
- Durbin-Watson 统计量为 1.794，接近2，说明残差中不存在自相关问题。
- Omnibus 和 Jarque-Bera 测试的 \(p\) 值都接近0，表明残差不符合正态分布，这可能影响了某些统计测试的有效性，尤其是在小样本数据中。
回归结果的经济学解释
在经济学研究中，分析者往往特别关注模型中的系数（\(\beta_i\)），因为这些系数能够揭示不同变量之间的关系，如大小、方向和显著性。这些关系有助于解释经济理论中的假设或验证特定的经济行为。每个系数的解释性质能够表达自变量对因变量的预期变化量，即当其他条件保持不变时，自变量每变化一个单位，因变量平均预期会如何变化。
对系数的解释
大小：系数的大小告诉我们，当自变量增加一个单位时，因变量预期会增加或减少多少。在经济学中，了解这种影响的规模是至关重要的，因为它有助于量化政策变动或市场条件变化对经济行为的影响。
方向：系数的符号（正或负）表明变量间的关系是正向还是负向。例如，如果某个变量的系数为正，则表明这两个变量是正相关的；如果系数为负，则两者是负相关的。
显著性：统计显著性（通常通过p值来衡量）告诉我们系数是否在统计上显著，即我们观察到的关系是否有可能仅仅是由于抽样误差。p值较小（通常小于0.05）意味着我们有足够的证据认为该系数在总体中不为零。
以前面的回归分析为例，我们研究了教育水平（educ）对工资（wage）的影响。根据回归结果，教育的系数是0.572，且非常显著（p值 < 0.001）。
这就回答了我们最初的问题：在其他因素不变的情况下，每增加一年的教育，工资平均预期增加0.572单位。
多元回归练习题
在一项研究中，研究人员收集了一批学生的高中成绩（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和ACT成绩的系数，包括意义和显著性。
- 思考题：讨论模型的适用性和限制。
# 读取数据
gpa1_data = pd.read_csv("data/gpa1.csv")
# 显示列名
gpa1_data.columns
Index(['age', 'soph', 'junior', 'senior', 'senior5', 'male', 'campus',
'business', 'engineer', 'colGPA', 'hsGPA', 'ACT', 'job19', 'job20',
'drive', 'bike', 'walk', 'voluntr', 'PC', 'greek', 'car', 'siblings',
'bgfriend', 'clubs', 'skipped', 'alcohol', 'gradMI', 'fathcoll',
'mothcoll'],
dtype='object')
固定效应
更进一步，我们可以使用 linearmodels来完成基准回归，并且加入固定效应
安装：
在Jupyter Notebook的python单元格中执行如下代码，注意开头有个感叹号：
!pip install 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 |
首先把表示个体和时间的列（准备要控制的固定效应）都转为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}
\]
要控制固定效应，只要在公式中添加 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
注意到 1. 2个“FE”列的最下方，显示这2个模型已经控制了个体和时间固定效应。 2. Cov. Est.行， 显示FE-Cluster模型使用了聚类标准误。