Мы не идем простыми путями или о применении P-сплайнов в регрессии

Обзор пакета SOP

С чем и как работаем

База данных для работы — данные по округам США (за 2017 год).

Задача: сравнить точность построенных моделей регрессии с использованием P-сплайнов при прогнозировании средней величины дохода в округе (переменная Income)

Алгоритм работы:

С помощью алгоритмов отбора признаков (возьмем, к примеру, недавно придуманный алгоритм, имитирующий охоту горбатых китов; на русском — https://russianblogs.com/article/89241016039/) ранжируем признаки по порядку важности.

Потом разделяем выборку на тренировочную и тестовую (80/20), строим уравнения регрессии (обычной и со сплайнами) с 1 и 2 переменными по тренировочной выборке, оцениваем степень точности прогнозов по тестовой выборке.

Минутка теории

Р-сплайн — модификация В-сплайнов, впервые представленная в 1996 году в статье Поля Эйлера и Брайана Маркса (https://projecteuclid.org/journals/statistical-science/volume-11/issue-2/Flexible-smoothing-with-B-splines-and-penalties/10.1214/ss/1038425655.full). В-сплайн представляет собой конструкцию вида

где n — степень полинома

Соответственно, сумма таких сплайнов позволяет подогнать почти любую функцию при достаточном количестве узлов разбиения t. Задача аппроксимации В-сплайна, таким образом, сводится к минимизации функционала

относительно вектора параметров а (здесь уже n — количество полиномов, m — число наблюдений в совокупности). Однако тут мы попадаем в классическую ловушку построения полиномиальных моделей — идеальная аппроксимация ведет к переобучению, что на графике отображается примерно вот так

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

Оригинальная статья Эйлера и Маркса, во-первых, содержит рекомендации по упрощению процедуры оптимизации этого выражения (путем преобразования второго слагаемого), во-вторых (и самых главных) по-другому ставит задачу, переводя ее в вероятностную плоскость — давайте подберем такие значения а, для которых вероятность появления штрафной функции будет минимальна (это и есть Р-сплайны — сплайны со штрафом). Остальные 28 страниц этой статьи посвящены решению этой задачи и комментариям рецензентов, нам важно зафиксировать постановку задачи и тот факт, что задача решалась для одномерного случая, и то, что параметр сглаживания предлагалось подбирать вручную с помощью AIC-критерия.

В 2014 (https://link.springer.com/article/10.1007/s11222-014-9464-2) и 2018 году (https://arxiv.org/abs/1801.07278) мультинациональная группа исследователей под руководством Поля Эйлера представила две статьи, в которых изложена модификация метода аппроксимации Р-сплайнами, реализованная для случая нескольких переменных и с автоматическим вычислением оптимального коэффициента сглаживания. В 2021 году этот алгоритм был реализован на языке программирования R в пакете SOP.

Расчеты в R

Первая часть — получаем список переменных, отсортированных по убыванию важности (берем 20 подвыборок по 80 % от основной выборки, проводим отбор признаков по подвыборке и потом считаем, сколько раз переменные оказались выбраны в модель)

library(tidyverse) library(ModelMetrics) library(openxlsx) library(readr) library(gvlma) Base <- read_csv(".../acs2017_census_tract_data.csv") Base_add <- as.data.frame(na.omit(Base[,c(5:37)]))  glimpse(Base_add) # 11 - 14 переменные: зависимые tr.index <- sample(1:74001, 74001*0.8)  # Отбираем переменные mod <- 20 library("FSinR") evaluator <- filterEvaluator('determinationCoefficient') # Чем мы мерим woa_search <- whaleOptimization() set.seed(100) Res<-as.data.frame(matrix(0, nrow = mod, ncol = 7)) colnames(Res) <- c("R2", "MAE", "MSE", "Approximation", "Appr_5_perc","Appr_15_perc", "Time") Res_per_1<-as.data.frame(matrix(0, nrow = mod, ncol = 29)) for(k in 1:mod) { print(k)   tr.index <- sample(1:74001, 74001*0.8)    Base_tr <- Base_add[tr.index, c(1:10,14:33)]   Base_test <- Base_add[-tr.index, c(1:10,14:33)]   model_1<-woa_search(Base_tr, 'Income', evaluator)   Res_per_1[k,] <- model_1$bestFeatures } colnames(model_1$bestFeatures)[order(apply(Res_per_1,2,sum), decreasing = TRUE)]

Таким образом, две переменные, наиболее сильно влияющие на средний доход в округе — доля бедных (как неожиданно) и доля работающих.

Построим сначала простую линейную модель для переменной Poverty

Теперь построим модель с P-сплайном. Важно — синтаксис построения модели отличается, в записи формулы должна быть обязательно переменная вида f(name_var_1, …) или ad(name_var), где f(name_var_1, …) обозначает построение одно- или многомерного сплайна без адаптивного сглаживания, а ad(name_var) — построение сплайна от конкретной переменной с адаптивным сглаживанием.

Построим модель без адаптивного сглаживания

m1_f <- sop(Income ~ f(Poverty), data = Base_tr_1) summary(m1_f)

Коэффициенты, кстати, получились совершенно иными.Построим теперь модель с адаптивным сглаживанием.

m1_ad <- sop(Income ~ ad(Poverty), data = Base_tr_1) m1_ad <- sop(Income ~ ad(Poverty, pord = 4, degree = 5), data = Base_tr_1) summary(m1_ad)

В некоторых случаях (как в первой строке выше) при заданных параметрах возникают проблемы с расчетами. Тогда можно повысить степень полинома у сплайна (параметр degree) и повысить степень производной в штрафующем коэффициенте (параметр pord). Тогда автоматически и усложнится само уравнение регрессии — в нашем случае оно получилось кубическим.

Результаты модели с 2 переменными, коэффициенты которой посчитаны с помощью МНК, выглядят так:

Модель без адаптивного сглаживания с 2 переменными строится так

m2_f <- sop(Income ~ f(Poverty,Professional), data = Base_tr_1) summary(m2_f)

Вполне можно одну переменную использовать и саму по себе, без сплайнов, а вторую — со сплайнами

m2_f_1 <- sop(Income ~ Poverty + f(Professional), data = Base_tr_1) summary(m2_f_1)

Также можно построить модель с адаптивным сглаживанием разных порядков для разных переменных

m2_ad <- sop(Income ~ ad(Poverty, pord = 4, degree = 5)+ad(Professional, pord = 3, degree = 4), data = Base_tr_1) summary(m2_ad)

Сводим результаты всех моделей в одну табличку

Видим, что модель с 2 переменными без адаптивного сглаживания обладает наименьшей средней ошибкой, более 15 % наблюдений на тестовой выборке предсказаны с погрешностью менее 5 %, а почти 42 % наблюдений — с погрешностью менее 15 %. Соответственно, эта модель значительно качественнее, чем обычные линейные модели, а сам метод P-сплайнов может эффективно применяться в случае невыполнения условий Гаусса-Маркова и смещенности оценок МНК, а также может использоваться и для подбора видов полиномиальных зависимостей (все материалы можно найти на https://github.com/acheremuhin/P-spline).