Python 中面板数据的多重共线性
Multicollinearity in Panel Data in Python
我习惯使用 Stata 或 R 来做线性回归模型，但我 运行将更多的工作流程分配给 Python。
这两个程序的有用之处在于它们直观地知道您不关心线性模型中的所有实体或时间固定效应，因此在估计面板模型时，它们会从中删除多重共线性虚拟变量模型（报告他们丢弃了哪些模型）。
虽然我知道以这种方式估计模型并不理想，并且应该小心回归到 运行（等），但这在实践中很有用，因为这意味着您可以先看到结果，然后担心假人的一些细微差别（特别是因为你不关心完全饱和的固定效应模型中的假人）。
让我举个例子。以下需要 linearmodels 并加载数据集并尝试 运行 面板回归。它是他们 documentation.
示例的修改版本
# Load the data (requires statsmodels and linearmodels)
import statsmodels.api as sm
from linearmodels.datasets import wage_panel
import pandas as pd
data = wage_panel.load()
year = pd.Categorical(data.year)
data = data.set_index(['nr', 'year'])
data['year'] = year
print(wage_panel.DESCR)
print(data.head())
# Run the panel regression
from linearmodels.panel import PanelOLS
exog_vars = ['exper','union','married']
exog = sm.add_constant(data[exog_vars])
mod = PanelOLS(data.lwage, exog, entity_effects=True, time_effects=True)
fe_te_res = mod.fit()
print(fe_te_res)
这会产生以下错误：
AbsorbingEffectError:
The model cannot be estimated. The included effects have fully absorbed
one or more of the variables. This occurs when one or more of the dependent
variable is perfectly explained using the effects included in the model.
但是，如果您通过将相同的数据导出到 Stata 在 Stata 中进行估算，运行宁:
data.drop(columns='year').to_stata('data.dta')
然后 运行在你的 stata 文件中设置等效项（加载数据后）：
xtset nr year
xtreg lwage exper union married i.year, fe
这将执行以下操作：
> . xtreg lwage exper union married i.year, fe
note: 1987.year omitted because of collinearity
Fixed-effects (within) regression Number of obs = 4360
Group variable: nr Number of groups = 545
R-sq: within = 0.1689 Obs per group: min = 8
between = 0.0000 avg = 8.0
overall = 0.0486 max = 8
F(9,3806) = 85.95
corr(u_i, Xb) = -0.1747 Prob > F = 0.0000
------------------------------------------------------------------------------
lwage | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
exper | .0638624 .0032594 19.59 0.000 .0574721 .0702527
union | .0833697 .0194393 4.29 0.000 .0452572 .1214821
married | .0583372 .0183688 3.18 0.002 .0223235 .0943509
|
year |
1981 | .0496865 .0200714 2.48 0.013 .0103348 .0890382
1982 | .0399445 .019123 2.09 0.037 .0024521 .0774369
1983 | .0193513 .018662 1.04 0.300 -.0172373 .0559398
1984 | .0229574 .0186503 1.23 0.218 -.0136081 .0595229
1985 | .0081499 .0191359 0.43 0.670 -.0293677 .0456674
1986 | .0036329 .0200851 0.18 0.856 -.0357456 .0430115
1987 | 0 (omitted)
|
_cons | 1.169184 .0231221 50.57 0.000 1.123851 1.214517
-------------+----------------------------------------------------------------
sigma_u | .40761229
sigma_e | .35343397
rho | .57083029 (fraction of variance due to u_i)
------------------------------------------------------------------------------
请注意，stata 从回归中任意删除了 1987，但仍然 运行。有没有办法在 linearmodels 或 statsmodels 中获得类似的功能？
我能想到的唯一方法是手动。
示例数据
import pandas as pd
import numpy as np
import statsmodels.api as sm
from linearmodels.datasets import wage_panel
from linearmodels.panel import PanelOLS
data = wage_panel.load()
首先，我们将跟随 Stata 的脚步，为每年的固定效应生成虚拟变量，并省略第一个值，按字典顺序排序（通过 drop_first=True 参数完成）。使用 np.unique 然后获取标签很重要，因为这也是一种排序。不用statsmodels加常量，自己动手就行：
data = wage_panel.load()
data = pd.concat([data, pd.get_dummies(data.year, drop_first=True)], axis=1)
exog_vars = ['exper','union','married'] + np.unique(data.year)[1::].tolist()
data = data.set_index(['nr', 'year'])
exog = data[exog_vars].assign(constant=1)
如果我们尝试 运行 这个模型，它会因为完全共线性而失败。因为我们正在做 within-regression 我们不能只测试 exog 的共线性，我们需要首先在面板内去均值，因为去均值的自变量的共线性才是最重要的。我会在这里复制一份，以免弄乱我们将在最终回归中使用的 exog。
exog2 = exog.copy()
exog2 = exog2 - exog2.groupby(level=0).transform('mean') + exog2.mean()
我们现在可以看到存在完全共线性；当我们回归 within de-meaned exog 变量对所有其他变量时，我们得到一些回归的完美 R 平方：
for col in exog2:
print(col, sm.OLS(exog2[col], exog2.drop(columns=col)).fit().rsquared)
#exper 1.0
#union 0.004326532427527674
#married 0.18901677578724896
#1981 1.0
#1982 1.0
#1983 1.0
#1984 1.0
#1985 1.0
#1986 1.0
#1987 1.0
现在 Stata 或其他软件包如何决定删除哪个变量对我来说是个谜。在这种情况下，您可能更愿意放弃其中一个年度虚拟变量而不是其他变量。我们就选择1987吧，这样我们最终可以得到和Stata一样的结果。
exog2 = exog2.drop(columns=1987)
for col in exog2:
print(col, sm.OLS(exog2[col], exog2.drop(columns=col)).fit().rsquared)
#exper 0.48631
#union 0.0043265
#married 0.1890
#1981 0.34978
#1982 0.28369
#1983 0.2478680
#1984 0.2469180
#1985 0.2846552
#1986 0.35067075
#constant 0.94641
所以我们已经处理了共线性并且可以回到模型。由于我们手动包含了年度固定效应，因此我们可以从模型中删除 time_effects。
mod = PanelOLS(data.lwage, exog.drop(columns=1987), entity_effects=True, time_effects=False)
print(mod.fit())
PanelOLS Estimation Summary
================================================================================
Dep. Variable: lwage R-squared: 0.1689
Estimator: PanelOLS R-squared (Between): -0.0882
No. Observations: 4360 R-squared (Within): 0.1689
Date: Sat, Mar 09 2019 R-squared (Overall): 0.0308
Time: 00:59:14 Log-likelihood -1355.7
Cov. Estimator: Unadjusted
F-statistic: 85.946
Entities: 545 P-value 0.0000
Avg Obs: 8.0000 Distribution: F(9,3806)
Min Obs: 8.0000
Max Obs: 8.0000 F-statistic (robust): 85.946
P-value 0.0000
Time periods: 8 Distribution: F(9,3806)
Avg Obs: 545.00
Min Obs: 545.00
Max Obs: 545.00
Parameter Estimates
==============================================================================
Parameter Std. Err. T-stat P-value Lower CI Upper CI
------------------------------------------------------------------------------
exper 0.0639 0.0033 19.593 0.0000 0.0575 0.0703
union 0.0834 0.0194 4.2887 0.0000 0.0453 0.1215
married 0.0583 0.0184 3.1759 0.0015 0.0223 0.0944
1981 0.0497 0.0201 2.4755 0.0133 0.0103 0.0890
1982 0.0399 0.0191 2.0888 0.0368 0.0025 0.0774
1983 0.0194 0.0187 1.0369 0.2998 -0.0172 0.0559
1984 0.0230 0.0187 1.2309 0.2184 -0.0136 0.0595
1985 0.0081 0.0191 0.4259 0.6702 -0.0294 0.0457
1986 0.0036 0.0201 0.1809 0.8565 -0.0357 0.0430
constant 1.1692 0.0231 50.566 0.0000 1.1239 1.2145
==============================================================================
这与 Stata 输出匹配。确实没有任何理由放弃 1987，我们可以选择其他内容，但至少我们可以看到结果与 xtreg 匹配。
我习惯使用 Stata 或 R 来做线性回归模型，但我 运行将更多的工作流程分配给 Python。
这两个程序的有用之处在于它们直观地知道您不关心线性模型中的所有实体或时间固定效应，因此在估计面板模型时，它们会从中删除多重共线性虚拟变量模型（报告他们丢弃了哪些模型）。
虽然我知道以这种方式估计模型并不理想，并且应该小心回归到 运行（等），但这在实践中很有用，因为这意味着您可以先看到结果，然后担心假人的一些细微差别（特别是因为你不关心完全饱和的固定效应模型中的假人）。
让我举个例子。以下需要 linearmodels 并加载数据集并尝试 运行 面板回归。它是他们 documentation.
# Load the data (requires statsmodels and linearmodels)
import statsmodels.api as sm
from linearmodels.datasets import wage_panel
import pandas as pd
data = wage_panel.load()
year = pd.Categorical(data.year)
data = data.set_index(['nr', 'year'])
data['year'] = year
print(wage_panel.DESCR)
print(data.head())
# Run the panel regression
from linearmodels.panel import PanelOLS
exog_vars = ['exper','union','married']
exog = sm.add_constant(data[exog_vars])
mod = PanelOLS(data.lwage, exog, entity_effects=True, time_effects=True)
fe_te_res = mod.fit()
print(fe_te_res)
这会产生以下错误：
AbsorbingEffectError: The model cannot be estimated. The included effects have fully absorbed one or more of the variables. This occurs when one or more of the dependent variable is perfectly explained using the effects included in the model.
但是，如果您通过将相同的数据导出到 Stata 在 Stata 中进行估算，运行宁:
data.drop(columns='year').to_stata('data.dta')
然后 运行在你的 stata 文件中设置等效项（加载数据后）：
xtset nr year
xtreg lwage exper union married i.year, fe
这将执行以下操作：
> . xtreg lwage exper union married i.year, fe
note: 1987.year omitted because of collinearity
Fixed-effects (within) regression Number of obs = 4360
Group variable: nr Number of groups = 545
R-sq: within = 0.1689 Obs per group: min = 8
between = 0.0000 avg = 8.0
overall = 0.0486 max = 8
F(9,3806) = 85.95
corr(u_i, Xb) = -0.1747 Prob > F = 0.0000
------------------------------------------------------------------------------
lwage | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
exper | .0638624 .0032594 19.59 0.000 .0574721 .0702527
union | .0833697 .0194393 4.29 0.000 .0452572 .1214821
married | .0583372 .0183688 3.18 0.002 .0223235 .0943509
|
year |
1981 | .0496865 .0200714 2.48 0.013 .0103348 .0890382
1982 | .0399445 .019123 2.09 0.037 .0024521 .0774369
1983 | .0193513 .018662 1.04 0.300 -.0172373 .0559398
1984 | .0229574 .0186503 1.23 0.218 -.0136081 .0595229
1985 | .0081499 .0191359 0.43 0.670 -.0293677 .0456674
1986 | .0036329 .0200851 0.18 0.856 -.0357456 .0430115
1987 | 0 (omitted)
|
_cons | 1.169184 .0231221 50.57 0.000 1.123851 1.214517
-------------+----------------------------------------------------------------
sigma_u | .40761229
sigma_e | .35343397
rho | .57083029 (fraction of variance due to u_i)
------------------------------------------------------------------------------
请注意，stata 从回归中任意删除了 1987，但仍然 运行。有没有办法在 linearmodels 或 statsmodels 中获得类似的功能？
我能想到的唯一方法是手动。
示例数据
import pandas as pd
import numpy as np
import statsmodels.api as sm
from linearmodels.datasets import wage_panel
from linearmodels.panel import PanelOLS
data = wage_panel.load()
首先，我们将跟随 Stata 的脚步，为每年的固定效应生成虚拟变量，并省略第一个值，按字典顺序排序（通过 drop_first=True 参数完成）。使用 np.unique 然后获取标签很重要，因为这也是一种排序。不用statsmodels加常量，自己动手就行：
data = wage_panel.load()
data = pd.concat([data, pd.get_dummies(data.year, drop_first=True)], axis=1)
exog_vars = ['exper','union','married'] + np.unique(data.year)[1::].tolist()
data = data.set_index(['nr', 'year'])
exog = data[exog_vars].assign(constant=1)
如果我们尝试 运行 这个模型，它会因为完全共线性而失败。因为我们正在做 within-regression 我们不能只测试 exog 的共线性，我们需要首先在面板内去均值，因为去均值的自变量的共线性才是最重要的。我会在这里复制一份，以免弄乱我们将在最终回归中使用的 exog。
exog2 = exog.copy()
exog2 = exog2 - exog2.groupby(level=0).transform('mean') + exog2.mean()
我们现在可以看到存在完全共线性；当我们回归 within de-meaned exog 变量对所有其他变量时，我们得到一些回归的完美 R 平方：
for col in exog2:
print(col, sm.OLS(exog2[col], exog2.drop(columns=col)).fit().rsquared)
#exper 1.0
#union 0.004326532427527674
#married 0.18901677578724896
#1981 1.0
#1982 1.0
#1983 1.0
#1984 1.0
#1985 1.0
#1986 1.0
#1987 1.0
现在 Stata 或其他软件包如何决定删除哪个变量对我来说是个谜。在这种情况下，您可能更愿意放弃其中一个年度虚拟变量而不是其他变量。我们就选择1987吧，这样我们最终可以得到和Stata一样的结果。
exog2 = exog2.drop(columns=1987)
for col in exog2:
print(col, sm.OLS(exog2[col], exog2.drop(columns=col)).fit().rsquared)
#exper 0.48631
#union 0.0043265
#married 0.1890
#1981 0.34978
#1982 0.28369
#1983 0.2478680
#1984 0.2469180
#1985 0.2846552
#1986 0.35067075
#constant 0.94641
所以我们已经处理了共线性并且可以回到模型。由于我们手动包含了年度固定效应，因此我们可以从模型中删除 time_effects。
mod = PanelOLS(data.lwage, exog.drop(columns=1987), entity_effects=True, time_effects=False)
print(mod.fit())
PanelOLS Estimation Summary
================================================================================
Dep. Variable: lwage R-squared: 0.1689
Estimator: PanelOLS R-squared (Between): -0.0882
No. Observations: 4360 R-squared (Within): 0.1689
Date: Sat, Mar 09 2019 R-squared (Overall): 0.0308
Time: 00:59:14 Log-likelihood -1355.7
Cov. Estimator: Unadjusted
F-statistic: 85.946
Entities: 545 P-value 0.0000
Avg Obs: 8.0000 Distribution: F(9,3806)
Min Obs: 8.0000
Max Obs: 8.0000 F-statistic (robust): 85.946
P-value 0.0000
Time periods: 8 Distribution: F(9,3806)
Avg Obs: 545.00
Min Obs: 545.00
Max Obs: 545.00
Parameter Estimates
==============================================================================
Parameter Std. Err. T-stat P-value Lower CI Upper CI
------------------------------------------------------------------------------
exper 0.0639 0.0033 19.593 0.0000 0.0575 0.0703
union 0.0834 0.0194 4.2887 0.0000 0.0453 0.1215
married 0.0583 0.0184 3.1759 0.0015 0.0223 0.0944
1981 0.0497 0.0201 2.4755 0.0133 0.0103 0.0890
1982 0.0399 0.0191 2.0888 0.0368 0.0025 0.0774
1983 0.0194 0.0187 1.0369 0.2998 -0.0172 0.0559
1984 0.0230 0.0187 1.2309 0.2184 -0.0136 0.0595
1985 0.0081 0.0191 0.4259 0.6702 -0.0294 0.0457
1986 0.0036 0.0201 0.1809 0.8565 -0.0357 0.0430
constant 1.1692 0.0231 50.566 0.0000 1.1239 1.2145
==============================================================================
这与 Stata 输出匹配。确实没有任何理由放弃 1987，我们可以选择其他内容，但至少我们可以看到结果与 xtreg 匹配。