Поговорим не про еду, а про разведочный анализ данных (exploratory data analysis, EDA) который является обязательной прелюдией перед любым суровым ML.

Будем честны, процесс довольно занудный, и чтобы выцепить хоть какие-то значимые инсайты про наши данные — требуется потратить достаточное количество времени активно используя любимую библиотеку визуализации.

А теперь представим что мы довольно ленивы (но любопытны) и будем следовать этому постулату всю эту статью.

Исходя из этого зададим себе вопрос: нет ли в природе такого хитрого инструмента который бы позволил просто нажать CTRL+ENTER в любимой IDE и вывести на одном лишь экране (без прокруток вниз и бесчисленных микроскопических фасетов) целостную картину с полезной информацией про наш датасет?

При этом держим в голове другую мысль — если такой инструмент и существует, то он не заменит классический EDA, но станет для нас сильным подспорьем для тех случаев когда не стоит часами упарываться по визуализации чтобы быстро подчеркнуть основные паттерны в наших данных.

Структура этой статьи:

1. Небольшой препроцессинг
2. Визуализация информативности предикторов
3. Дискретизация переменных
4. Correlationfunnel
5. Ranked Cross-Correlations
6. easyalluvial

Закончим с вводной и возьмем за основу практический пример.


Наиболее удобным в качестве примера мне показался Титаник, его размеры не слишком малы как Iris, он имеет малоинформативные переменные, он хорошо изучен и имеет понятные предикторы и что немаловажно — историческую основу.

Тем более, что я нашел на Хабре статью где автор провел достаточно скрупулезный EDA этого набора данных и на основе картинок продемонстрировал найденные выводы. Это будет свое рода наш Baseline.

Ссылка на статью с громким названием для нашего «Baseline_EDA»:
Титаник на Kaggle: вы не дочитаете этот пост до конца.

Чтобы не заморачиваться со скачиванием/чтением csv из сети сразу цепляем из CRAN оригинальный набор данных

install.packages("titanic") 
data("titanic_train",package="titanic")

Краткий препроцессинг

Данный пример настолько изъезжен в сети препроцессингом вдоль и поперек, что особо обсасывать эту тему не буду, делаю базовые вещи: извлекаю из имени гоноратив (титул) как важный предиктор, по нему делаю медианное заполнение пропусков в возрасте.

library(tidyverse)
titanic_train %>% str

d <- titanic_train %>% as_tibble %>%
  mutate(title=str_extract(Name,"\\w+\\.") %>% str_replace(fixed("."),"")) %>%
  mutate(title=case_when(title %in% c('Mlle','Ms')~'Miss', # нормализуем вариации
                         title=='Mme'~ 'Mrs',
                         title %in% c('Capt','Don','Major','Sir','Jonkheer', 'Col')~'Sir',
                         title %in% c('Dona', 'Lady', 'Countess')~'Lady',
                         TRUE~title)) %>%
  mutate(title=as_factor(title),
         Survived=factor(Survived,levels = c(0,1),labels=c("no","yes")),
         Sex=as_factor(Sex),
         Pclass=factor(Pclass,ordered = T)) %>%
  group_by(title) %>% # ниже - заполняем пропуски медианами по титулу
  mutate(Age=replace_na(Age,replace = median(Age,na.rm = T))) %>% ungroup

# посмотрим на распределение титулов по полу чтобы убедиться что все в порядке
table(d$title,d$Sex) 

Не все йогурты одинаково полезны…

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

Метрикой оценки «полезности» переменной являются freqRatio (соотношение частот самого популярного значения относительно второго значения по частотности) и percentUnique (мощность или cardinality — доля уникального числа значений от общего числа значений)
Детальную справку можно увидеть из пакета caret
?caret::nearZeroVar

(feat.scan <- caret::nearZeroVar(x = d,saveMetrics = T) %>% rownames_to_column("featName") %>% as_tibble)

