Обзор пакета 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 - число наблюдений в совокупности). Однако тут мы попадаем в классическую ловушку построения полиномиальных моделей - идеальная аппроксимация ведет к переобучению, что на графике отображается примерно вот так

Картинка взята с https://ppt-online.org/22644

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

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

Справа - проверка на выполнимость условий Гаусса-Маркова с помощью пакета gvlma

Теперь построим модель с 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).

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