Инфраструктура для обучения моделей на R: rsample и recipes

При создании моделей машинного обучения довольно часто концентрируются на настройке гиперпараметров, в то время как предварительная обработка данных и создание новых признаков (feature engineering) может быть гораздо более продуктивных вложением сил и времени. В R настраивать гиперпараметры можно множеством способов, популярным вариантом по историческим причинам является использование пакета caret. Что касается этапа преобразования и создания признаков, то частой ошибкой является его выполнение на всем наборе данных, а это чревато переобучением. Такие операции, как преобразование Бокса-Кокса, биннинг количественных переменных, mean/median target encoding, импутация средним (или более продвинутые методы заполнения пропусков) и даже обычная нормализация являются обучаемыми: их параметры должны оцениваться на обучающей части выборки и применяться и на ней, и на валидационной, и на тестовой выборке. Также целесообразно пробовать разные методы преобразования, импутации и т.д. в составе “конвейера” обучения моделей (модель - это сам алгоритм + предобработка данных), то есть добавлять для этих этапов гиперпараметры, проверяемые в ходе grid search-а.
В Python все описанное реализовано к классе sklearn.model_selection.GridSearchCV: в него можно передать сетку гиперпараметров вида

param_grid = {
    # гиперпараметры импутации и биннинга
    'tr__num__imp__strategy': ['mean', 'median', 'constant'],
    'tr__num_bin__bi__n_bins': [2, 3],
    'tr__cat__imp__strategy': ['most_frequent', 'constant'],
    # гиперпараметры модели, в данном случае логрег с регуляризацией
    'lr__C': [.001, .01, .1]
}

… и конвейер операций, выполняемых внутри ресемплов по описанной выше схеме (с обучением на обучающей выборке и применением на обучающей и валидационной):

# создаем конвейер для количественных переменных,
# которые не будут подвергнуты биннингу
num_pipe = Pipeline([
    ('imp', SimpleImputer()),
    ('yeo_john', PowerTransformer(method='yeo-johnson', standardize=True))
])

# создаем конвейер для количественных переменных,
# которые будут подвергнуты биннингу
num_bin_pipe = Pipeline([
    ('imp', SimpleImputer()),
    ('power', FunctionTransformer(power, validate=False)),
    ('bi', KBinsDiscretizer(encode='onehot-dense'))
])

# создаем конвейер для категориальных переменных
cat_pipe = Pipeline([
    ('imp', SimpleImputer()),
    ('ohe', OneHotEncoder(sparse=False, handle_unknown='ignore'))
])


# создаем список трехэлементных кортежей, в котором
# первый элемент кортежа - название конвейера с
# преобразованиями для определенного типа признаков
transformers = [('num', num_pipe, num_columns),
                ('num_bin', num_bin_pipe, num_bin_columns),
                ('cat', cat_pipe, cat_columns)]


# передаем список трансформеров в ColumnTransformer
transformer = ColumnTransformer(transformers=transformers)


# задаем итоговый конвейер
ml_pipe = Pipeline([('tr', transformer), 
                    ('lr', LogisticRegression(solver='lbfgs', max_iter=200))])

В R подобный функционал можно реализовать при помощи пакетов rsample и recipes. В составе документации есть достаточно комплексные примеры, например, Recipes with rsample и Grid Search Tuning of Keras Models. Но решения задачи подбора (суб)оптимальных вариантов предварительной обработки в ходе единого процесса обучения моделей показано не было; этот недостаток мы и будем восполнять.
Воспользуемся набором данных Vodafone_missing.csv (source):

library(rsample)

## Loading required package: tidyr

library(recipes)

## Loading required package: dplyr

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

## 
## Attaching package: 'recipes'

## The following object is masked from 'package:stats':
## 
##     step

library(purrr)
library(yardstick)

## For binary classification, the first factor level is assumed to be the event.
## Set the global option `yardstick.event_first` to `FALSE` to change this.

library(glmnet)

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following object is masked from 'package:tidyr':
## 
##     expand

## Loading required package: foreach

## 
## Attaching package: 'foreach'

## The following objects are masked from 'package:purrr':
## 
##     accumulate, when

## Loaded glmnet 2.0-18

data <- read.csv2("Vodafone_missing.csv", sep = ";", 
                  stringsAsFactors = FALSE,
                  na.strings = c("NA", ""))
data$churn <- factor(data$churn, levels = c(0, 1))

str(data)