Мне наиболее удобно мониторить переменные в двумерной плоскости (прологарифмировав обе оси чтобы не случился overplotting точек в одну маленькую кучу из-за точек-выбросов).
Никогда не задавался вопросом — является ли этот шаг EDA, но пока писал эту статью задумался: мы же сейчас проводим разведочный анализ некой полезности предикторов, их визуальную оценку, так почему ж это не EDA?

# install.packages("ggrepel")
library(ggrepel)
ggplot(feat.scan,aes(x=percentUnique,y=freqRatio,label=featName,col=featName))+ geom_point(size=2)+
  geom_text_repel(data = feat.scan,size=5)+scale_x_log10()+scale_y_log10()+theme_bw()

Малоинформативными считаем предикторы-выбросы либо по мощности (ось Х) либо по соотношению частот (ось Y) и соответственно откладываем в сторону:
PassengerId; Name; Ticket; Cabin

useless.feature <- c("PassengerId","Name","Ticket","Cabin")
d <- d %>% select_at(vars(-useless.feature))

Эта вселенная дискретна

Для того чтобы понимать как нижеперечисленные библиотеки подготоавливают данные — в этом разделе покажем на небольших примерах что происходит в этих библиотеках на этапе подготовки данных.

На первом шаге необходимо привести все данные к единому типу — часто данные в одном наборе могут быть и категориальными и числовыми, причем числа могут иметь выбросы а категориальные данные — редкие категории.

Для конвертации непрерывных переменных в категориальные можно разложить наши числа по бинам с определённым периодом дискретизации.

Простейший пример разложения на 5 бинов:

iris %>% as_tibble %>% mutate_if(is.numeric,.funs = ggplot2::cut_number,n=5)

Для получения силы и направленности взаимосвязей отдельных элементов среди предикторов используется второй прием - one hot encoding

library(recipes)
iris %>% as_tibble %>% mutate_if(is.numeric,cut_number,n=5) %>% 
  recipe(x = .) %>% step_dummy(all_nominal(),one_hot = T) %>%  prep %>% juice %>% glimpse

Вместо 5 предикторов у нас их теперь 23, зато бинарных:



В общем то на этом трюки по преобразованию заканчиваются, но с этих этапов как раз начинается работа 2-х из 3-х библиотек для нашего «неклассического» EDA.

Далее я знакомлю с функциональностью 3-х библиотек визуализации:

1. Correlationfunnel — показывает влияние отдельных значений предикторов на таргет (т.е. можно назвать это EDA supervized learning)
2. Lares — показывает влияние отдельных значений предикторов на другие отдельные значения других предикторов (т.е. можно назвать это EDA unsupervized learning)
3. easyalluvial — показывает совокупную взаимосвязь сгруппированных значений топ «X» предикторов на таргет (т.е. можно назвать это EDA supervized learning)

Видно что функциональность у них разная, поэтому, демонстрируя эти библиотеки, я буду цитировать выводы автора из статьи нашего «Baseline_EDA» в зависимости от вышеописанной функциональности этого пакета. (Например если автор показывает зависимость возраста на выживаемость то вставлю такую цитату в Correlationfunnel, если возраст от класса — то в Lares и т.д.)

На сцене первая библиотека.

correlationfunnel

correlationfunnel is to speed up Exploratory Data Analysis (EDA)


В виньетке к библиотеке неплохо описана методология, приведу фрагмент расчета корреляции по бинарным значениям



Библиотека предполагает наличие таргета (зависимой переменной) в наших данных и сразу на одной картинке показывает силу и направление взаимосвязи а так же ранжирует по убыванию этой силы образуя визуальную воронку (собственно отсюда и название).

Встроенные в библиотеку функции бинаризации позволяют так же сводить в Others малочисленные категории.

Так как библиотека не работает с целочисленными переменными — преобразуем их в numeric и вернемся к нашему Титанику.

#install.packages("correlationfunnel")
library(correlationfunnel)
d <- d %>% mutate_if(is.integer,as.numeric)
d %>% binarize(n_bins = 5,thresh_infreq = .02,one_hot = T) %>% # бинаризация втроенной функцией
  correlate(target = Survived__yes) %>% plot_correlation_funnel() # "interactive = T" - plotly!

