Join
H0H1: про HR-аналитику
@h0h1_hr_analytics
HR аналитика
May 25

Регрессиада. Часть 4. За пределами корреляции

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

DAG (directed acyclic graph) – это направленный ациклический граф, центральный инструмент в причинно-следственном анализе, предложенный и разработанный Иудеем Перлом [1]. К первоисточнику отсылаю самых смелых, так как там много сложной математики и философии. Мы же будем полагаться на переработку этой темы Ричардом МакЭлритом [2].

Введение

Руководство компании «Линейные уравнения» хочет понять, является ли участие в проекте кадрового резерва причиной перевода сотрудника на руководящую должность, когда такая возможность возникает. Известно, что для того, чтобы стать участником кадрового резерва, необходимо иметь высокие оценки за предыдущее ревью. Также известно, что переводы на руководящую должность могут происходить и без участия в кадровом резерве, но при условии наличия высоких оценок за ревью.

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

В нашем распоряжении 1000 наблюдений за длительный период. Данные выглядят следующим образом:

Глядя на сводку, допустимо предположить, что все сотрудники с высокими оценками могли бы получить повышение вне зависимости от участия в кадровом резерве, так как перевод может состояться при достижении одного из условия, а сам факт высокой оценки — это предпосылка включения в резерв. Как же изолировать каузальный эффект участия в кадровом резерве? Разобраться во всём этом нам помогут DAG и байесовская регрессия.

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

Прежде, чем переходить к DAG, давайте определимся с каким видом регрессионной модели мы столкнулись в этот раз. Факт перевода на руководящую должность – это бинарная переменная, которая может принимать только два значения: перевод состоялся или не состоялся; да или нет; 0 или 1 и т.д. Для того чтобы работать с таким типом зависимой переменной нам понадобится (биномиальная) логистическая регрессия.

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

Лично я нахожу интересным график логистической функции и её производной. Из которого видно:

  • все значения y ограничены на отрезке от 0 до 1;
  • функция нелинейна;
  • скорость изменения не равномерна.

Практика

Рисуем DAG

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

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

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

В нашем распоряжении следующие переменные:

  • rotation – факт перевода: перевод состоялся или не состоялся;
  • pool – факт участия в кадровом резерве: да или нет;
  • mark – оценка по результатам ревью от 1 до 5.

Построим DAG.

# Подключаем необходимые библиотеки
library(tidyverse)  
library(dagitty) 
library(ggdag)

# Определяем DAG на основе известных условий
dag <- ggdag::dagify(
  rotation ~ pool,
  rotation ~ mark,
  pool ~ mark,
  
  exposure = 'pool',
  outcome = 'rotation'
)

# Визуализируем
p <- ggdag::ggdag(dag, node_size = 20) + theme_dag()

p$layers[[3]]$mapping <- 
  aes(colour = case_when(
    name == "rotation" ~ "Зависимая",
    TRUE ~ "Независимая"
  ))

p + geom_dag_edges() + 
  scale_color_manual(values = c("Независимая" = "#288ba8", "Зависимая" = "#F0A70B")) +
  theme(legend.position = "bottom", 
        legend.title = element_text(size = 11))

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

Определившись с DAG, мы готовы приступить к построению модели логистической регрессии, но вначале загрузим данные и посмотрим на них.

Строим модель логистической регрессии

Загрузим данные с моего GitHub и взглянем на них.

# Загрузка данных
url <- "https://raw.githubusercontent.com/alexander-botvin/h0h1_about_hr_analytics/main/Regressions/example_three.csv"
data <- read.csv(url)

# Структура данных
str(data)

'data.frame':	1000 obs. of  4 variables:
 $ X       : int  1 2 3 4 5 6 7 8 9 10 ...
 $ mark    : int  4 5 3 5 1 4 2 2 1 3 ...
 $ pool    : chr  "Нет" "Нет" "Нет" "Да" ...
 $ rotation: chr  "Перевод состоялся" "Перевод не состоялся" "Перевод не состоялся" "Перевод состоялся" ...

На основе сводки видно, что данные требуют определённого преобразования:

  1. Переменные должны быть перекодированы как факторы. В том числе оценки mark мы тоже интерпретируем как уровни переменной, а не как целые числа.
  2. Для классического фреймворка было бы достаточно этого, но для байесовского фреймворка Ричарда МакЭлрита переменная, отражающая перевод на должность, должна быть строго представлена в виде 0 или 1.

В этот момент может возникнуть вопрос: если зависимая переменная числовая, то почему не использовать множественную линейную регрессию? Ответ прост, нормальное распределение, которое находится под капотом множественной линейной регрессии не имеет никаких ограничений на принимаемые значений, в то время как нас интересует результат равный 0 или 1.

Выполним соответствующие преобразования данных.

# Трансформируем переменные в факторы и вручную выставим уровни
data$pool <- factor(data$pool, c('Нет', 'Да'))
data$mark <- factor(data$mark, c('1', '2', '3', '4', '5'))
data$rotation <- factor(data$rotation, c('Перевод не состоялся', 'Перевод состоялся'))

# Зависимая переменная в виде 0 или 1
data$rotation_num <- as.numeric(data$rotation) - 1

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

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

Также стоит отметить, что в этот раз все наши переменные категориальные. Мы обсуждали такие переменные и как они могут быть обработаны на уровне моделей во второй статье серии [6]. Если мы работаем во фреймворке rethinking, то мы используем индексные переменные, так что для каждого индекса высчитывается свой бета-коэффициент. Например, наша переменная оценки имеет 5 уровней, значит наш коэффициент будет определяться так: b_mark[mark], где название коэффициента может быть произвольным, в то время как ссылка на индексную переменную должна совпадать с одноимённым полем в датафрейме.

Зададим модель:

library(rethinking)

model <- ulam(
  alist(
    rotation_num ~ dbinom(1, p),
    logit(p) <- a + b_pool[pool] + b_mark[mark],
    
    a ~ dnorm(0, 10),
    b_pool[pool] ~ dnorm(0, 1),
    b_mark[mark] ~ dnorm(0, 1)    
  ),
  data = data,
  chains = 4, cores = 4, 
  iter = 5000, warmup = 1000,
  messages = FALSE, refresh = 0
)

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

# Обратите внимание, что когда мы работаем с категориальными переменными,
# то мы должны использовать дополнительный аргумент depth=2
summary_model <- precis(model, depth=2)

# Экспонируем коэффициент
summary_df <- as.data.frame(summary_model)[,c('mean', '5.5%', '94.5%')]
summary_df[] <- lapply(summary_df, function(x) if (is.numeric(x)) exp(x) else x)

# Собираем оригинальные названия уровней переменных
mark_labels <- paste("b_mark[", 1:length(levels(data$mark)), "]: ", levels(data$mark), sep = "")
pool_labels <- paste("b_pool[", 1:length(levels(data$pool)), "]: ", levels(data$pool), sep = "")
labels <- c(mark_labels, pool_labels)

# Сопоставляем названия переменных
rownames(summary_df)[4:8] <- mark_labels
rownames(summary_df)[2:3] <- pool_labels

Мы готовы построить график.

library(ggplot2)

ggplot(summary_df, aes(y = rownames(summary_df), x = mean)) +
  geom_point() + 
  geom_errorbarh(aes(xmin = `5.5%`, xmax = `94.5%`), height = 0.5, size = 1) +
  labs(x = "Estimate", y = "Term") +
  scale_y_discrete(limits = rev(rownames(summary_df))) +
  theme_minimal() + 
  theme(legend.position = "bottom", panel.grid.major.x  = element_blank() , panel.grid.minor.x  = element_blank()) + geom_vline(xintercept  = 1, colour = "red")

Поскольку коэффициенты были экспонированы, то граница выставлена в значение 1, так как exp(0) = 1

Мы видим, что участие в кадровом резерве имеет бета-коэффициент около 3.2, а правая граница 89% интервала превышает 10! То есть, в среднем, участие в кадровом резерве повышает вероятность перевода на (3.2 – 1) * 100% = 220%! Также мы видим, что высокие оценки тоже играют свою роль: оценка равная 5 дает однозначный вклад в вероятность перевода, оценка равная 4 имеет некоторую степень неопределенности, так как левый хвост пересекает 1.

Итак, глядя на коэффициенты модели, готовы ли мы утверждать, что участие в кадровом резерве — это причина переводов? Нет, ещё не готовы.

Делаем причинно-следственный вывод

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

  • Фактическое состояние — это то, что мы имеем. К примеру, результат нашей модели и/или прогноз на её основе.
  • Контрфактическое состояние — что было бы, если. В нашем кейсе мы можем задаться вопросом: что было бы с переводами, если бы все люди, которые были в кадровом резерве, не приняли бы участия в нём?

Начнем с того, что посчитаем оба состояния.

# Создадим функцию, которая будет генерировать прогноз
simulation <- function(model, pool, mark){
  p <- link(model, data = data.frame(pool = pool, mark = mark))
  return(apply(p, 2, mean))
}

# Симулируем фактическое состояние
factual <- simulation(model, data$pool, data$mark)

# Симулируем контфрактическое состояние,
# будто бы в кадровом резерве никто не состоял
data$new_pool <- factor('Нет', levels = levels(data$pool))
counterfactual  <- simulation(model, data$new_pool, data$mark)

# Для удобства добавим обе вероятности в исходный датфрейм
data$factual <- factual
data$counterfactual <- counterfactual

Построим график распределения вероятности получить перевод для этих двух сценариев.

# Перейдем от широкого формала таблицы к длинном
library(reshape2)
long_data <- melt(data[, c('factual', 'counterfactual')])

ggplot(long_data, aes(value, fill = variable)) +
  geom_density(alpha = 0.5) +
  xlim(c(0, 1)) +

  scale_fill_manual(labels = c("Фактическое", "Контрфактическое"), values = c("red", "blue")) +
  labs(x = "Вероятность перевода",
       y = "Плотность",
       fill = "Состояние") +
  theme_minimal() + 
  theme(legend.position = "bottom", panel.grid.major.x  = element_blank() , panel.grid.minor.x  = element_blank(), panel.grid.major.y  = element_blank() , panel.grid.minor.y  = element_blank())

График говорит о том, что если бы все люди, которые на самом деле были участниками кадрового резерва, в какой-то альтернативной реальности ими не были, то вероятность получить повышение оказалась бы для них гораздо ниже (обычно решение об отнесении к классу определяется границей в 0.5). Так на основе нашего графика мы бы сказали, что не один из сотрудников не получил бы перевод в контрфактическом состоянии, не смотря на высокие оценки.

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

Для тех, кому хочется больше подробностей, на тему причинно-следственного вывода, то вновь и вновь отсылаю к МакЭлриту [2] и Гельману [7].

Заключение

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

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

Ссылки

  1. Judea Pearl. Causality Models Reasoning and Inference
  2. Richard McElreath. Statistical Rethinking
  3. Ботвин А.Ю. Регрессиада. Часть 3. Мир полон неопределенности
  4. Ботвин А.Б. Конфаундеры. Труба.
  5. http://www.machinelearning.ru/wiki/index.php?title=Биномиальное_распределение
  6. Ботвин А.Ю. Регрессиада. Часть 2. Метрики лгут
  7. Andrew Gelman. Regression and Other Stories
@alexander.botvin
May 25, 08:32
0 views
0 reactions