## 'data.frame':    1000 obs. of  12 variables:
##  $ region : chr  "Region 2" "Region 3" "Region 3" "Region 2" ...
##  $ tenure : int  13 11 68 NA 23 41 45 38 45 68 ...
##  $ age    : int  44 33 52 33 30 NA 22 35 NA 41 ...
##  $ marital: chr  "mar" "mar" "mar" NA ...
##  $ address: int  9 7 24 12 9 17 NA 5 7 21 ...
##  $ income : int  64 136 116 NA 30 78 19 76 NA NA ...
##  $ employ : int  5 5 29 NA 2 16 NA 10 31 22 ...
##  $ retire : chr  "no" "no" "no" "no" ...
##  $ gender : chr  "f" "f" NA NA ...
##  $ reside : int  2 6 NA 1 4 1 5 3 5 NA ...
##  $ custcat: chr  "cat 1" "cat 4" NA "cat 1" ...
##  $ churn  : Factor w/ 2 levels "0","1": 2 2 1 2 1 1 2 1 1 1 ...

head(data)

##     region tenure age marital address income employ retire gender reside
## 1 Region 2     13  44     mar       9     64      5     no      f      2
## 2 Region 3     11  33     mar       7    136      5     no      f      6
## 3 Region 3     68  52     mar      24    116     29     no   <NA>     NA
## 4 Region 2     NA  33    <NA>      12     NA     NA     no   <NA>      1
## 5 Region 2     23  30     mar       9     30      2     no      f      4
## 6 Region 2     41  NA   unmar      17     78     16     no      m      1
##   custcat churn
## 1   cat 1     1
## 2   cat 4     1
## 3    <NA>     0
## 4   cat 1     1
## 5   cat 3     0
## 6   cat 3     0

Как видим, у нас есть как количественные, так и категориальные переменные, а в данных довольно много пропусков.
Разобьем весь набор данных на обучающий и тестовый в пропорции 70/30:

set.seed(42)
data_split <- initial_split(data, prop = 0.3, strata = "churn")
train_data <- training(data_split)
test_data <- testing(data_split)

Отберем переменные, подлежащие разной предварительной обработке:

cat_columns <- names(data)[sapply(data, is.character)]
num_columns <- setdiff(names(data), c(cat_columns, "churn"))
num_bin_columns <- c("tenure", "age")

Создадим схему перекрестной проверки со стратифицированной по целевой переменной с разбивкой на 5 блоков:

set.seed(42)
cv_splits <- vfold_cv(train_data, v = 5, strata = "churn")

Строить будем модель логистической регрессии (пакет glmnet), проверяя три значения гиперпараметра lambda (сила регуляризции) и три значения гиперпараметра alpha (пропорция l1- и l2-штрафов). Также проверим два варианта импутации пропусков и два варианта биннинга (разбивка на 2 и на 3 категории):

param_grid <- expand.grid(imp_strategy = c("mean_imp", "median_imp"),
                          n_bins = c(2, 3), 
                          alpha = c(0, 0.5, 1),
                          lambda = c(0.1, 0.01, 0.001)
                          )

Нам понадобится функция, создающая “рецепт”, которая будет параметризована по типу импутации и по вариантам биннинга:

create_recipe <- function(data = train_data,
                          imp_strategy = "mean_imp", 
                          n_bins = 3) {
 
  # Импутация константным значением не реализована
  imp_methods <- list(mean_imp = step_meanimpute, 
                      median_imp = step_medianimpute)

  rec_obj <- recipe(churn ~ ., data = data) %>%
    # Импутация средним или медианой для всех количественных
    imp_methods[[imp_strategy]](one_of(num_columns)) %>%
    # Преобразование Йео-Джонсона для всех количественных
    step_YeoJohnson(one_of(num_columns)) %>%
    # Стандартизация
    step_center(one_of(num_columns)) %>%
    step_scale(one_of(num_columns)) %>%
    
    # Возводим в квадрат num_bin_columns (импутация для них уже выполнена)
    # Отдельной функции для этого нет (нужно писать кастомный step), 
    # поэтому реализуем вручную через mutate
    # Такие трансформации обычно делают вне ресемплов
    step_mutate(tenure = tenure^2, age = age^2) %>%
    # Выполняем биннинг
    step_discretize(one_of(num_bin_columns), 
                    options = list(cuts = n_bins, keep_na = FALSE)) %>%
    # Преобразовываем метки бинов в one-hot
    step_dummy(one_of(num_bin_columns), one_hot = TRUE) %>%
    
    # Импутация самым частым значением для всех категориальных
    step_modeimpute(one_of(cat_columns)) %>%
    # one-hot для всех категориальных
    step_dummy(one_of(cat_columns), one_hot = TRUE)
  
  return(rec_obj)
}

Функция для обучения модели, возвращающая значение выбранной метрики (в данном случае accuracy):