По оси Х у нас сила и направление корреляции, по оси Y наши предкторы, отранжированные по убыванию. Первым сверху всегда отражается таргет т.к. у него самая сильная корреляция с самим собой (-1;1).

Давайте проверим как выводы по этому графику пересекаются с выводами автора нашего «Baseline_EDA».

Следующий график подтверждает теорию, что чем выше класс каюты пассажира — тем больше шансы выжить. (Под «выше»" я имею ввиду обратный порядок, т.к. первый класс выше чем второй и, тем более, третий.)

Воронка показывает что класс является третьим по силе корреляции предиктором и действительно у 3 класса обратная корреляция, у 1го — сильная положительная.

Сравним шансы выжить у мужчин и женщин. Данные подтверждают теорию, высказанную ранее.

(В целом, уже можно сказать, что основными факторами модели будет пол пассажира)

Воронка показывает, что пол пассажира является 2-м по степени корреляции, женский пол коррелирован с выживанием, мужской — со смертью.

Также можно проверить гипотезу, что выживают более молодые, т.к. они быстрее двигаются, лучше плавают и т.д.

Как видно, явная зависимость здесь не просматривается.

Воронка действительно говорит о слабой значимости этого предиктора (напомню что гоноратив/титул содержит в себе возраст именно поэтому возраст не столь значим), но даже тут воронка показывает что больше шансов выжить у категорий «минус бесконечность — 20 лет» (т.е. дети) и 30-38 (состоятельные люди, возможно 1 класс).

Давайте введём такой показатель как Процент выживаемости и посмотрим на его зависимость от групп, которые получились на предыдущем этапе

(группы у автора — имеется в виду титул).

Воронка полностью подтверждает найденные выводы автора

Теперь посмотрим на информацию, которую можно получить из количества родственников на корабле.

Очень похоже, что на выживаемость отрицательно влияет как отсутствие родственников, так и большое из количество.

SibSP в воронке явно говорит о том же.

Ну и конечно помимо выводов автора здесь можно увидеть и другие закономерности, удовольствие от созерцания оставлю читателю

Lares

Find Insights with Ranked Cross-Correlations



Автор данной библиотеки пошел еще дальше — он показывает зависимости не только на таргет, но и всех на все.

Ranked Cross-Correlations not only explains relationships of a specific target feature with the rest but the relationship of all values in your data in an easy to use and understand tabular format.

It automatically converts categorical columns into numerical with one hot encoding (1s and 0s) and other smart groupings such as “others” labels for not very frequent values and new features out of date features.

По ссылке выше можно посмотреть пример где автор скармливает датасет звездных войн своему пакету и демонстрирует найденные зависимости, я на его странице залип, оч.понравилось.

Попробуем на нашем примере.

# Осторожно, тянет довольно много зависимых пакетов:
# devtools::install_github("laresbernardo/lares")
library(lares)
corr_cross(df = d,top = 30)

Помимо пересечения с выводами на основе цитат в Correlationfunnell приведем отдельные цитаты которые здесь можем увидеть безотносительно таргета:

Также можно обнаружить и другие закономерности. Между возрастом и классом существует отрицательная корреляция, что, скорее всего, связано с более возрастные пассажиры чаще могли себе позволить более дорогую каюту.

В цитате выше автор делает такой вывод по корреляционному анализу 2-х полей в совокупности, у нас же с учетом One-Hot-Encoding это видно по сильной положительной корреляции между Age+P_Class_1.

Кроме того, стоимость билета и класс тесно связаны (высокий коэффициент корреляции), что вполне ожидаемо.

Третья строка сверху: Fare+P_Class_1

Помимо пересечения с выводами автора тут можно подчеркнуть так же много чего еще интересного, так же оставлю удовольствие от созерцания для читателя.

Помимо опционального выбора топ Х самых сильных инсайтов можно так же отразить всю картину и место этих значимых точек в общей массе

corr_cross(df = d,type=2)

easyalluvial

Data exploration with alluvial plots



Здесь автор так же как и в 2-х предыдущих пакетах выполняет на старте бинаризацию числовых переменных, однако дальше его пути с теми библиотеками расходятся: вместо {One-HotEncoding + корреляция} библиотека раскладывает топ Х самых интересных предикторов (пользователь решает сам — какие передать) по значениям, формируя потоки, цвет которых зависим от таргета, а ширина потока от числа наблюдений в этом потоке.

Числовые переменные раскладываются на категории HH (High High), MH(Medium High), M (Medium), ML (Medium Low), LL (Low Low)

Для начала возьмём наиболее значимые предикторы на основе графика из correlationfunnel:

cor.feat <- c("title","Sex","Pclass","Fare")

Далее делаем график

# install.packages("easyalluvial")
library(easyalluvial)

al <- d %>% select(Survived,cor.feat) %>% 
  alluvial_wide(fill_by = "first_variable")
add_marginal_histograms(p = al,data_input = d,keep_labels = F)

Для цитат автора перерисуем график с использованием соответствующих предикторов

cor.feat <- c("Sex","Pclass","Age")
al <- d %>% select(Survived,cor.feat) %>% 
  alluvial_wide(fill_by = "first_variable")
add_marginal_histograms(p = al,data_input = d,keep_labels = F)

Например, по следующему графику прекрасно видно, что основные группы выживших — это женщины первого и второго класса всех возрастов.

График помимо этого показывает что выжившие женщины 3 класса так же не малочисленная группа

А среди мужчин выжили все мальчики моложе 15 лет кроме третьего класса обслуживания и небольшая доля мужчин более старшего возраста и в основном из первого класса.

Сказанное подтверждается, но опять таки видим потоки выживших мужчин 3 класса в категории возраста LL, ML.

Все выше было про пакет «easyalluvial», однако автор написал второй пакет «parcats» который поверх plotly делает вышеприведенный график интерактивным (как в заголовке этого раздела).
Это дает возможность не только видеть tooltip-контекст, но и переориентировать потоки для лучшего визуального восприятия. (к сожалению пока библиотека не очень оптимизирована и на титанике у меня подтормаживает)

# install.packages("parcats")
library(parcats)
cor.feat <- c("title","Sex","Pclass","Fare")
a <- d %>% select(Survived,cor.feat) %>% 
  alluvial_wide(fill_by = "first_variable")
parcats(p = a,marginal_histograms = T,data_input = d)

Бонус

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

Ссылка на статью автора: Visualise model response with alluvial plots
Причем особенность в том что из всех библиотек, которые я видел, максимальное на одном графике объяснялся отклик черного ящика не более чем в 2 мерной системе координат (по одной на каждый предиктор), цветом объяснялся отклик.

Библиотека easyalluvial позволяет делать такое более чем по 2-м предикторам одновременно (лучше конечно не увлекаться).

Для примера обучим на нашем массиве данных случайный лес и отразим объяснение случайного леса по 3м предикторам.

library(ranger)
m <- ranger(formula = Survived~.,data = d,mtry = 6,min.node.size = 5, num.trees = 600,
            importance = "permutation")
library(easyalluvial)
(imp <- importance(m) %>% as.data.frame %>% easyalluvial::tidy_imp(imp = .,df=d)) # фрейм важности предикторов из модели
# генерим N-мерное пространство предикторов с комбинациями (к сожалению в том.числе и невозможными!) их значений
dspace <- get_data_space(df = d,imp,degree = 3) 
# получаем отклик по пространству
pred = predict(m, data = dspace)
alluvial_model_response(pred$predictions, dspace, imp, degree = 3)

Так же у автора есть коннектор к CARET-моделям (не знаю насколько это актуально сейчас учитывая tidymodels)

library(caret)
trc <- trainControl(method = "none")
m <- train(Survived~.,data = d,method="rf",trControl=trc,importance=T)
alluvial_model_response_caret(train = m,degree = 4,bins=5,stratum_label_size = 2.8) 

Заключение

Еще раз повторюсь что не призываю к замене классического EDA, но согласитесь — приятно когда есть альтернатива, позволяющая сэкономить кучу времени, особенно учитывая что человеки от природы достаточно ленивы, а это, как известно, двигатель прогресса :)