Заметки по R: Мультиколлинеарность
24.03.2016
Загружаем нужные пакеты:
library("knitr") # создание отчётов
# opts_chunk$set(cache=TRUE)
library("ggplot2") # графики
library("glmnet") # LASSO
library("car") # compute vif()
library("MASS") # ridge regression
library("quantreg") # quantile regressionВозьмем встроенные в R данные по скорости автомобилей и длине тормозного пути. Это данные 1920-х годов. Скорость в милях в час, миля — это примерно \(1.61\) километра. Длина тормозного пути — в футах. Один фут — это примерно \(30\) сантиметров.
ggplot(data = cars, aes(x = speed, y = dist)) + geom_point() +
labs(x = "Скорость машины, миль в час",
y = "Длина тормозного пути, футов",
title = "Длина тормозного пути") 
Построим полиномиальную регрессию тормозного пути от скорости:
m.poly3 <- lm(dist ~ speed + I(speed^2) + I(speed^3), data = cars)Предвосхищаем вопрос: Что значит I(...)? Буква I со скобками заставляет R воспринимать выражение внутри скобок в естественном математическом смысле, а не в специальном “регрессионном” смысле. Например:
y ~ x + z— регрессия \(y\) на \(x\) и \(z\)y ~ I(x + z)— регрессия \(y\) на одну переменную \(x+z\)y ~ (x + z + w)^2— регрессия \(y\) на \(x\), \(z\), \(w\) и все попарные произведения, т.е. на \(xz\), \(xw\) и \(wz\).y ~ I((x + z + w)^2)— регрессия на одну переменную \((x+z+w)^2\)y ~ (x + z + w)^3— регрессия \(y\) на \(x\), \(z\), \(w\), \(xz\), \(xw\), \(zw\), \(xzw\)y ~ I((x + z + w)^3)— регрессия на одну переменную \((x+z+w)^3\)
Смотрим на отчёт по регрессии:
summary(m.poly3)##
## Call:
## lm(formula = dist ~ speed + I(speed^2) + I(speed^3), data = cars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.670 -9.601 -2.231 7.075 44.691
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -19.50505 28.40530 -0.687 0.496
## speed 6.80111 6.80113 1.000 0.323
## I(speed^2) -0.34966 0.49988 -0.699 0.488
## I(speed^3) 0.01025 0.01130 0.907 0.369
##
## Residual standard error: 15.2 on 46 degrees of freedom
## Multiple R-squared: 0.6732, Adjusted R-squared: 0.6519
## F-statistic: 31.58 on 3 and 46 DF, p-value: 3.074e-11Мы видим, что:
ни один коэффициент не значим
регрессия в целом значима
Это один из признаков мультиколлинеарности.
Мультиколлинеарность — ситуация, когда все предпосылки теоремы Гаусса-Маркова выполнены, но дисперсия оценок является слишком высокой по твоему мнению
Другие признаки мультиколлинеарности:
Нестабильность коэффициентов при удалении одного наблюдения;
Высокие корреляции регрессоров;
Показатели в духе VIF и CI.
Корреляционная матрица регрессоров:
X <- model.matrix(~ 0 + speed + I(speed^2) + I(speed^3), data = cars)
cor(X)## speed I(speed^2) I(speed^3)
## speed 1.0000000 0.9794765 0.9389237
## I(speed^2) 0.9794765 1.0000000 0.9884061
## I(speed^3) 0.9389237 0.9884061 1.0000000Расчитаем индекс обусловленности, condition index, \[ CI=\sqrt{\frac{\lambda_{max}}{\lambda_{min}}} \] где \(\lambda\) — собственное число матрицы \(X'X\).
Вроде бы в некоторых источниках считают собственные числа у корреляционной или ковариационной матрицы регрессоров (???) (Лена Вакуленко) (Некоторые источники не указаны несколько дней. Вы можете улучшить эту статью :)
X <- model.matrix( ~ speed + I(speed^2) + I(speed^3), data = cars)
XX <- t(X) %*% X
eigen <- eigen(XX)
eigen$values## [1] 2.080900e+09 1.586193e+05 7.315488e+01 2.719130e-01CI <- sqrt(max(eigen$values) / min(eigen$values))
CI## [1] 87480.38Наверняка есть пакет, в котором CI уже реализован. Но иногда дольше искать нужный пакет, чем написать функцию самому. Тем более тут нагляднее формула видна.
Википедия говорит, что \(CI>30\) признак мультиколлинеарности. Но это не строгая граница. Мультиколлинеарность — это субъективное недовольство высокими стандартными ошибками коэффициентов.
Рассчитаем коэффициент вздутия дисперсии, variance inflation factor \[ VIF_j=\frac{1}{1-R_j^2} \] где \(R_j^2\) — это коэффициент \(R^2\) в регрессии \(j\)-ой объясняющей переменной на остальные.
Популярная граница для VIF — 10. У нас:
vif(m.poly3)## speed I(speed^2) I(speed^3)
## 274.1131 1408.0884 483.0570NB: нужно уметь рассчитывать vif ручками, зная коэффициенты R^2 во вспомогательных регрессиях. Например:
m.vif1 <- lm(speed ~ I(speed^2) + I(speed^3), data = cars)
r2.1 <- summary(m.vif1)$r.squared
vif1 <- 1 / (1 - r2.1)
vif1## [1] 274.1131Забудем про имеющуюся мультиколлинеарностьи попробуем заняться прогнозированием…
Из модели можно вытащить ковариационную матрицу оценок коэффициентов
vcov(m.poly3)## (Intercept) speed I(speed^2) I(speed^3)
## (Intercept) 806.8612229 -186.19855321 12.875757330 -0.2736117808
## speed -186.1985532 46.25543452 -3.346993631 0.0733088858
## I(speed^2) 12.8757573 -3.34699363 0.249882783 -0.0055981567
## I(speed^3) -0.2736118 0.07330889 -0.005598157 0.0001276477Можно построить прогноз длины тормозного пути для машины с скоростью \(speed=10\) миль в час.
new <- data.frame(speed = 10)
predict(m.poly3, newdata=new, interval = 'confidence')## fit lwr upr
## 1 23.79228 15.8941 31.69045Если нужно вытащить стандартную ошибку для \(\hat{y}\), то можно внутрь команды predict добавить опцию se.fit=TRUE.
На всякий случай, есть два типа доверительных интервалов при прогнозировании. Доверительный интервал для мат. ожидания будущего игрека — опция 'confidence' в R. И прогнозный интервал для индивидуального значения игрека — опция 'prediction' в R.
Попробуем перейти к центированным переменным.
cars$spc <- cars$speed - mean(cars$speed)
m.poly3c <- lm(dist ~ spc + I(spc^2) + I(spc^3), data = cars)
summary(m.poly3c)##
## Call:
## lm(formula = dist ~ spc + I(spc^2) + I(spc^3), data = cars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.670 -9.601 -2.231 7.075 44.691
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 39.75032 2.86281 13.885 < 2e-16 ***
## spc 3.32577 0.84227 3.949 0.000268 ***
## I(spc^2) 0.12399 0.07120 1.741 0.088298 .
## I(spc^3) 0.01025 0.01130 0.907 0.368919
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.2 on 46 degrees of freedom
## Multiple R-squared: 0.6732, Adjusted R-squared: 0.6519
## F-statistic: 31.58 on 3 and 46 DF, p-value: 3.074e-11Вау! Теперь появились значимые коэффициенты!!!!!
Но это всего лишь иллюзия. У этих оценок коэффициентов другая ковариационная матрица.
vcov(m.poly3c)## (Intercept) spc I(spc^2) I(spc^3)
## (Intercept) 8.195662560 0.325738550 -0.134011045 -0.0061108467
## spc 0.325738550 0.709415046 -0.016836556 -0.0082955869
## I(spc^2) -0.134011045 -0.016836556 0.005069365 0.0002991650
## I(spc^3) -0.006110847 -0.008295587 0.000299165 0.0001276477И если построить доверительный интервал для прогнозов, то ничего не изменится.
new$spc <- 10 - mean(cars$speed)
predict(m.poly3c, newdata = new, interval = 'confidence')## fit lwr upr
## 1 23.79228 15.8941 31.69045Модель с центрированными переменными — это всего лишь по-другому записанная исходная модель.
Кстати, нужно уметь строить доверительный интервал для прогноза “ручками” имея на руках оценки коэффициентов и оценку ковариационной матрицы.
Как же бороться с мультиколлинеарностью?
- Смириться с высокими стандартными ошибками. Высокие? Ну и чёрт с ними!
- Увеличить количество наблюдений.
- Задать другой вопрос.
“I checked it very thoroughly,” said the computer, “and that quite definitely is the answer. I think the problem, to be quite honest with you, is that you’ve never actually known what the question is.”
Douglas Adams, Hitchhiker’s Guide to the Galaxy.
Обычная регрессия отвечает на такой вопрос:
Чему равны несмещенные оценки с наименьшей дисперсией в модели \(E(y)=\beta_1+\beta_2 x +\beta_3 z\)?
Можно задаться целью получить смещенные оценки. Почему нам интересны смещенные оценки? Дао учит нас, что любая смещенная оценка коэффициента \(\beta\) есть несмещенная оценка какого-то коэффициента \(\gamma\).
Поэтому можно:
- Убрать из регрессии некоторые переменные
- В частности, можно оставить в регрессии “главные компоненты”, искусственные переменные, которые содержат большую часть информации обо всех исходных переменных
- При минимизации RSS ввести штраф за слишком большие коэффициенты \(\hat{\beta}\). Два популярных способа ввести штраф — Ridge Regression и LASSO
- Выяснить как зависит медиана (или другой квантиль \(y\)) от объясняющих переменных
Убрать из регрессии некоторые переменные
В регрессии \(y\) на \(x\) и \(z\) коэффициент при \(x\) показывает насколько изменится \(y\) при изменении \(x\) на единичку и фиксированном \(z\). В регрессии \(y\) на \(x\) коэффициент при \(x\) показывает насколько изменится \(y\) при изменении \(x\) на единичку и нефиксированном \(z\).
m.line <- lm(dist ~ speed, data = cars)
summary(m.line)##
## Call:
## lm(formula = dist ~ speed, data = cars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -29.069 -9.525 -2.272 9.215 43.201
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -17.5791 6.7584 -2.601 0.0123 *
## speed 3.9324 0.4155 9.464 1.49e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.38 on 48 degrees of freedom
## Multiple R-squared: 0.6511, Adjusted R-squared: 0.6438
## F-statistic: 89.57 on 1 and 48 DF, p-value: 1.49e-12Метод главных компонент
В частности можно попытаться оставить только несколько переменных, которые “впитывают” в себя большую часть изменчивости всех исходных переменных. Этот метод называется метод главных компонент. Подробнее отдельно.
Ridge regression
Теоретически RR позволяет при удаче получить оценки \(\hat{\beta}\) с меньшей, чем у OLS величиной среднеквадратичной ошибкой, MSE.
Алгоритм Ridge regression минимизирует функцию \[ \sum_{i=1}^n (y_i-\hat{y}_i)^2+\lambda \sum_{j=1}^k \hat{\beta}_j^2 \]
По сравнению с обычной регрессией введен “штраф” \(\lambda\) за большие значения \(\hat{\beta}_j\).
Эту задачу можно решить в явном виде: \[ \hat{\beta}_{RR}=(X'X+\lambda I)^{-1}X'y$ \]
По другому RR можно сформулировать так: минимизируем обычный RSS при ограничении \(\sum_{j=1}^k \hat{\beta}_j^2 \leq c\). Выписав функцию Лагранжа для этой задачи минимизации получим формулировку задачи RR с лямбдой.
Делаем RR в R :)
lambdas <- seq(0.1, 10, by = 0.1) # для этих лямбд будем делать RR
m.rr <- lm.ridge(dist ~ speed + I(speed^2) + I(speed^3), lambda = lambdas, data = cars)Посмотрим, как зависят оценки коэффициентов, от штрафа за высокие значения \(\hat{\beta}\):
head(coef(m.rr))## speed I(speed^2) I(speed^3)
## 0.1 -3.01537827 2.527598 -0.030849757 0.003107425
## 0.2 -1.22232251 2.075884 0.002367543 0.002367979
## 0.3 -0.49028466 1.899616 0.015025886 0.002089223
## 0.4 -0.06934523 1.803796 0.021696131 0.001944379
## 0.5 0.21734664 1.742467 0.025810896 0.001856488
## 0.6 0.43327898 1.699150 0.028599871 0.001797998plot(m.rr)
График не оцень информативен, т.к. масштаб коэффициентов существенно отличается.
Не вдаваясь в подробности скажем, что есть много алгоритмов, определяющих, какой \(\lambda\) выбрать. Помимо кросс-валидации популярны еще:
m.rr$kHKB # оценка Hoerl, Kennard and Baldwin## [1] 0.03674321m.rr$kLW # оценка Lawless and Wang ## [1] 0.527701Достоинства RR: * Оценка RR всегда есть. Даже в случае жесткой мультиколлинеарности. Даже если регрессоров больше чем, наблюдений. * Теорема: существует \(\lambda\), при котором \(MSE(\hat{\beta}_{RR})\) строго меньше, чем \(MSE(\hat{\beta}_{OLS})\). Напомню, что MSE — mean squared error, это \(MSE(\hat{\beta})=E((\hat{\beta}-\beta)^2)\). Оценки \(\hat{\lambda}_{HKB}\) и \(\hat{\lambda}_{LW}\) пытаются поймать это неизвестное \(\lambda\) минимизирующее \(MSE\).
Недостатки RR: * Нужно выбирать \(\lambda\). Конечно, есть куча алгоритмов, как это сделать. * Насколько мне известно, удачных методов постоения доверительных интервалов для \(\beta\) нет.
LASSO
LASSO очень похож на RR, отличие состоит в том, что в LASSO минимизируется функция \[ \sum_{i=1}^n (y_i-\hat{y}_i)^2+\lambda \sum_{j=1}^k |\hat{\beta}_j| \]
Заметим сразу, что функция “мерзкая”: в ней есть модуль, поэтому она не дифференциируема. Поэтому явной готовой формулы для \(\hat{\beta}_{LASSO}\) ожидать не приходится.
По другому LASSO можно сформулировать так: минимизируем обычный RSS при ограничении \(\sum_{j=1}^k |\hat{\beta}_j| \leq c\). Выписав функцию Лагранжа для этой задачи минимизации получим формулировку задачи LASSO с лямбдой.
Достоинства LASSO: * В угловых решениях некоторые \(\hat{\beta}_j\) в точности оказываются равными нулю. Т.е. LASSO без всякой проверки гипотез в каком-то смысле делит коэффициенты на “существенные” и “несущественные” * Оценка LASSO всегда есть. Даже в случае жесткой мультиколлинеарности. Даже если регрессоров больше чем, наблюдений.
Недостатки LASSO: * Нужно выбирать \(\lambda\). Конечно, есть куча алгоритмов, как это сделать. * Насколько мне известно, удачных методов постоения доверительных интервалов для \(\beta\) нет.
LASSO в R. Сначала надо ручками сваять матрицу регрессоров без единичного столбца.
X <- model.matrix(~ 0 + speed + I(speed^2) + I(speed^3), data = cars)
y <- cars$dist
fit <- glmnet(X, y, alpha = 0)
plot(fit)
Из результатов оценивания LASSO можно вытащить:
head(fit$lambda) # вектор лямбд## [1] 20817.22 18967.88 17282.83 15747.47 14348.50 13073.82head(fit$a0) # вектор оценок свободных коэффициентов## s0 s1 s2 s3 s4 s5
## 42.98000 42.82094 42.80550 42.78857 42.77000 42.74964И оценки остальных коэффициентов для каждого лямбда
fit$beta[, 1:5] # покажем только первые 5## 3 x 5 sparse Matrix of class "dgCMatrix"
## s0 s1 s2 s3
## speed 3.972130e-36 5.268060e-03 5.779479e-03 6.340312e-03
## I(speed^2) 1.302714e-37 1.727694e-04 1.895413e-04 2.079337e-04
## I(speed^3) 4.956383e-39 6.573278e-06 7.211391e-06 7.911156e-06
## s4
## speed 6.955284e-03
## I(speed^2) 2.281014e-04
## I(speed^3) 8.678467e-06Выбираем оптимальное лямбда с помощью кросс-валидации (без подробностей)
cv.fit <- cv.glmnet(X, y)
cv.fit$lambda.min## [1] 0.1503172