logreg_fit_pred <- function(split, 
                            data = train_data,
                            imp_strategy = "mean_imp", 
                            n_bins = 3, 
                            lambda = 0.1, 
                            alpha = 0) {
  
  rec_obj <- create_recipe(data = data, 
                           imp_strategy = imp_strategy,
                           n_bins = n_bins)
  rec_obj <- prepper(split_obj = split, rec_obj, retain = FALSE)
  
  train <- bake(rec_obj, 
                 new_data = analysis(split),
                 everything())
  val <- bake(rec_obj, 
              new_data = assessment(split),
              everything())
  
  train_x <- 
    train %>%
    select(-churn) %>%
    as.matrix()
  train_y <- train$churn
  
  val_x <- val %>%
    select(-churn) %>%
    as.matrix()
  val_y <- val$churn
  
  model <- glmnet(x = train_x, 
                  y = train_y, 
                  family = "binomial",
                  lambda = lambda,
                  alpha = alpha)
  
  preds <- predict(model, 
                    type = "response",
                    newx = val_x)
  preds <- ifelse(preds > 0.5, 1, 0)
  preds <- factor(preds, levels = c(0, 1))
  result <- accuracy(data.frame(preds = preds, val_y = val_y), 
                     preds, 
                     val_y)
  
  return(result)
}

logreg_fit_pred(cv_splits$splits[[1]], data = train_data)

## # A tibble: 1 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.721

Функция, обучающая модели для всех комбинаций значений гиперпараметров:

across_grid <- function(split, grid) {
  metrics <- lapply(seq_len(nrow(grid)),
                    function(i) logreg_fit_pred(split, 
                                                data = train_data,
                                                grid[i, "imp_strategy"],
                                                grid[i, "n_bins"],
                                                grid[i, "lambda"],
                                                grid[i, "alpha"]))
  result <- do.call(rbind, metrics)
  return(cbind(grid, result))
}

across_grid(cv_splits$splits[[1]], param_grid)

##    imp_strategy n_bins alpha lambda  .metric .estimator .estimate
## 1      mean_imp      2   0.0  0.100 accuracy     binary 0.6885246
## 2    median_imp      2   0.0  0.100 accuracy     binary 0.6885246
## 3      mean_imp      3   0.0  0.100 accuracy     binary 0.7213115
## 4    median_imp      3   0.0  0.100 accuracy     binary 0.7213115
## 5      mean_imp      2   0.5  0.100 accuracy     binary 0.7213115
## 6    median_imp      2   0.5  0.100 accuracy     binary 0.7213115
## 7      mean_imp      3   0.5  0.100 accuracy     binary 0.7213115
## 8    median_imp      3   0.5  0.100 accuracy     binary 0.7213115
## 9      mean_imp      2   1.0  0.100 accuracy     binary 0.7213115
## 10   median_imp      2   1.0  0.100 accuracy     binary 0.7213115
## 11     mean_imp      3   1.0  0.100 accuracy     binary 0.7213115
## 12   median_imp      3   1.0  0.100 accuracy     binary 0.7213115
## 13     mean_imp      2   0.0  0.010 accuracy     binary 0.6721311
## 14   median_imp      2   0.0  0.010 accuracy     binary 0.6721311
## 15     mean_imp      3   0.0  0.010 accuracy     binary 0.7213115
## 16   median_imp      3   0.0  0.010 accuracy     binary 0.7049180
## 17     mean_imp      2   0.5  0.010 accuracy     binary 0.6721311
## 18   median_imp      2   0.5  0.010 accuracy     binary 0.6721311
## 19     mean_imp      3   0.5  0.010 accuracy     binary 0.7377049
## 20   median_imp      3   0.5  0.010 accuracy     binary 0.7377049
## 21     mean_imp      2   1.0  0.010 accuracy     binary 0.6885246
## 22   median_imp      2   1.0  0.010 accuracy     binary 0.6885246
## 23     mean_imp      3   1.0  0.010 accuracy     binary 0.7540984
## 24   median_imp      3   1.0  0.010 accuracy     binary 0.7540984
## 25     mean_imp      2   0.0  0.001 accuracy     binary 0.6721311
## 26   median_imp      2   0.0  0.001 accuracy     binary 0.6721311
## 27     mean_imp      3   0.0  0.001 accuracy     binary 0.7213115
## 28   median_imp      3   0.0  0.001 accuracy     binary 0.7213115
## 29     mean_imp      2   0.5  0.001 accuracy     binary 0.6721311
## 30   median_imp      2   0.5  0.001 accuracy     binary 0.6721311
## 31     mean_imp      3   0.5  0.001 accuracy     binary 0.7213115
## 32   median_imp      3   0.5  0.001 accuracy     binary 0.7213115
## 33     mean_imp      2   1.0  0.001 accuracy     binary 0.6721311
## 34   median_imp      2   1.0  0.001 accuracy     binary 0.6721311
## 35     mean_imp      3   1.0  0.001 accuracy     binary 0.7213115
## 36   median_imp      3   1.0  0.001 accuracy     binary 0.7213115

