Свежий пакет "conquer".
Базовым пакетом для квантильной регрессии является пакет quantreg, про который очень хорошо было написано вот тут.
Но сегодня мы поговорим о пакете, который был презентован 1 ноября.
Немного математики
Если есть выборка размера n, и классическая матрица признаки-ответы (X-Y), то задача квантильной регрессии для квантиля τ записывается как
![](https://habrastorage.org/getpro/habr/upload_files/dda/347/1e5/dda3471e5b2ffb16a9e1d92f68186e7b.png)
Оценки коэффициентов β вычисляются по формуле
![](https://habrastorage.org/getpro/habr/upload_files/3c8/582/8be/3c85828bec03a8e242c0aa67a4636791.png)
![А это - функция потерь квантильной регрессии А это - функция потерь квантильной регрессии](https://habrastorage.org/getpro/habr/upload_files/983/bf4/ca2/983bf4ca290aad7df4b39829a1d62e66.png)
Развитие методов регуляризации в обычной регрессии логичным образом перешло и в квантильную регрессию. Постановка задачи регуляризации в квантильной регрессии такая же, как и в обычной:
![Где L(f) – функция потерь, P(f) – штрафная функция, λ – коэффициент Где L(f) – функция потерь, P(f) – штрафная функция, λ – коэффициент](https://habrastorage.org/getpro/habr/upload_files/4c9/c63/941/4c9c63941e6ee0cbfc0c1bc0be97652a.png)
Более полно эту задачу можно расписать как
![где Ω = [a;b] – интервал, охватывающий все точки, f(x) – функция регрессии где Ω = [a;b] – интервал, охватывающий все точки, f(x) – функция регрессии](https://habrastorage.org/getpro/habr/upload_files/852/655/536/85265553637194c3ea1fcf7f8abbcb62.png)
Решение этого функционала сводится к задаче квадратичного программирования. В 1994 году было показано, что можно убрать квадрат из штрафной функции, заменив ее задачей
![](https://habrastorage.org/getpro/habr/upload_files/9ba/c0e/12b/9bac0e12bfc9ffbc4c99f5fc46c04e9f.png)
Разница между получающимися функциями представлена на рисунке ниже
![Цитируется по: Koenker. R., Chernozhukov, V., He, X., Peng, L. Handbook of Quantile Regression. CRC Press. 2018. ISBN 978-1-4987-2528-6, с.47 Цитируется по: Koenker. R., Chernozhukov, V., He, X., Peng, L. Handbook of Quantile Regression. CRC Press. 2018. ISBN 978-1-4987-2528-6, с.47](https://habrastorage.org/getpro/habr/upload_files/727/dcb/939/727dcb9393e804a5ffeb1df13200d803.png)
Фактически, последняя формула повторяет постановку задачи лассо-регрессии, в которой берется модуль от некоторого оценщика. Еще одно ее преимущество – вычислительная простота, поскольку конкретные значения коэффициентов можно получить в результате решения задачи линейного программирования.
Последующее развитие постановки задачи привело ее к виду
![](https://habrastorage.org/getpro/habr/upload_files/a72/448/7e3/a724487e3884c222e539eae8ad69728b.png)
где степень производной k предполагалось выбирать в зависимости от свойств данных:
- если функция f(x) предполагает наличие мгновенных скачков, то выбирается k=1
- если функция f(x) непрерывна, но имеет точки перегиба, то k=2
- если функция f(x) непрерывна и дифференцируема, то k=3
Далее усилия теоретиков были направлены на определение оптимального значения параметра λ. Как правило, для этого используется алгоритм перекрестной проверки (кросс-валидация с разбиением на несколько групп)
В сентябре 2021 года коллектив китайских ученых выпустил препринт статьи (https://arxiv.org/abs/2109.05640), в котором предложил способ борьбы с систематическими ошибками, которые дает лассо-регуляризация в квантильной регрессии
Постановка задачи нахождения коэффициентов квантильной регрессии записывается теперь так:
![](https://habrastorage.org/getpro/habr/upload_files/2a1/b12/03f/2a1b1203fc7ea2b0464c9b2654f6eafc.png)
где qλ – некоторая функция от λ
Кроме того, первое слагаемое функции потерь заменяется ими на ядро K(), и предлагается общий итеративный алгоритм нахождения коэффициентов квантильной регрессии. Именно он и используется в новом пакете conquer
Немного расчетов
Для расчетов возьмем базу данных Schooling из пакета Ecdat. В ней представлены данные о зарплате американских рабочих. Опрос был проведен в 1976 году.
В качестве зависимой переменной выберем заработную плату, независимых – количество лет образования, возраст, результаты теста iq и опыт работы.
Простую модель квантильной регрессии можно построить с помощью команды conquer(), которая содержит следующие параметры:
X, Y – матрица независимых и вектор зависимого признаков
tau – квантиль (при tau = 0.5 получается медианная регрессия)
kernel - определяет, какое ядро сглаживания K() используется (выбор из следующий значений: "Gaussian", "logistic", "uniform", "parabolic", "triangular")
ci – параметр, определяющий необходимость расчета доверительных интервалов коэффициентов регрессии
h, checkSing, tol, iteMax, alpha – параметры алгоритма расчета
library(Ecdat)
library(tidyverse)
library(conquer)
Base1<-na.omit(Ecdat::Schooling)
glimpse(Base1)
Y <- Base1$wage76
X <- as.matrix(Base1[,c(7,9,25,28)])
mod_1<-conquer(X,Y,tau=0.9)
mod_1$coeff
![](https://habrastorage.org/getpro/habr/upload_files/2f9/34e/968/2f934e9680f36fdb3f7a5efb3154dd78.png)
В итоге выводятся коэффициенты регрессии в порядке следования параметров, первое число – значение свободного члена
Также по результатам расчетов можно вывести:
ite - Число итераций до сходимости.
residual - Вектор полученных моделью значений
perCI, pivCI, normCI – Доверительные интервалы коэффициентов регрессии
Посчитаем модель с другим ядром и выведем доверительные интервалы
![](https://habrastorage.org/getpro/habr/upload_files/6c6/5ed/3f3/6c65ed3f38a20c959ab0feac299a5036.png)
Как видно, коэффициенты, во-первых, отличаются не сильно, и, во-вторых, относительно устойчивы (доверительные интервалы не включают 0)
Функция conquer.cv.reg – позволяет реализовать квантильную регрессию с применением методов регуляризации (за это отвечает параметр penalty, которое может принимать одно из трех значений: "lasso", "scad", "mcp") и кросс-валидации (параметр kfolds отвечает за число блоков). Посчитаем квантильную регрессию с лассо-регуляризацией (придется подождать, считает дольше)
![](https://habrastorage.org/getpro/habr/upload_files/d83/f68/5d6/d83f685d6caa06255e7d1dbca08a230b.png)
Как видно, регуляризация посчитала, что последнюю переменную можно не учитывать.
Также крайне удобна функция conquer.process. Она позволяет посчитать сразу целый блок квантильных регрессий для разных значений квантиля:
tau_t <- seq(0.1, 0.9, by = 0.05)
mod_3 <-conquer.process(X,Y,tauSeq = tau_t)
colnames(mod_3$coeff) <- tau_t
t(mod_3$coeff)
![](https://habrastorage.org/getpro/habr/upload_files/89e/8e9/668/89e8e9668157511e92aa8963c7fc2f62.png)
Собственно, это показывает, что для разных квантилей коэффициенты регрессии могут сильно отличаться
Кроме того, с небольшой модификацией можно использовать эти функции для отбора признаков. Посмотрим, как будут меняться коэффициенты функции при разных значениях λ
lambda_v <- 10^(seq(-4, -1, by = 0.1))
Res <- matrix(0, nrow = 31, ncol = 5)
for(i in 1:31) {
mod_4 <- conquer.reg(X,Y,lambda = lambda_v[i],tau = 0.9, penalty = c("lasso"))
Res[i,]<-mod_4$coeff}
rownames(Res)<-lambda_v
Res
![](https://habrastorage.org/getpro/habr/upload_files/ab4/a9c/62a/ab4a9c62ab9444fafcd9b59d43a187e3.png)
Эту динамику значений вполне можно использовать как характеристику важности переменных. Таким образом, ранжирование переменных по значимости будет выглядеть так: возраст → количество лет учебы → результаты теста iq → опыт работы.
kenoma
Взглянуть то взглянули, а можно хотя бы пару слов о том, чем же этот подход лучше\хуже остальных и в каких классах задач есть выгода от его применения?
acheremuhin Автор
В целом, он позволяет моделировать зависимые переменные со сложным / меняющимся распределением. Условно, при х = 10 зависимая переменная распределена нормально, при х = 20 - экспоненциально. Если вас интересуют конкретные примеры, то укажите конкретную отрасль, чтобы говорить более предметно
kenoma
В вашем примере параметр х скрытый или явный?
acheremuhin Автор
В моем условном примере - явный. А, например, в статье, там где модель mod_3, верхняя строчка (соответствующая tau = 0.1) - это модель регрессии для значения, отсекающего 10 % самых низких зарплат, а нижняя строчка (соответствующая tau = 0.9) - это модель регрессии для значения, отсекающего 10 % самых высоких зарплат