Мой канал - https://t.me/tidy_mng
Последнее время пакет tidymodels активно развивается в направлении задач машинного обучения.
Почему tidymodels? В R есть некоторая несогласованность между библиотеками, что является результатом разных подходов в их разработке.
Несколько лет назад Макс Кун разработал пакет caret, целью которого было создать единую платформу для моделей машинного обучения, существующих в R. Caret был прекрасен во многих отношениях, но далек от идеала. Но это был прекрасный старт. В связи с этим RStudio пригласила Макса Куна для разработки “аккуратной” версии данного пакета. В итоге, мы получили tidymodels.
Что такое tidymodels?
Tidymodels, как и tidyverse, состоит из нескольких библиотек. К основным можно отнести:
rsample— для разделения выборки;recipes— для предварительной обработки;parsnip— для определения модели;yardstick— для оценки модели;tune— для процедуры настройки параметров;workflows— для объединения предварительных процедур.
Чтобы установить весь набор библиотек, которые содержатся в tidymodels — выполните команду:
install.packages("tidymodels") #устанавливаем библеотеку
library(tidymodels) #запускаем библиотеку── Attaching packages ──────────────────────────────────────── tidymodels 1.0.0 ──
✔ broom 1.0.0 ✔ recipes 1.0.1
✔ dials 1.0.0 ✔ rsample 1.1.0
✔ dplyr 1.0.9 ✔ tibble 3.1.8
✔ ggplot2 3.3.6 ✔ tidyr 1.2.0
✔ infer 1.0.2 ✔ tune 1.0.0
✔ modeldata 1.0.0 ✔ workflows 1.0.0
✔ parsnip 1.0.0 ✔ workflowsets 1.0.0
✔ purrr 0.3.4 ✔ yardstick 1.0.0Набор данных
Мы будем использовать набор данных, в котором содержится информация о статусе диабета женщин племени Пима.
Информация о данных:
pregnant — кол-во беременностей;
glucose — концентрация глюкозы в крови;
pressure — артериальное давление;
triceps — толщина кожной складки трицепса;
insulin — 2-часовой сывороточный инсулин;
mass — индекс массы тела;
pedigree — функция родословной;
age — возраст;
diabetes — переменная класса (тест на диабет)
#Загрузим набор данных из библиотеки mlbench
library(mlbench)
data(PimaIndiansDiabetes)
#Для простоты работы, переименум наш набор данных
diabetes_orig <- as_tibble(PimaIndiansDiabetes)
diabetes_orig# A tibble: 768 × 9
pregnant glucose pressure triceps insulin mass pedigree age diabetes
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct>
1 6 148 72 35 0 33.6 0.627 50 pos
2 1 85 66 29 0 26.6 0.351 31 neg
3 8 183 64 0 0 23.3 0.672 32 pos
4 1 89 66 23 94 28.1 0.167 21 neg
5 0 137 40 35 168 43.1 2.29 33 pos
6 5 116 74 0 0 25.6 0.201 30 neg
7 3 78 50 32 88 31 0.248 26 pos
8 10 115 0 0 0 35.3 0.134 29 neg
9 2 197 70 45 543 30.5 0.158 53 pos
10 8 125 96 0 0 0 0.232 54 pos
# … with 758 more rowsПоверхностный анализ показывает, что в данных много нулей. Стоит отметить, что, например, индекс массы тела или толщина кожной складки трицепса не могут быть равны нулю. Это означает, что пропущенные значения в нашем наборе данных, записывались как “0”.
Построим гистограмму по triceps, чтобы убедиться в нашем предположении:
ggplot(diabetes_orig) +
geom_histogram(aes(x = triceps))
Эта проблема актуальная также для параметров: glucose, pressure, insulin, mass. Давайте заменим все нули в этих данных на NA.
diabetes_clean <- diabetes_orig %>%
mutate(across(c(glucose:mass), na_if, 0))# A tibble: 768 × 9
pregnant glucose pressure triceps insulin mass pedigree age diabetes
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct>
1 6 148 72 35 NA 33.6 0.627 50 pos
2 1 85 66 29 NA 26.6 0.351 31 neg
3 8 183 64 NA NA 23.3 0.672 32 pos
4 1 89 66 23 94 28.1 0.167 21 neg
5 0 137 40 35 168 43.1 2.29 33 pos
6 5 116 74 NA NA 25.6 0.201 30 neg
7 3 78 50 32 88 31 0.248 26 pos
8 10 115 NA NA NA 35.3 0.134 29 neg
9 2 197 70 45 543 30.5 0.158 53 pos
10 8 125 96 NA NA NA 0.232 54 pos
# … with 758 more rowsРазделим данные на тестовую и обучающую выборки
Для начала, разделим наш набор данных на обучающую и тестовые выборки. Обучающая выборка будет использоваться для подгонки нашей модели и настройки ее параметров, а тестовая — для оценки производительности нашей модели.
Это разделение мы можем выполнить автоматически с помощью функции inital_split() (из пакета rsample), которая создает специальный объект “split”.
set.seed(234589)
#разделим на обучение (75%) и на тест (25%)
diabetes_split <- initial_split(diabetes_clean,
prop = 3/4)
diabetes_split<Training/Testing/Total>
<576/192/768>Выводимый результат diabetes_split, сообщает нам, сколько наблюдений у нас есть в обучающем наборе, тестовом наборе и сколько в общем наблюдений: <train/test/total>.
Обучающие и тестовые наборы могут быть извлечены из объекта diabetes_split с помощью функций training() и testing(). Хотя, на самом деле, мы не будем использовать эти объекты в конвейере (нам понадобится сам объект diabetes_split).
#извлекам обучающие и тестовые наборы
diabetes_train <- training(diabetes_split)
diabetes_test <- testing(diabetes_split)В дальнейшем мы захотим выполнить некоторые настройки параметров. Для этого нам понадобится метод кросс-валидации. Для создания такого набора данных мы воспользуемся функцией vfold_cv().
#создадим объект кросс-валидации из обучающего набора
diabetes_cv <- vfold_cv(diabetes_train)
diabetes_cv# 10-fold cross-validation
# A tibble: 10 × 2
splits id
<list> <chr>
1 <split [518/58]> Fold01
2 <split [518/58]> Fold02
3 <split [518/58]> Fold03
4 <split [518/58]> Fold04
5 <split [518/58]> Fold05
6 <split [518/58]> Fold06
7 <split [519/57]> Fold07
8 <split [519/57]> Fold08
9 <split [519/57]> Fold09
10 <split [519/57]> Fold10Определим “рецепт”
"Рецепты" позволяют указать роль каждой переменной в качестве результата или предиктора, а также любые шаги предварительной обработки, которые мы хотим провести (например, нормализация, PCA и т. д.).
Создание “рецепта” состоит из двух частей:
Указать формулу (
recipe()): необходимо указать переменную результата и переменные предиктора.Указать шаги предварительной обработки (
step_zzz()): такие как: работа с пропущенными значениями, создание фиктивных переменных, масштабирование и многое другое.
Например, мы можем сделать следующий “рецепт”:
#определяем recipe
diabetes_recipe <-
#который состоит из формулы (исход ~ предикторы)
recipe(diabetes ~ pregnant + glucose + pressure + triceps +
insulin + mass + pedigree + age,
data = diabetes_clean) %>%
#и некоторые этапы предварительной обработки
step_normalize(all_numeric()) %>%
step_impute_knn(all_predictors())Иные функции пакета вы можете найти здесь.
В приведённом выше примере мы использовали функции all_numeric() и all_predictors() в качестве аргументов предварительной обработки. Такие функции называются - «ролевыми селекторами» и указывают, что мы хотим применить шаг ко "всем числовым" переменным или "всем предикторным переменным". Список всех потенциальных селекторов ролей можно найти, набрав ?selections в консоли.
Обратите внимание, что мы работаем с исходным объектом данных diabetes_clean, а не diabetes_train или diabetes_split. В целом, мы могли бы использовать любой из них, так как все, что "рецепт" "берет" - это имена и роли переменных результата и предиктора. Позже мы применим этот рецепт к конкретным наборам данных.
Если мы распечатаем сводку по объекту diabetes_recipe, то убедимся в выше сказанном:
diabetes_recipeRecipe
Inputs:
role #variables
outcome 1
predictor 8
Operations:
Centering and scaling for all_numeric()
K-nearest neighbor imputation for all_predictors()Если вы хотите извлечь предварительно обработанный набор данных, то вы можете сначала использовать prep(), а на следующем шаге, при помощи функции juice(), извлечь результат. На самом деле данный шаг не обязателен, так как он будет выполнен "под капотом" при подгонке модели.
diabetes_train_preprocessed <- diabetes_recipe %>%
# применяем рецепт к обучающим данным
prep(diabetes_train) %>%
# извлекаем предварительно обработанный набор данных
juice()
diabetes_train_preprocessedA tibble: 576 × 9
pregnant glucose pressure triceps insulin mass pedigree age diabetes
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct>
1 1.23 -0.390 0.262 0.892 -0.317 -0.656 0.502 -0.201 pos
2 0.0447 1.09 -0.0616 -0.0348 -0.219 -0.168 -0.400 0.301 neg
3 1.82 1.91 -0.224 0.595 0.146 0.379 -0.813 0.301 neg
4 -1.14 -0.0948 -0.710 -0.591 -0.522 -0.311 3.76 -1.04 neg
5 -0.548 0.233 0.424 0.706 0.239 1.56 2.25 -0.201 pos
6 0.637 -0.0620 -1.84 -0.683 0.190 -0.771 2.53 -0.0335 pos
7 -0.844 -1.64 -2.01 -1.05 -0.629 -1.73 -0.445 -0.953 neg
8 -0.548 -0.423 -1.68 -0.313 -0.735 0.00505 -0.460 -0.953 neg
9 -1.14 0.463 -0.386 1.17 0.796 1.41 -0.320 -0.786 pos
10 1.53 1.41 0.424 0.150 0.877 0.0482 -0.968 0.969 pos
# … with 566 more rowsУкажем модель
Пока, мы только разделили наши данные на тестовую и обучающую выборки и указали шаги предварительной обработки данных при помощи пакета recipe. Следующее, что мы хотим сделать — нашу модель. Для этого мы воспользуемся пакетом parsnip.
Parsnip предлагает унифицированный интерфейс для огромного количества моделей, существующих в R. Это значит, что мы сможем при помощи одной строки кода строить различные модели, придерживаясь одной логики.
Есть несколько основных компонентов, которые нам понадобятся для построения модели:
Тип модели: какую модель мы хотим. Для этого необходимо выбрать одну из нескольких функций, например:
rand_forest()для случайного леса,logistic_reg()для логистической регрессии,svm_poly()для сингулярное разложения. С полным перечнем функций можно ознакомиться здесь.Аргументы: значения параметров модели (указывааются последовательно в разных моделях), задаем с помощью
set_args().Движок: базовый пакет, из которого должна быть получена модель (например,
rangerдля реализации случайного леса), задаем с помощьюset_engine().Режим: тип прогнозирования, так как различные пакеты могут выполнять как классификацию, так и регрессию. Задаем с помощью
set_mode().
Например, если мы хотим подогнать модель случайного леса, реализованную пакетом ranger, для целей классификации, и хотим указать параметр mtry (количество случайно выбранных переменных, которые будут учитываться при каждом разбиении деревьев), то нам нужно задать следующую спецификацию модели:
rf_model <-
# указываем, что модель представляет собой случайный лес
rand_forest() %>%
# указываем `mtry`
set_args(mtry = tune()) %>%
# выберите движок/пакет, лежащий в основе модели
set_engine("ranger", importance = "impurity") %>%
# выбрать режим регрессии или классификации
set_mode("classification") Если вы хотите в дальнейшем посмотреть на оценки важности переменных в вашей модели, необходимо будет воспользоваться аргументом importance в set_engine. Для пакета ranger важными вариантами являются impurity или permutation.
В качестве другого примера воспользуемся пакетом glm для логистической регрессии:
lr_model <-
# указать, что что модель - логистическая регрессия
logistic_reg() %>%
# выбрать пакет, который лежит в основе модели
set_engine("glm") %>%
# выбрать режим регрессии или классификации
set_mode("classification") Обратите внимание, что данный код плохо подходит для конечной модели. Во-первых, как и recipe, он просто описывает модель. Во-вторых, в этой спецификации модели ничего не относится к набору данных о диабете.
Кроме того, установка параметра tune() означает, что этот параметр будет настроен позже на этапе настройки конвейера (то есть, будет выбрано значение параметров, обеспечивающих наилучшую производительность). Но вы можете настроить этот параметр сами, если хотите контролировать модель, например, сделав это так: set_args(mtry = 4).
Соберем все в один рабочий процесс
Теперь мы готовы объединить модель и “рецепты” в один рабочий процесс. Это вы можете сделать при помощи workflow() из пакета workflows , куда передадите данные о модели и “рецепте”.
# установливаем рабочий процесс
rf_workflow <- workflow() %>%
# добавляем "рецепт"
add_recipe(diabetes_recipe) %>%
# добавляем модель
add_model(rf_model)Обратите внимание, что мы еще не реализовали этапы предварительной обработки в “рецепте” и не подогнали модель. Мы написали основу. Только после того, как мы настроим параметры или подгоним модель, мы сможем получить результаты.
Настройка параметров
Поскольку у нас есть параметр mtry, нам необходимо настроить его (т. е. выбрать значение, обеспечивающее наилучшую производительность). Если у вас нет параметров для настройки, вы можете пропустить этот шаг.
Обратите внимание, что для настройки мы будем использовать объект кросс-валидации diabetes_cv. Для этого мы указываем диапазон значений mtry, которые мы хотим попробовать, а затем добавляем параметр tune_grid() в наш рабочий процесс, используя пакет tune. Нас интересует два показателя: accuracy и roc_auc (пакет yardstick).
# укажите, какие значения вы хотите попробовать
rf_grid <- expand.grid(mtry = c(3, 4, 5))
# извлекаем результаты
rf_tune_results <- rf_workflow %>%
tune_grid(resamples = diabetes_cv, # объект кросс-валидации
grid = rf_grid, # сетка значений, котрую мы пробуем
metrics = metric_set(accuracy, roc_auc) # показатели, которые нас интересуют
)Вы можете настроить несколько параметров одновременно, использовав функцию expand.grid(). Например, expand.grid(mtry = c(3, 4, 5), trees = c(100, 500))
Всегда интересно посмотреть на результаты кросс-валидации. collect_metrics() — удобная функция, которую можно использовать для извлечения метрик, рассчитанных для объекта.
# посмотрим на результат
rf_tune_results %>%
collect_metrics()# A tibble: 6 × 7
mtry .metric .estimator mean n std_err .config
<dbl> <chr> <chr> <dbl> <int> <dbl> <chr>
1 3 accuracy binary 0.762 10 0.0108 Preprocessor1_Model1
2 3 roc_auc binary 0.840 10 0.0116 Preprocessor1_Model1
3 4 accuracy binary 0.767 10 0.0130 Preprocessor1_Model2
4 4 roc_auc binary 0.840 10 0.0128 Preprocessor1_Model2
5 5 accuracy binary 0.766 10 0.0108 Preprocessor1_Model3
6 5 roc_auc binary 0.837 10 0.0118 Preprocessor1_Model3И accuracy, и roc_auc, при mtry = 4 дают наилучшие результаты.
Завершение рабочего процесса
Мы хотим добавить в наш рабочий процесс слой, который соответствует настроенному параметру, т.е. установить такой mtry , который дал наилучшие результаты. Если вы не настроили никаких параметров, вы можете пропустить этот шаг.
Мы можем извлечь наилучшее значение для метрик точности, применив функцию select_best() к объекту tune.
param_final <- rf_tune_results %>%
select_best(metric = "accuracy")
param_final# A tibble: 1 × 2
mtry .config
<dbl> <chr>
1 4 Preprocessor1_Model2Затем мы можем добавить этот параметр в рабочий процесс с помощью функции finalize_workflow().
rf_workflow <- rf_workflow %>%
finalize_workflow(param_final)Оценим модель на тестовом наборе
Теперь, когда мы определили наш рецепт, нашу модель и настроили параметры модели, мы готовы приступить к подгонке окончательной модели.
Поскольку вся информация содержится в объекте рабочего процесса (rf_workflow), мы применим функцию last_fit() к rf_workflow в рамках разделенного объекта diabetes_split. Это позволит автоматически обучить модель, указанную в рабочем процессе, используя данные для обучения и произвести оценку на основе тестового набора.
rf_fit <- rf_workflow %>%
last_fit(diabetes_split)
rf_fit# Resampling results
# Manual resampling
# A tibble: 1 × 6
splits id .metrics .notes .predictions .workflow
<list> <chr> <list> <list> <list> <list>
1 <split [576/192]> train/test split <tibble> <tibble> <tibble> <workflow>Обратите внимание, что создаваемый объект fit является объектом tibble со списком столбцов.
Эта замечательная особенность tidymodels упрощает взаимодействие с tidyverse. Хотя для реального использования этой связи требуется умение работать с purrr. Если вы хотите избежать пакета purrr, то существуют функции, которые позволяют достаточно легко извлекать необходимую нам информацию.
Так как мы используем объект diabetes_split, то метрики оцениваются только в тестовом наборе. Теперь, когда мы применим функцию collect_metrics() (напомним, что мы использовали ее при настройке параметров), она извлечет производительность конечной модели (поскольку rf_fit теперь состоит из одной конечной модели), примененной к тестовому набору.
test_performance <- rf_fit %>% collect_metrics()
test_performance# A tibble: 2 × 4
.metric .estimator .estimate .config
<chr> <chr> <dbl> <chr>
1 accuracy binary 0.792 Preprocessor1_Model1
2 roc_auc binary 0.853 Preprocessor1_Model1В целом производительность очень хорошая: accuracy 0.79 и AUC 0.85.
Вы также можете извлечь прогнозы тестового набора с помощью функции collect_predictions(). Обратите внимание, что в объекте прогнозов — 192 строки, которые соответствуют количеству наблюдений в тестовом наборе (это некоторое доказательство того, что они основаны на тестовом наборе, а не на тренировочном).
test_predictions <- rf_fit %>% collect_predictions()
test_predictions# A tibble: 192 × 7
id .pred_neg .pred_pos .row .pred_class diabetes .config
<chr> <dbl> <dbl> <int> <fct> <fct> <chr>
1 train/test split 0.365 0.635 3 pos pos Preprocessor1_…
2 train/test split 0.213 0.787 9 pos pos Preprocessor1_…
3 train/test split 0.765 0.235 11 neg neg Preprocessor1_…
4 train/test split 0.511 0.489 13 neg neg Preprocessor1_…
5 train/test split 0.594 0.406 18 neg pos Preprocessor1_…
6 train/test split 0.356 0.644 26 pos pos Preprocessor1_…
7 train/test split 0.209 0.791 32 pos pos Preprocessor1_…
8 train/test split 0.751 0.249 50 neg neg Preprocessor1_…
9 train/test split 0.961 0.0385 53 neg neg Preprocessor1_…
10 train/test split 0.189 0.811 55 pos neg Preprocessor1_…
# … with 182 more rowsПоскольку это всего лишь обычный tibble, мы можем сгенерировать матрицу ошибок.
test_predictions %>%
conf_mat(truth = diabetes, estimate = .pred_class) Truth
Prediction neg pos
neg 107 17
pos 23 45Мы также можем построить график распределения предсказанных вероятностей для каждого класса.
test_predictions %>%
ggplot() +
geom_density(aes(x = .pred_pos, fill = diabetes),
alpha = 0.5)
Вы можете использовать пакет purrr для извлечения прогнозов с помощью функции pull(). Следующий код делает почти то же самое, что и collect_predictions().
test_predictions <- rf_fit %>% pull(.predictions)
test_predictions# A tibble: 192 × 6
.pred_neg .pred_pos .row .pred_class diabetes .config
<dbl> <dbl> <int> <fct> <fct> <chr>
1 0.365 0.635 3 pos pos Preprocessor1_Model1
2 0.213 0.787 9 pos pos Preprocessor1_Model1
3 0.765 0.235 11 neg neg Preprocessor1_Model1
4 0.511 0.489 13 neg neg Preprocessor1_Model1
5 0.594 0.406 18 neg pos Preprocessor1_Model1
6 0.356 0.644 26 pos pos Preprocessor1_Model1
7 0.209 0.791 32 pos pos Preprocessor1_Model1
8 0.751 0.249 50 neg neg Preprocessor1_Model1
9 0.961 0.0385 53 neg neg Preprocessor1_Model1
10 0.189 0.811 55 pos neg Preprocessor1_Model1
# … with 182 more rowsУстановка и использование окончательной модели
В предыдущем разделе оценивалась модель, обученная на тренировочных данных. Но как только вы определили свою окончательную модель, вы хотите обучить ее на своем полном наборе данных, а затем использовать ее для прогнозирования. Для этого вам нужно применить функцию fit() в вашем рабочем процессе.
final_model <- fit(rf_workflow, diabetes_clean)Объект final_model содержит несколько вещей, включая объект ranger, обученный с параметрами, установленными с помощью рабочего процесса, содержащегося в rf_workflow, на основе данных в diabetes_clean
final_model══ Workflow [trained] ════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: rand_forest()
── Preprocessor ──────────────────────────────────────────────────────────────────
2 Recipe Steps
• step_normalize()
• step_impute_knn()
── Model ─────────────────────────────────────────────────────────────────────────
Ranger result
Call:
ranger::ranger(x = maybe_data_frame(x), y = y, mtry = min_cols(~4, x), importance = ~"impurity", num.threads = 1, verbose = FALSE, seed = sample.int(10^5, 1), probability = TRUE)
Type: Probability estimation
Number of trees: 500
Sample size: 768
Number of independent variables: 8
Mtry: 4
Target node size: 10
Variable importance mode: impurity
Splitrule: gini
OOB prediction error (Brier s.): 0.1583874 Если бы мы хотели предсказать статус диабета для новой женщины из племени Пима, мы бы обратились к функции predict().
Создадим пример, для новой женщины.
new_woman <- tribble(~pregnant, ~glucose, ~pressure, ~triceps, ~insulin, ~mass, ~pedigree, ~age,
2, 95, 70, 31, 102, 28.2, 0.67, 47)
new_woman# A tibble: 1 × 8
pregnant glucose pressure triceps insulin mass pedigree age
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 2 95 70 31 102 28.2 0.67 47Прогнозируемый диабетический статус для новой женщины - “отрицательный”.
predict(final_model, new_data = new_woman)# A tibble: 1 × 1
.pred_class
<fct>
1 neg Оценка важности переменных
Если вы хотите извлечь оценки важности из своей модели, вам нужно применить функцию extract_fit_parsnip(), а затем извлечь объект fit, который содержится в выводе:
ranger_obj <- extract_fit_parsnip(final_model)$fit
ranger_objRanger result
Call:
ranger::ranger(x = maybe_data_frame(x), y = y, mtry = min_cols(~4, x), importance = ~"impurity", num.threads = 1, verbose = FALSE, seed = sample.int(10^5, 1), probability = TRUE)
Type: Probability estimation
Number of trees: 500
Sample size: 768
Number of independent variables: 8
Mtry: 4
Target node size: 10
Variable importance mode: impurity
Splitrule: gini
OOB prediction error (Brier s.): 0.1583874 Затем вы можете извлечь оценки важности при помощи объекта $variable.importance
ranger_obj$variable.importancepregnant glucose pressure triceps insulin mass pedigree age
16.40687 79.68408 17.08361 22.10685 52.27195 42.60717 30.12246 33.19040 