Подключим пакеты:

library("quantreg") # квантильная регрессия
library("ggplot2") # графики
library("dplyr") # манипуляции с данными

Обычная регрессия предполагает, что условное среднее переменной \(y\) линейно зависит от регрессора \(x\):

\[ E(y_i|x_i) = \beta_1 + \beta_2 x_i \]

Медианная регрессия предполагает, что от регрессоров линейно зависит условная медиана зависимой переменной:

\[ Med(y_i|x_i)=\beta_1+\beta_2 x_i \]

Оценки медианной регрессии находятся путём минимизации суммы модулей. Если, как и в обычной регрессии, обозначить \(\hat y_i = \hat\beta_1 + \hat\beta_2 x_i\), то целевая функция имеет вид:

\[ \sum |y_i - \hat y_i| \to \min_{\hat\beta} \]

Обобщением медианной регрессии является квантильная регрессия. Квантильная регрессия предполагает, что условная квантиль порядка \(\tau\) переменной \(y\) линейно зависит от регрессора \(x\):

\[ q_{\tau}(y_i|x_i)=\beta_1+\beta_2 x_i \]

Здесь \(q_{\tau}\) - это квантиль порядка \(\tau\).

Задача оптимизации для квантильной регрессии имеет вид:

\[ \sum w_i|y_i - \hat y_i| \to \min_{\hat\beta}, \] где веса \(w_i\) показывают, насколько большой штраф получается в зависимости от того, перескочил ли прогноз \(\hat y_i\) фактическое значение \(y_i\)

\[ w_i =\begin{cases} \tau, \text{ если } \hat y_i < y_i \\ 1-\tau, \text{ если } \hat y_i \geq y_i \end{cases} \]

Возьмём встроенный в R набор данных по стоимости бриллиантов и перейдем к логарифму цены бриллианта и логарифму массы:

diams <- mutate(diamonds, 
                log_price = log(price),
                log_carat = log(carat))

Оценим обычную и несколько квантильных регрессий:

model_ols <- lm(data = diams, log_price ~ log_carat)
model_q <- rq(data = diams, log_price ~ log_carat,
              tau = c(0.01, 0.1, 0.5, 0.9, 0.99))

Посмотрим на описания всех пяти моделей

summary(model_q)
## 
## Call: rq(formula = log_price ~ log_carat, tau = c(0.01, 0.1, 0.5, 0.9, 
##     0.99), data = diams)
## 
## tau: [1] 0.01
## 
## Coefficients:
##             Value     Std. Error t value   Pr(>|t|) 
## (Intercept)   7.73780   0.02101  368.23594   0.00000
## log_carat     1.49805   0.02482   60.35184   0.00000
## 
## Call: rq(formula = log_price ~ log_carat, tau = c(0.01, 0.1, 0.5, 0.9, 
##     0.99), data = diams)
## 
## tau: [1] 0.1
## 
## Coefficients:
##             Value      Std. Error t value    Pr(>|t|)  
## (Intercept)    8.13418    0.00219 3706.09415    0.00000
## log_carat      1.69432    0.00326  519.66516    0.00000
## 
## Call: rq(formula = log_price ~ log_carat, tau = c(0.01, 0.1, 0.5, 0.9, 
##     0.99), data = diams)
## 
## tau: [1] 0.5
## 
## Coefficients:
##             Value      Std. Error t value    Pr(>|t|)  
## (Intercept)    8.43377    0.00156 5420.28232    0.00000
## log_carat      1.65124    0.00236  699.07418    0.00000
## 
## Call: rq(formula = log_price ~ log_carat, tau = c(0.01, 0.1, 0.5, 0.9, 
##     0.99), data = diams)
## 
## tau: [1] 0.9
## 
## Coefficients:
##             Value      Std. Error t value    Pr(>|t|)  
## (Intercept)    8.80055    0.00306 2875.63314    0.00000
## log_carat      1.71922    0.00313  548.93883    0.00000
## 
## Call: rq(formula = log_price ~ log_carat, tau = c(0.01, 0.1, 0.5, 0.9, 
##     0.99), data = diams)
## 
## tau: [1] 0.99
## 
## Coefficients:
##             Value     Std. Error t value   Pr(>|t|) 
## (Intercept)   9.15211   0.01928  474.77492   0.00000
## log_carat     1.83616   0.02236   82.11804   0.00000

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

summary(model_q, se = "boot")

Можно на графике посмотреть, как меняются коэффициенты при изменении квантиля:

plot(model_q)

На графике видно, что у самых дорогих бриллиантов цена более чувствительна к массе, чем у самых дешёвых. Для бриллиантов в среднем ценовом диапазоне эластичность цены по массе стабильна.

Изобразим на одном графике линию обычной регрессии и линии квантильных регрессий.

base <- ggplot(data = diams, aes(x = log_carat, y = log_price)) + 
  geom_point() + ggtitle("Зависимость цены от массы бриллианта") +
  xlab("Логарифм массы в каратах") + ylab("Логарифм цены в долларах")
base + stat_smooth(method = "lm", se = FALSE) +
  stat_smooth(method = "rq", se = FALSE, method.args = list(tau = 0.01), col = "red") +
  stat_smooth(method = "rq", se = FALSE, method.args = list(tau = 0.99), col = "red")

Без опции se = FALSE при методе rq вылезает ошибка.

Todo:

Почиташки: