Регрессиада. Часть 2. Метрики лгут

В прошлой статье этой серии [1] мы познакомились с семейством моделей регрессий, уделив особое внимание простой линейной регрессии. Сегодня мы продолжим наше знакомство, перейдя к множественной линейной регрессии. Но, как и всегда, мы делаем это не ради самого инструмента, а ради понимания его смысла и бизнес-применения.

Логика, лежащая в основе этой статьи, заимствована из моей рабочей практики. Она основана на ситуации, когда HR-метрики могут вводить в заблуждение, но глубокая аналитика позволяет избежать этого.

Введение

Руководство вымышленной компании «Линейные уравнения» ищет ответ на вопрос о том какие факторы влияют на эффективность сотрудников. Существует распространённое мнение среди менеджеров и HR, что выпускники столичных и сибирских ВУЗов отличаются по уровню эффективности. Известно, что всех действующих сотрудников можно разделить на две категории по этим группам ВУЗов. Менеджмент компании просит команду HR-аналитиков построить дашборд, в котором будет доступна динамика эффективности сотрудников по месяцам с разбивкой на ВУЗы, которые были ими окончены.

HR-аналитики подготавливают график с метрикой для дашборда. Также график дополняется индикаторами, согласно которым, средняя годовая эффективность выпускников московских ВУЗов составила 93%, в то время как их коллег из Сибири — 83%.

График с HR-дашборда компании "Линейные уравнения"

Внимательно посмотрите на график с метрикой. Посмотрели? Теперь ответьте на вопрос: какие выпускники более эффективные? Думаю, что большая часть читателей ответит также как менеджмент компании «Линейные уравнения», что более эффективны выпускники Московских ВУЗов. Так и в компании «Линейные уравнения», проанализировав этот график готовы к созданию стратегии привлечения новых сотрудников исключительно из московских ВУЗов.

Но есть одна проблема – это полностью неверный вывод. И вот почему: потому что я создал эти данные и мне известно какой вывод верный. Для анализа этой проблемы далее мы сфокусируемся на множественной линейной регрессии – инструменте, который позволит нам не попасть в ловушку. Мы рассмотрим её математические детали под капотом, научимся строить модели в R и сделаем на основе результата единственной верный вывод. Давайте начинать!

Совсем немного математики

Множественная линейная регрессия, с которой мы познакомимся и будем работать в этой статье, очень похожа на простую [1], с точки зрения формулы мы всего лишь добавляем больше переменных в модель и рассчитываем бета-коэффициенты для каждой.

Что меняется, точнее усложняется, так это интерпретация коэффициентов. Эндрю Гельман выделяет два вида интерпретации бета-коэффициентов в зависимости от цели [2]:

Предиктивная интерпретация – как в среднем изменяется зависимая переменная при сравнении двух объектов, которые отличаются на 1 единицу по соответствующей переменной, при условии, что остальные переменные неизменны.
Контрфактическая интерпретация - как в среднем изменяется зависимая переменная на уровне одного объекта, если изменить соответствующую переменную на 1 единицу по соответствующей переменной, при условии, что остальные переменные неизменны. Этот вид интерпретации возникает в причинно-следственном анализе, про который мы поговорим в четвертой статье.

В этом же блоке давайте поговорим ещё и про то, как обрабатывать в моделях нечисловые переменные. К примеру, как складывать или умножать переменную, в которой содержится текстовое название группы ВУЗа? Существуют разные способы как обработать такой случай, мы рассмотрим два:

переменные-индикаторы (дамми-переменные) – базовый уровень переменной исключается, остальные уровни принимают значение 0 или 1. Иными словами, количество бета-коэффициентов = количество категорий переменной - 1. Значения переменной напрямую участвуют в выражении регрессии. Классические регрессии используют этот вариант.
индексные переменные – уровни переменной раскладываются на целые числа от 1 до + ∞, но при этом напрямую не участвуют в выражении регрессии, а используются как количество уровней переменной, по которым происходит расчёт бета-коэффициентов. В этом случае количество бета-коэффициентов = количеству уровней переменной. Активно применяются в байесовских регрессиях по фреймворку МакЭлрита [3].

Давайте для понимая концепции визуализируем её на примере нашей переменной с ВУЗами. Если переменная ВУЗ = [“Московский ВУЗ”, “Сибирский ВУЗ”], тогда:

Практика

Процесс генерации данных

