Загружаем нужные пакеты:
library(knitr)
opts_chunk$set(cache = TRUE)
library(ggplot2) # графики
library(glmnet) # LASSO
library(car) # compute vif()
library(MASS) # ridge regression
library(quantreg) # quantile regression
# install.packages('glmnet') # может быть нужно установить пакеты...
Возьмем встроенные в R данные по скорости автомобилей и длине тормозного пути. Это данные 1920-х годов. Скорость в милях в час, миля — это примерно \( 1.61 \) километра. Длина тормозного пути - в футах. Один фут - это примерно 30 сантиметров.
h <- cars
plot(h)
Построим полиномиальную регрессию тормозного пути от скорости
m.poly3 <- lm(dist ~ speed + I(speed^2) + I(speed^3), data = h)
Предвосхищаем вопрос: Что значит 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~(x+z+w)^3 — регрессия \( y \) на \( x \), \( z \), \( w \), \( xz \), \( xw \), \( zw \), \( xzw \). Смотрим на отчёт по регрессии
summary(m.poly3)
##
## Call:
## lm(formula = dist ~ speed + I(speed^2) + I(speed^3), data = h)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.67 -9.60 -2.23 7.08 44.69
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -19.5050 28.4053 -0.69 0.50
## speed 6.8011 6.8011 1.00 0.32
## I(speed^2) -0.3497 0.4999 -0.70 0.49
## I(speed^3) 0.0103 0.0113 0.91 0.37
##
## Residual standard error: 15.2 on 46 degrees of freedom
## Multiple R-squared: 0.673, Adjusted R-squared: 0.652
## F-statistic: 31.6 on 3 and 46 DF, p-value: 3.07e-11
Мы видим, что:
Это один из признаков мультиколлинеарности.
Мультиколлинеарность — ситуация, когда все предпосылки теоремы Гаусса-Маркова выполнены, но дисперсия оценок является слишком высокой по твоему мнению
Другие признаки мультиколлинеарности:
Корреляционная матрица регрессоров:
X <- model.matrix(~0 + speed + I(speed^2) + I(speed^3), data = h)
cor(X)
## speed I(speed^2) I(speed^3)
## speed 1.0000 0.9795 0.9389
## I(speed^2) 0.9795 1.0000 0.9884
## I(speed^3) 0.9389 0.9884 1.0000
Расчитаем индекс обусловленности, condition index, \[ CI=\sqrt{\frac{\lambda_{max}}{\lambda_{min}}} \] где \( \lambda \) — собственное число матрицы \( X'X \).
Вроде бы в некоторых источниках считают собственные числа у корреляционной или ковариационной матрицы регрессоров (???) (Лена Вакуленко) (Некоторые источники не указаны несколько дней. Вы можете улучшить эту статью :)
X <- model.matrix(~speed + I(speed^2) + I(speed^3), data = h)
XX <- t(X) %*% X
eigen <- eigen(XX)
eigen$values
## [1] 2.081e+09 1.586e+05 7.315e+01 2.719e-01
CI <- sqrt(max(eigen$values)/min(eigen$values))
CI
## [1] 87480
Наверняка есть пакет, в котором 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.1 1408.1 483.1
NB: нужно уметь рассчитывать vif ручками, зная коэффициенты R2 во вспомогательных регрессиях. Например:
m.vif1 <- lm(speed ~ I(speed^2) + I(speed^3), data = h)
r2.1 <- summary(m.vif1)$r.squared
vif1 <- 1/(1 - r2.1)
vif1
## [1] 274.1
Из модели можно вытащить ковариационную матрицу оценок коэффициентов
vcov(m.poly3)
## (Intercept) speed I(speed^2) I(speed^3)
## (Intercept) 806.8612 -186.19855 12.875757 -0.2736118
## speed -186.1986 46.25543 -3.346994 0.0733089
## I(speed^2) 12.8758 -3.34699 0.249883 -0.0055982
## I(speed^3) -0.2736 0.07331 -0.005598 0.0001276
Можно построить прогноз длины тормозного пути для машины с скоростью \( speed=10 \) миль в час.
new <- data.frame(speed = 10)
predict(m.poly3, newdata = new, interval = "confidence")
## fit lwr upr
## 1 23.79 15.89 31.69
Если нужно вытащить стандартную ошибку для \( \hat{y} \), то можно внутрь команды predict добавить опцию se.fit=TRUE.
На всякий случай, есть два типа доверительных интервалов при прогнозировании. Доверительный интервал для мат. ожидания будущего игрека — опция 'confidence' в R. И прогнозный интервал для индивидуального значения игрека — опция 'prediction' в R.
Попробуем перейти к центированным переменным.
h$spc <- h$speed - mean(h$speed)
m.poly3c <- lm(dist ~ spc + I(spc^2) + I(spc^3), data = h)
summary(m.poly3c)
##
## Call:
## lm(formula = dist ~ spc + I(spc^2) + I(spc^3), data = h)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.67 -9.60 -2.23 7.08 44.69
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 39.7503 2.8628 13.89 < 2e-16 ***
## spc 3.3258 0.8423 3.95 0.00027 ***
## I(spc^2) 0.1240 0.0712 1.74 0.08830 .
## I(spc^3) 0.0103 0.0113 0.91 0.36892
## ---
## 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.673, Adjusted R-squared: 0.652
## F-statistic: 31.6 on 3 and 46 DF, p-value: 3.07e-11
Вау! Теперь появились значимые коэффициенты!!!!!
Но это всего лишь иллюзия. У этих оценок коэффициентов другая ковариационная матрица.
vcov(m.poly3c)
## (Intercept) spc I(spc^2) I(spc^3)
## (Intercept) 8.195663 0.325739 -0.1340110 -0.0061108
## spc 0.325739 0.709415 -0.0168366 -0.0082956
## I(spc^2) -0.134011 -0.016837 0.0050694 0.0002992
## I(spc^3) -0.006111 -0.008296 0.0002992 0.0001276
И если построить доверительный интервал, то ничего не изменится.
new$spc <- 10 - mean(h$speed)
predict(m.poly3c, newdata = new, interval = "confidence")
## fit lwr upr
## 1 23.79 15.89 31.69
Модель с центрированными переменными — это всего лишь по-другому записанная исходная модель.
Кстати, нужно уметь строить доверительный интервал для прогноза “ручками” имея на руках оценки коэффициентов и оценку ковариационной матрицы.
“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 \).
Поэтому можно:
В регрессии \( y \) на \( x \) и \( z \) коэффициент при \( x \) показывает насколько изменится \( y \) при изменении \( x \) на единичку и фиксированном \( z \). В регрессии \( y \) на \( x \) коэффициент при \( x \) показывает насколько изменится \( y \) при изменении \( x \) на единичку и нефиксированном \( z \).
m.line <- lm(dist ~ speed, data = h)
summary(m.line)
##
## Call:
## lm(formula = dist ~ speed, data = h)
##
## Residuals:
## Min 1Q Median 3Q Max
## -29.07 -9.53 -2.27 9.21 43.20
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -17.579 6.758 -2.60 0.012 *
## speed 3.932 0.416 9.46 1.5e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.4 on 48 degrees of freedom
## Multiple R-squared: 0.651, Adjusted R-squared: 0.644
## F-statistic: 89.6 on 1 and 48 DF, p-value: 1.49e-12
В частности можно попытаться оставить только несколько переменных, которые “впитывают” в себя большую часть изменчивости всех исходных переменных. Этот метод называется метод главных компонент. Подробнее отдельно.
Теоретически 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 = h)
Посмотрим, как зависят оценки коэффициентов, от штрафа за высокие значения \( \hat{\beta} \):
head(coef(m.rr))
## speed I(speed^2) I(speed^3)
## 0.1 -3.01538 2.528 -0.030850 0.003107
## 0.2 -1.22232 2.076 0.002368 0.002368
## 0.3 -0.49028 1.900 0.015026 0.002089
## 0.4 -0.06935 1.804 0.021696 0.001944
## 0.5 0.21735 1.742 0.025811 0.001856
## 0.6 0.43328 1.699 0.028600 0.001798
plot(m.rr)
График не оцень информативен, т.к. масштаб коэффициентов существенно отличается.
Не вдаваясь в подробности скажем, что есть много алгоритмов, определяющих, какой \( \lambda \) выбрать. Помимо кросс-валидации популярны еще:
m.rr$kHKB # оценка Hoerl, Kennard and Baldwin
## [1] 0.03674
m.rr$kLW # оценка Lawless and Wang
## [1] 0.5277
Достоинства RR:
Недостатки RR:
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:
Недостатки LASSO:
LASSO в R. Сначала надо ручками сваять матрицу регрессоров без единичного столбца.
X <- model.matrix(~0 + speed + I(speed^2) + I(speed^3), data = cars)
y <- h$dist
fit <- glmnet(X, y)
plot(fit)
Из результатов оценивания LASSO можно вытащить:
head(fit$lambda) # вектор лямбд
## [1] 20.82 18.97 17.28 15.75 14.35 13.07
head(fit$a0) # вектор оценок свободных коэффициентов
## s0 s1 s2 s3 s4 s5
## 42.98 39.95 37.19 34.67 32.38 30.29
И оценки остальных коэффициентов для каждого лямбда
fit$beta[, 1:5] # покажем только первые 5
## 3 x 5 sparse Matrix of class "dgCMatrix"
## s0 s1 s2 s3 s4
## speed . . . . .
## I(speed^2) . 0.01146 0.0219 0.03141 0.04008
## I(speed^3) . . . . .
Выбираем оптимальное лямбда с помощью кросс-валидации (без подробностей)
cv.fit <- cv.glmnet(X, y)
cv.fit$lambda.min
## [1] 0.137
Пытается получить состоятельные оценки коэффициентов у уравнении:
\( q_{\alpha}(y)=\beta_1+\beta_2 x +\beta_3 z \)
Здесь \( q_{\alpha} \) - это квантиль порядка \( \alpha \)
m.poly3q <- rq(dist ~ speed + I(speed^2) + I(speed^3), data = h, tau = c(0.1,
0.5, 0.9))
plot(m.poly3q)
В квантильной регрессии нет простых явных формул для оценок дисперсий коэффициентов. Один из популярных способов получить их — бутстрэп. Без подробностей.
summary(m.poly3q, se = "boot")
##
## Call: rq(formula = dist ~ speed + I(speed^2) + I(speed^3), tau = c(0.1,
## 0.5, 0.9), data = h)
##
## tau: [1] 0.1
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) -31.74688 50.15200 -0.63301 0.52986
## speed 7.41086 11.52656 0.64294 0.52346
## I(speed^2) -0.40313 0.83230 -0.48435 0.63043
## I(speed^3) 0.01057 0.01879 0.56234 0.57662
##
## Call: rq(formula = dist ~ speed + I(speed^2) + I(speed^3), tau = c(0.1,
## 0.5, 0.9), data = h)
##
## tau: [1] 0.5
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) -17.22751 31.52784 -0.54642 0.58742
## speed 5.70370 7.60831 0.74967 0.45727
## I(speed^2) -0.25463 0.57576 -0.44225 0.66038
## I(speed^3) 0.00761 0.01357 0.56031 0.57798
##
## Call: rq(formula = dist ~ speed + I(speed^2) + I(speed^3), tau = c(0.1,
## 0.5, 0.9), data = h)
##
## tau: [1] 0.9
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) -6.13904 66.79214 -0.09191 0.92717
## speed 4.16336 17.20828 0.24194 0.80990
## I(speed^2) -0.04768 1.27408 -0.03743 0.97031
## I(speed^3) 0.00388 0.02860 0.13578 0.89259
Изобразим на одном графике линию обычной регрессии и линии квантильных регрессий.
# стащить из коэнкэра клёвый график :)