Contents

Семинар

I. Мультиколлинеарность.

Строгая мультиколлинераность появляется при нарушении условия ТГМ о том, что матрица \(X\) является матрицей полного ранга и возникает при наличии линейной зависимости между регрессорами. Однако на практике проблемы возникают и при не строгой мультиколлинеарности (квазимультиколлинеарности), при которой \(\det(X'X) \approx 0\).

В теории: строгая мультиколлинеарность приведёт к тому, что оценки МНК будут получаться неединственными, нестрогая мультиколлинеарность не нарушит ТГМ.

На практике: мультиколлинеарность приводит к тому, что:

  • оценки \(\hat{\beta}\) становятся неустойчивыми

  • оценки коэффициентов имеют высокие стандартные ошибки

  • доверительные интервалы становятся широкими

  • коэффициенты меняют знак или становятся незначимыми при небольших изменениях данных

  • t-статистики отдельных коэффициентов малы, но F-статистика на значимость в целом велика.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.formula.api as smf

diamonds = sns.load_dataset('diamonds')

# Создадим строгую мультиколлинеарность
X = diamonds[['price', 'x', 'y']]
X['sum_xy'] = X['x'] + X['y']
model1 = smf.ols(data = X, formula='price ~ x + y + sum_xy').fit()

# Слегка изменим данные
X1 = X.copy()
X1['x'] = X1['x'] + 0.001
model2 = smf.ols(data = X1, formula='price ~ x + y + sum_xy').fit()

# Посмотрим на разницу между коэффициентами
print(model1.params)
print(model2.params)
print('Суммарная квадратичная разность: %f' % (np.sum(np.square(model1.params - model2.params))))

<ipython-input-80-5cb0f7e5d949>:3: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  X['sum_xy'] = X['x'] + X['y']

Intercept   -14105.076664
x             1863.157850
y             -814.091949
sum_xy        1049.065901
dtype: float64
Intercept   -2.326418e+11
x            2.326418e+14
y            2.326418e+14
sum_xy      -2.326418e+14
dtype: float64
Суммарная квадратичная разность: 162366708248113373337532497920.000000

# Для сравнения: модель без мультиколлинеарности
Z = diamonds[['price', 'x', 'y']]
model3 = smf.ols(data = Z, formula='price ~ x + y').fit()

# Слегка изменим данные
Z1 = Z.copy()
Z1['x'] = Z1['x'] + 0.001
model4= smf.ols(data = Z1, formula='price ~ x + y').fit()

# Посмотрим на разницу между коэффициентами
print(model3.params)
print(model4.params)
print('Суммарная квадратичная разность: %f' % (np.sum(np.square(model3.params - model4.params))))

Intercept   -14105.076664
x             2912.223752
y              234.973952
dtype: float64
Intercept   -14107.988887
x             2912.223752
y              234.973952
dtype: float64
Суммарная квадратичная разность: 8.481047

# Исследуем ковариационную матрицу коэффициентов
model2.summary()
OLS Regression Results
Dep. Variable: price R-squared: 0.782
Model: OLS Adj. R-squared: 0.782
Method: Least Squares F-statistic: 6.462e+04
Date: Mon, 16 Nov 2020 Prob (F-statistic): 0.00
Time: 22:59:20 Log-Likelihood: -4.8265e+05
No. Observations: 53940 AIC: 9.653e+05
Df Residuals: 53936 BIC: 9.653e+05
Df Model: 3
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept -2.326e+11 2.87e+11 -0.811 0.418 -7.95e+11 3.3e+11
x 2.326e+14 2.87e+14 0.811 0.418 -3.3e+14 7.95e+14
y 2.326e+14 2.87e+14 0.811 0.418 -3.3e+14 7.95e+14
sum_xy -2.326e+14 2.87e+14 -0.811 0.418 -7.95e+14 3.3e+14
Omnibus: 18702.724 Durbin-Watson: 0.442
Prob(Omnibus): 0.000 Jarque-Bera (JB): 143426.007
Skew: 1.464 Prob(JB): 0.00
Kurtosis: 10.433 Cond. No. 8.90e+14


Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The smallest eigenvalue is 1.4e-23. This might indicate that there are
strong multicollinearity problems or that the design matrix is singular.

Индикаторы мультиколлинеарности.

  1. Очевидный: посмотреть на \(\mathrm{scorr}(X)\):

np.corrcoef(X1.T)

array([[1.        , 0.88443516, 0.8654209 , 0.88042812],
       [0.88443516, 1.        , 0.97470148, 0.99354016],
       [0.8654209 , 0.97470148, 1.        , 0.99376929],
       [0.88042812, 0.99354016, 0.99376929, 1.        ]])

  1. Condition number.

Напоминание: \(\|\cdot\|\) – это норма матрицы, если

  1. \(\|A\| \geq 0\), и если \(A = 0\), то \(\|A\| = 0\).

  2. \(\|a A\| = |a| \|A\|\)

  3. \(\|A+B\| \leq \|A\| + \|B\|\)

Некоторые матричные нормы также удовлетворяют свойству субмультипликативности: \(\Vert A B \Vert \leq \Vert A \Vert \Vert B \Vert\).

Condtition number определяется как \(\mathrm{K}(A) = \Vert A \Vert \Vert A^{-1} \Vert\). Можно показать, что в случае \(l_2\)-нормы для симметричной положительно определённой матрицы \(A\) верно:

\[ K(A) = \frac{\lambda_{\max}(A)}{\lambda_{\min}(A)} \]
np.linalg.cond(X.T @ X)

8.396003693050826e+18

  1. Variance Inflation Factor (VIF).

\[ \mathrm{VIF(X_k)} = \dfrac{1}{1 - R^2_j}, \]

где \(R^2_j\) – коэффициент детерминации в регрессии \(X_j\) на остальные факторы.

Если VIF превышает некоторое пороговое значение (5/10/12\(\ldots\)), то \(X_j\) сильно коллинеарен прочим регрессорам.

from statsmodels.stats.outliers_influence import variance_inflation_factor

variance_inflation_factor(model.model.exog, 1)

inf

variance_inflation_factor(model.model.exog, 2)

inf

variance_inflation_factor(model.model.exog, 3)

inf

Борьба с мультиколлинеарностью.

  1. Ничего не делать: оценки остаются несмещёнными и эффективными среди линейных по \(y\) несмещённых оценок.

  2. Выбросить часть регрессоров. В этом случае жертвуем несмещённостью оставшихся.

  3. Регуляризация. В этом случае жертвуем несмещённостью коэффициентов.

  4. Иногда – PCA.

  5. Получить больше наблюдений.

II. Отбор регрессоров.

  1. Метод последовательного исключения:

    1. Оценить модель.

    2. Выбросить наиболее незначимый регрессор.

    3. Переоценить модель.

    4. Повторять шаги B-C до тех пор, пока все регрессоры не станут значимыми.

  2. Метод последовательного включения:

    1. Оценить регрессию \(y\) на каждый из \(X_1\), \(\ldots\), \(X_k\). Выбрать самый значимый: \(X_q\).

    2. Оценить регерссию \(y\) на \(X_q\) и на каждый из \(X_1\), \(\ldots\), \(X_{q-1}\), \(X_{q+1}\), \(\ldots\), \(X_k\). Выбрать сымый значимый.

    3. Повторять до тех пор, пока среди включаемых коэффициентов не останется незначимых.

  1. Если есть \(k\) регрессоров и хотим узнать, нужно ли включать все или можно включить \(m<k\) регрессоров: F-тест (см. позапрошлый семинар).

  1. Если есть \(k\) регрессоров и хотим узнать, не нужно ли было включить ещё какие-то переменные, которых у нас нет: RESET-тест Рамсея.

\[\begin{split} \begin{cases} H_0: y = X\beta + u, \\ H_1: \text{есть пропущенные переменные}, \end{cases} \end{split}\]

Схема:

  1. Оценить регрессию \(\hat{y} = X\hat{\beta}\).

  2. Оценить регрессию \(y_i = \beta_0 + \beta_1X_1 + \ldots + \beta_kX_k + a_1\hat{y}_i^2 + \ldots + a_m\hat{y}_i^m + u_i\).

  3. Проверить гипотезу $\( \begin{cases} H_0: \forall i : a_i = 0, \\ H_1: \sum_i a_i^2 > 0, \end{cases} \)$

например, при помощи F-теста.

from statsmodels.stats.diagnostic import linear_reset

model = smf.ols('np.log(price) ~ depth + table + x + y + z', data = diamonds).fit()

linear_reset(model, use_f = True)
# H0 отвергается на 5% уровне значимости!

<class 'statsmodels.stats.contrast.ContrastResults'>
<F test: F=array([[7697.00812222]]), p=0.0, df_denom=5.39e+04, df_num=2>
Семинар