Мы подходим к теме, которой обычно не уделяется достаточно внимания. На психологическом факультете, где я проходил курс по статистике, как и позднее на экономическом, этот аспект, например, не рассматривался. Обычно практика аналитика, даже того, который строит регрессии, выглядит так: берутся имеющиеся данные и сразу же строится модель. Но этот шаг точно не должен быть первым, честно говоря, первый шаг мы обсудим в завершающей четвертой статье серии, сейчас же мы поговорим про шаг под номером два – это процесс генерации данных.

Прежде чем использовать какой-либо аналитический инструмент, необходимо понять закономерности, по которым были сгенерированы данные. На основе этого предположения мы проводим симуляцию и анализируем результаты модели. Если модель верно определяет заложенные в данных природные закономерности, тогда и только тогда её можно применять к реальным данным. В своей книге Ричард МакЭлрит [3] ярко иллюстрирует идею на примере послеэффектного смещения. Мы же воспользуемся этим подходом для того, чтобы показать, как более вдумчивая аналитика, позволяет избежать ошибочного вывода, к которому могут подвести метрики.

Давайте сгенерируем данные в R, которые рассматривали на графике выше. Хитрость кейса состоит в том, что эффективность зависит от стажа, а вовсе не от ВУЗа, который был окончен. Однако сотрудники, принятые из разных учебных заведений, имеют разный средний стаж, то есть мы искусственно создаем ситуацию конфаундера под названием труба [4].

# Устанавливаем зерно для воспроизводимости результатов
set.seed(753)
 
# Количество наблюдений
n <- 24

# Месяцы для каждой группы
months <- rep(1:(n/2), times=2)

# Группы ВУЗов
groups <- rep(c("Сибирский ВУЗ", "Московский ВУЗ"), each=(n/2))


# Генерируем разный опыт для каждой группы на основе нормального распределения
experience <- c(rnorm(n/2, mean=3, sd=1.5), # для сибирских ВУЗов
           		rnorm(n/2, mean=6, sd=1.5)) # для московских ВУЗов

# Генерируем эффективность, 
# где интерсепт = 75, бета-коэффициент для опыта = 3 и добавляем ошибку
efficiency <- rnorm(n, 75 + 3 * experience, 2.5) + rnorm(n, 0, 1)

# Соединяем все данные в один датафрейм
data <- data.frame(months, groups, experience, efficiency)

Построение модели

Помимо аналитиков, которые построили дашборд с HR-метрикой по техзаданию, в компании «Линейные уравнения» были и аналитики, предпочитающие думать над задачами. И первой их мыслью было проверить влияние возможных конфаундеров. Аналитики уже знали по предыдущим исследованиям, что на эффективность в их компании одно из сильнейших воздействий оказывает стаж. В то время как линейную ассоциацию между парами переменных мы можем обнаружить с помощью корреляции, то для случая, где у нас более двух переменных она не подходит, но зато подходит множественная линейная регрессия.

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

# Подключаем библиотеку psych, которая содержит один из вариантов графика
library(psych)

# Исключаем переменную с номером месяца и строим график
subset(data, select = -months) |> 
  pairs.panels(hist.col = "#ffa600")

На главной диагонали распределение переменной, ниже диагонали графики рассеивания по парным сочетаниям переменных, выше диагонали коэффициенты линейной корреляции между парами

Мы видим, что корреляция тоже завела бы нас не туда, показывая, что коэффициент между группой ВУЗа и эффективностью = -0.76. Также нам видна ключевая корреляция между опытом и эффективностью равная = 0.91.

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

Также практикой, которая вызывает осуждение, является оценка распределения по гистограммам. Ричард МакЭлрит даже придумал для этого термин «гистомантия», у меня был пост с переводом этого [5] выдающегося отрывка из его книги [3].

Столь же спорны статистические критерии проверки нормальности, такие как тест Шапиро-Уилка [6], так как на больших выборках они всегда будут утверждать, что распределение не нормальное. Ещё можно полагаться на оценку эксцесса и асимметрии, сам я предпочитаю именно этот метод, но в нашем случае он тоже ничего хорошего не дает.

# Библиотека для анализа распределений
library(fitdistrplus) 

# Запуск функции
descdist(data$efficiency, boot = 1000)

x = квадарт ассиметрии, y = эксцесс. Синяя точка отражает наблюдаемые данные, желтые точки изъяты бутстрэпом