Остается лишь выполнить собственно перекрестную проверку:

result <- lapply(cv_splits$splits, across_grid, grid = param_grid)
library(data.table)

## 
## Attaching package: 'data.table'

## The following object is masked from 'package:purrr':
## 
##     transpose

## The following objects are masked from 'package:dplyr':
## 
##     between, first, last

result <- data.table::rbindlist(result, idcol = "split")
result[, .(acc = mean(.estimate)), 
       by = .(imp_strategy, n_bins, lambda, alpha)
       ][
         order(acc, decreasing = TRUE),]

##     imp_strategy n_bins lambda alpha       acc
##  1:     mean_imp      3  0.100   0.0 0.7344818
##  2:   median_imp      3  0.100   0.0 0.7344818
##  3:     mean_imp      3  0.001   0.0 0.7245346
##  4:   median_imp      3  0.001   0.0 0.7245346
##  5:     mean_imp      3  0.001   0.5 0.7245346
##  6:   median_imp      3  0.001   0.5 0.7245346
##  7:     mean_imp      2  0.100   0.5 0.7243123
##  8:   median_imp      2  0.100   0.5 0.7243123
##  9:     mean_imp      3  0.100   0.5 0.7243123
## 10:   median_imp      3  0.100   0.5 0.7243123
## 11:     mean_imp      2  0.100   1.0 0.7243123
## 12:   median_imp      2  0.100   1.0 0.7243123
## 13:     mean_imp      3  0.100   1.0 0.7243123
## 14:   median_imp      3  0.100   1.0 0.7243123
## 15:   median_imp      2  0.010   0.0 0.7243123
## 16:     mean_imp      2  0.010   0.0 0.7242012
## 17:   median_imp      2  0.001   0.0 0.7242012
## 18:     mean_imp      3  0.010   0.0 0.7214782
## 19:     mean_imp      3  0.010   0.5 0.7212559
## 20:   median_imp      3  0.010   0.5 0.7212559
## 21:     mean_imp      3  0.001   1.0 0.7212559
## 22:   median_imp      3  0.001   1.0 0.7211448
## 23:     mean_imp      2  0.010   1.0 0.7209225
## 24:   median_imp      2  0.010   1.0 0.7209225
## 25:     mean_imp      2  0.001   0.0 0.7208113
## 26:   median_imp      3  0.010   0.0 0.7181995
## 27:     mean_imp      2  0.100   0.0 0.7179772
## 28:   median_imp      2  0.010   0.5 0.7176438
## 29:     mean_imp      3  0.010   1.0 0.7176438
## 30:   median_imp      3  0.010   1.0 0.7176438
## 31:   median_imp      2  0.001   0.5 0.7174215
## 32:     mean_imp      2  0.001   1.0 0.7174215
## 33:   median_imp      2  0.100   0.0 0.7145874
## 34:     mean_imp      2  0.010   0.5 0.7142540
## 35:     mean_imp      2  0.001   0.5 0.7140317
## 36:   median_imp      2  0.001   1.0 0.7140317
##     imp_strategy n_bins lambda alpha       acc

Проверим выбранную модель (первая строка из списка выше) на тестовых данных:

rec_obj <- create_recipe(data = train_data, 
                         imp_strategy = "mean_imp", 
                         n_bins = 3)

rec_obj <- prep(rec_obj, training = train_data, retain = FALSE)

train <- bake(rec_obj, 
              new_data = train_data,
              everything())
train_x <- 
  train %>%
  select(-churn) %>%
  as.matrix()
train_y <- train$churn
  
test <- bake(rec_obj, 
             new_data = test_data,
             everything())
test_x <- 
  test %>%
  select(-churn) %>%
  as.matrix()
test_y <- test$churn

model <- glmnet(x = train_x, 
                y = train_y, 
                family = "binomial",
                lambda = 0.1,
                alpha = 0)
  
preds <- predict(model,
                 type = "response",
                 newx = test_x)
preds <- ifelse(preds > 0.5, 1, 0)
preds <- factor(preds, levels = c(0, 1))
result <- accuracy(data.frame(preds = preds, test_y = test_y), 
                   preds, 
                   test_y)
result

## # A tibble: 1 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.704