Мультиколлинеарность в панельных данных в Python - PullRequest
2 голосов
/ 09 марта 2019

Я привык использовать Stata или R для создания моделей линейной регрессии, но я перевожу больше рабочих процессов на Python.

Полезная вещь в этих двух программах заключается в том, что они интуитивно знают, что вам все равновсе связанные с сущностью или временем эффекты в линейной модели, поэтому при оценке панельных моделей они будут отбрасывать мультиколлинеарные манекены из модели (сообщая, какие из них отбрасываются).

Хотя я понимаю, что оценка моделей втакой способ не идеален, и нужно быть осторожным с регрессиями для запуска (и т. д.), это полезно на практике, потому что это означает, что вы можете сначала увидеть результаты, а потом беспокоиться о некоторых нюансах манекенов (особенно если выменя не волнуют пустышки в полностью насыщенной модели с фиксированными эффектами).

Позвольте привести пример.Следующее требует linearmodels и загружает набор данных и пытается запустить регрессию панели.Это модифицированная версия примера из их документации .

# 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: Модель не может быть оценена.Включенные эффекты полностью поглотили одну или несколько переменных.Это происходит, когда одна или несколько зависимых переменных полностью объяснены с использованием эффектов, включенных в модель.

Однако, если вы оцениваете в 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?

1 Ответ

2 голосов
/ 09 марта 2019

Единственный способ, которым я могу думать, это сделать вручную.

Образцы данных

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()

Сначала мы пойдем по стопам Статы, создавая макеты для каждого изгод фиксировал эффекты, и мы пропускаем первое значение, лексикографически отсортированное (выполняется с аргументом 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)

Если мы попытаемся запустить эту модель, она потерпит неудачу из-за идеальной коллинеарности.Поскольку мы делаем в рамках -регрессии, мы не можем просто проверить коллинеарность на экзоге, нам нужно сначала деинвестировать в панели, так как коллинеарность подразумеваемых независимых переменных имеет значение.Я сделаю копию здесь, чтобы не связываться с exog, который мы будем использовать в окончательной регрессии.

exog2 = exog.copy()
exog2 = exog2 - exog2.groupby(level=0).transform('mean') + exog2.mean()

Теперь мы можем видеть, что существует совершенная коллинеарность;когда мы регрессируем в пределах де-значения переменной 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.

...