График говорит, что по наблюдаемым данным у нас бета распределение, если делать бутстрэп [7], то тут возможно, что оно: бета, нормальное, равномерное, логнормальное или гамма… Но мы будем исходить из того, что оно нормальное и далее проверка остатков это подтвердит.

Переходим к построению множественной линейной регрессии, где эффективность ~ группа ВУЗа + стаж.

# Строим модель
model <- lm(efficiency ~ groups + experience, data=data)

# Анализируем результат
summary(model)


Call:
lm(formula = efficiency ~ groups + experience, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.5788 -1.7122 -0.4045  1.2251  5.4536 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          75.6790     3.5350  21.409 9.57e-16 ***
groupsСибирский ВУЗ  -1.0911     2.1592  -0.505    0.619    
experience            2.7890     0.5524   5.049 5.34e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.599 on 21 degrees of freedom
Multiple R-squared:  0.8562,	Adjusted R-squared:  0.8425 
F-statistic: 62.52 on 2 and 21 DF,  p-value: 1.432e-09

Во- первых, мы видим, что интерсепт = 76, что близко к исходному значению 75, которое мы сами задали. Бета-коэффициент переменной experience= 2.8, что также близко к ожидаемому значению 3. Также заметим, что наша модель объясняет 84% дисперсии.

Теперь мы смотрим на бета-коэффициент дамми-переменной «Сибирский ВУЗу», которая = -1.1, что можно было бы интерпретировать, как то, что в среднем уровень эффективности по-сибирски на 1.1% процента ниже, чем по-московски (именно так сравниваются дамми-переменные – относительно своего базового уровня). Но мы и этого с вами не сделаем, так как бросаем взгляд на колонку Pr(>|t|) и видим, что фактор ВУЗа статистически незначим!

Аналитики уже готовы оформить результат и рассказать бизнесу, что фактор ВУЗа никак не влияет на эффективность сотрудников, в том время как фактор опыта статистически значим и за каждый месяц опыта, уровень эффективности в среднем возрастает на 2.8%. Но поскольку аналитики были думающие, они понимали, что прежде, чем делать выводы – требуется провести диагностику модели.

Сдаем анализы

Как известно ничего хорошего не выйдет, если просто изобретать какие-либо девайсы, писать сайты, выдумывать новые алгоритмы, считать метрики и не проверять всё это. Тоже самое и с моделями регрессии – они тоже требуют проверки. Но какой? Список обязательных проверок зависит от вида модели, который используется. Для случая линейных моделей и вашего удобства, я свел всю основную информацию по этому вопросу в следующую таблицу:

Приступаем к проверкам:

1. Линейность. Для её оценки достаточно подняться по тексту и посмотреть на пэйрплот, точнее на нижнюю часть, в которой есть парные графики рассеивания. Мы видим, что наблюдаемые данные вполне описываются простой линей.

2. Независимость остатков, проверим критерием Дарбина-Уотсона [8]. Здесь все в порядке, p-value > 0.05, двигаемся дальше.

# Библиотека, содержащая тест
library(lmtest)

# Тест Дарбина-Уотсона
dwtest(model)


	Durbin-Watson test

data:  efficiency_lm
DW = 2.2949, p-value = 0.7437
alternative hypothesis: true autocorrelation is greater than 0

3. Нормальное распределение остатков проверим как графически, так и с помощью критерия. И убедимся, что все в полном порядке и по этому пункту.

# Библиотека, содержащая график QQ 
library(car)

# Квантильная проверка на нормальное распределениие
qqPlot(model$residuals)

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

# Критерий Шапиро-Уилка
shapiro.test(model$residuals)

	Shapiro-Wilk normality test

data:  efficiency_lm$residuals
W = 0.97174, p-value = 0.7099

4. Равномерная дисперсия. Проверим её с помощью критерия Бройша-Пагана [9].

# Вызов критерия Бройша-Пагана
bptest(model)

	studentized Breusch-Pagan test

data:  efficiency_lm
BP = 6.7284, df = 2, p-value = 0.03459

Статистика показывает нам, что значения p-value < 0.05, следовательно это требование нарушено. У нас с вами всего 24 наблюдения, предполагаем, что скорее всего есть какое-то влиятельное значение вызывающие эту ситуацию. Для поиска таких значений существует специальная метрика, называемая расстоянием Кука [10]. Строится эта метрика для каждого наблюдения.

library(ggplot2)

# Создадим датафрейм с метрикой
cooks_distances = data.frame(cd = cooks.distance(model),
                             n = 1:n)

# Визулизируем
ggplot(cooks_distances, aes(x=n, y=cd)) + 
  geom_col(fill = "#a8a8a8") +
  labs(title = "Расстояние Кука",
       x = "Номер наблюдения",
       y = "Значение") +
  scale_x_continuous(breaks = seq(1, n, by = 1)) +
  theme_classic() +
  theme(legend.position = "bottom",  
        legend.title = element_blank(),
        text = element_text(size = 12),
        legend.text = element_text(size = 12))

На норму этой метрики можно посмотреть по-разному, мы же просто исключим самое влиятельное относительно других наблюдение под номером 14. Если после этого пересчитать модель, то мы успешно пройдем данную проверку, при этом сам результат модели и вывод останутся прежними. В тексте я не даю пересчёта, но вы можете проверить самостоятельно.

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

# Рассчитаем коэффициент инфляции дисперсии
vif(model)

    groups experience 
  2.928427   2.928427

Прогнозы

Построить модель, проанализировать её результаты и проверить качество это всё ещё недостаточный набор действий. Так как самое важное это понять, какие прогнозы мы сможем делать на основе неё. Для этого тоже существуют разные подходы, мы пойдем следующим путем: создадим искусственный датасет, который будет содержать все возможные сочетания стажа и группы ВУЗов, сделаем на его основе прогноз, визуализируем результат.

# Создадим фейковый датасет, в котором все сочетания ВУЗов и стажа
fake_data <- data.frame(groups = rep(c("Сибирский ВУЗ", "Московский ВУЗ"), each=(10)), 
                        experience = rep(seq(1, 10, by = 1),2))

# Посчитаем прогноз средней эффективности на основе модели
fake_data$efficiency <- predict(model, fake_data)


# Визуализируем с помощью ggplot
ggplot(data = fake_data, aes(x=groups, y=efficiency, color=groups, group=groups)) +
  stat_boxplot(geom ='errorbar', color = c("#ffa600", "#a8a8a8"), size=2) + 
  geom_boxplot(size=2) +
  scale_color_manual(values = c("#ffa600", "#a8a8a8")) +
  labs(title = "Прогноз средней эффективности",
       x = "Группа ВУЗа",
       y = "Эффективность (%)") +
  theme_classic() +
  theme(legend.position="none",
        text = element_text(size = 12),
        legend.text = element_text(size = 12))

Прогноз показывает ведущую роль стажа, различие в 1% по группам ВУЗа, как мы помним статистически незначимо. Так мы видим, что кто-то бы к нам не пришел, в первый месяц работы его/её средняя эффективность будет равна 80%, а к 10 месяцу работы увеличится до ~110%.

Заключение

Читатель может возразить, что дилемму можно было бы решить, добавив на дашборд HR-метрику стажа. Сделать это было бы не лишним, правда есть ряд «но», почему этого недостаточно.

Во-первых, когда у нас всего две переменные, человеческий мозг еще способен уловить закономерности из графиков с несколькими категориями и данными за год. Но когда добавляется еще пять или десять переменных, понять фундаментальную зависимость становится невозможно.
Во-вторых, мы перекладываем ответственность за анализ на пользователя, у которого может быть совсем другой опыт, и разные пользователи могут прийти к разным выводам, на основе одних и тех же графиков и HR-метрик.
В-третьих, мы не улавливаем фундаментальные связи между факторами и не получаем коэффициентов, необходимых для прогнозирования и не строим сам прогноз.

Означает ли это, что HR-метрики бесполезны? Нет, не означает. Мы используем термометр для измерения температуры тела и тонометр для давления, но знаем, что этого недостаточно для полноценных и комплексных выводов о состоянии здоровья.

Пытливый читатель может указать на то, что если бы мы построили простую линейную регрессию только с переменной ВУЗ, мы бы пришли к неверному выводу, так как фактор оказался бы статистически значимым. В ответ на это я отмечу, что регрессии и другие модели, в отличие от простого расчета метрик, требуют осмысленного, глубокого и планируемого подхода к анализу. В четвертой статье серии мы рассмотрим фреймворк DAG, который помогает избежать ловушки выбора неправильных переменных или невключения правильных в модель регрессии.

Ссылки