Регрессиада. Часть 3. Мир полон неопределенности

В двух прошлых статьях [1][2], опираясь на классическую статистику, мы рассмотрели простую и множественную линейные регрессии. В этой и следующей статье мы рассмотрим не только другие виды регрессий, но и перейдем от классической статистики к байесовской. Байесовская статистика уже была предметом обсуждения в серии моих статей: [3], [4], [5].

Главные преимущества байесовского подхода:

Мы можем учесть ранее полученные знания или убеждения при построении модели.
Результаты модели представляют собой вероятностное распределение, а не бинарную оценку «значимо/незначимо», больше никакого p-value!
Позволяет делать причинно-следственные выводы, на основе результатов моделей, о чем мы будем говорить в следующей статье.
Способен обнаружить эффект на небольших по размеру выборках.

Работать с байесовской статистикой мы будем во фреймворке Ричарда МакЭлрита, подробно описанном в его книге [6]. Если у вас есть вопросы о том, как настроить среду на своем компьютере для работы в этом фреймворке, то воспользуйтесь этим руководством [7].

Введение

Компания «Линейные уравнения» открыла новое подразделение в одном из регионов страны. Это не первый опыт расширения деятельности компании, так как подобные мероприятия происходят время от времени. Перед отделом HR стоит задача обеспечить уровень найма в новом регионе не менее 50 сотрудников каждую неделю. Такой показатель найма является типичным для компании. По опыту известно, что вклад одного рекрутера в количество приёмов варьируется от 50 до 75%. Однако, поскольку в новом регионе компания «Линейные уравнения» проводит подбор персонала только первый месяц, на данный момент доступны данные лишь за четыре недели. С одной стороны, эти данные не подтверждают прогнозы, с другой стороны, известно, что не все необходимые инструменты, включая ATS и процессы, были внедрены в полной мере в новом подразделении. Аналитиков просят подготовить расчет, который покажет вероятность выполнения планов найма в размере 50 сотрудников еженедельно.

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

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

Для разных вариантов зависимых переменных и распределений существуют не просто разные виды регрессий, но разные семейства таких моделей. Одно из таких семейств – это обобщенные линейные модели (Generalized Linear Models, GLM) [8]. В отличие от моделей, которые мы рассматривали ранее, ключевой особенность GLM является функция связи, которая гарантирует, что выходные данные модели соответствуют типу и распределению зависимой переменной. Я не буду углубляться в детали, но всех, кого интересуют подробности я отсылаю к книге Питера Дана про GLM модели [9].

На моей схеме, которую вы уже видели в первой статье серии [1], я отметил GLM модели зеленым цветом. Сегодня мы рассмотрим случай, когда зависимая переменная счётная, а в следующей статье — когда она бинарная.

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

Формула регрессии Пуассона

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

Теперь мы не можем утверждать, что изменение переменной x на одну единицу приводит к изменению y на b единиц. Вместо этого, мы говорим, что при изменении переменной x на одну единицу переменная y меняется на b процентов, при условии, что другие переменные в модели неизменны. Мы ещё рассмотрим интерпретацию в практической части статьи.

Кроме того, переменная y должна следовать распределению Пуассона [10]. Это распределение моделируется всего лишь одним известным параметром лямбда – который отражает среднее значение. Из этого вытекает, что среднее и дисперсия определяется одним параметром лямбда, который означает, что они равны. В реальности дисперсия чаще больше среднего, что и называется сверхдисперсией.

Формула распределения Пуассона

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

Пропорциональность апостериорного распределения и произведения функции правдоподобия на априорное распределение

p(B|Y) – апостериорное распределение бета-коэффициентов;

L(B) – функция правдоподобия бета-коэффициентов (оцениваемая на основе наших данных);

p(B) – априорное распределение бета-коэффициентов.

Практика

Загрузим имеющиеся данные компании «Линейные уравнения» и взглянем на них [11]. Датасет содержит всего 4 наблюдения, соответствующих недельным отчётам. Переменные включают в себя следующую информацию:

recruiters – количество рекрутеров;

hires – количество приёмов.

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

Классическая модель

Прежде чем бросаться в объятия байесовской регрессии давайте сделаем плавный переход и построим регрессию Пуассона в рамках классического фреймворка, который вам уже хорошо знаком и понятен по двум предыдущим статьям серии. Делается это почти точно также, только вместо функции lm() мы используем функцию glm(), в которой указываем функцию связи, для случая регрессии Пуассона это poisson(link = "log").

# Строим регрессию Пуассона
classic_glm <- glm(hires ~ recruiters, poisson(link = "log"), data = data)


# Результаты модели
summary(classic_glm)

Call:
glm(formula = hires ~ recruiters, family = poisson(link = "log"), 
    data = data)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  2.97318    1.10704   2.686  0.00724 **
recruiters   0.05972    0.12222   0.489  0.62512   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 3.2577  on 3  degrees of freedom
Residual deviance: 3.0188  on 2  degrees of freedom
AIC: 28.387

Number of Fisher Scoring iterations: 4

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

# Экспонирование коэффициентов модели
exp(coef(classic_glm))

(Intercept)  recruiters 
  19.554019    1.06153

Интерсепт мы могли бы интерпретировать так: когда количество рекрутеров равно 0, то среднее количество приёмов составляет ~20 человек. С бета-коэффициентом при переменной recruiters, история другая. Здесь мы говорим, что при увеличении количества рекрутёров на 1 единицу, количество приёмов меняется на +6%, так как: (1.06 - 1) * 100% = 6%.

Согласитесь, что результат более, чем странный, так как вряд ли вклад рекрутёра в приёмы это всего лишь 6%, при этом коэффициент статистически незначим, что намекает на то, что рекрутеры вообще никакой роли не играют в процессе подбора! Очень спорное утверждение, не думаю, что аналитики согласятся его презентовать. Но ситуация объяснима, у нас очень мало данных, а также мы знаем, что ещё не все необходимые для разгона найма условия, воссозданы на должном уровне в новом подразделении.

Байесовские модели

Любая байесовская модель складывается из трех составляющих:

То, что нам известно заранее – априорное знание;
Наблюдаемые данные – правдоподобие;
Результат модели – апостериорное знание или компромисс между априорным знанием и правдоподобием.

Я нахожу эту идею очень красивой, и именно поэтому я влюблен в байесовскую статистику!

Существует множество различных способов и сред, в которых можно построить байесовскую регрессию, но мы будем делать это в R, и не просто в R, а именно во фреймворке Ричарда МакЭлрита.

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

Построим первую модель и на её примере подробно разберём, как это делать.

# Подключаем библиотеку rethinking
library(rethinking)

# Модель с неопределенным априорным представлением
wide_prior <- ulam(
  alist(
    hires ~ dpois(mu),
    log(mu) <- a + b_recruiters * recruiters,
    
    a ~ dnorm(0, 3),
    b_recruiters ~ dnorm(0, 3)
  ),
  data = data,
  chains = 4, cores = 4, 
  iter = 5000, warmup = 1000,
  messages = FALSE, refresh = 0
)

Давайте разбираться с содержимым:

ulam() – функция, используемая для построения модели. Внутри неё мы определяем как расчётные формулы, так и прочие параметры модели.
alist() – здесь располагается ядро модели.

Я буду читать содержимое alist() не по порядку, и далее вы поймете, почему так проще и нагляднее. Помните, я говорил, что внутри моделей мы задаем наши априорные знания о коэффициентах, а также сказал, что в первой модели мы ничего о них не знаем? Как же это выразить математически?

Представим, что интерсепт и бета-коэффициент нормально распределены со средним 0 и гигантским стандартным отклонением равным 3, что выражает высшую степень неопределенности. Тогда внутри alist() мы напишем:

a ~ dnorm(0, 3) – интерсепт.
b_recruiters ~ dnorm(0, 3) – бета-коэффициент для переменной recruiters.

Названия коэффициентов могут быть любыми.

Определив априорное представление нам нужно задать функцию правдоподобия. Выразим среднее количество приёмов mu, как линейную зависимость от интерсепта и вклада рекрутеров:

mu <- a + b_recruiters * recruiters.

Но поскольку далее мы будем моделировать распределение Пуассона, это среднее не может быть отрицательным. Решается это просто - обернем среднее функцией натурального логарифма:

log(mu) <- a + b_recruiters * recruiters.

Наша зависимая переменная счётная, и мы решили использовать пуассоновское распределение, чтобы смоделировать количество приёмов, где в качестве лямбды мы берем логарифм от mu.

hires ~ dpois(mu).

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

За пределами alist() мы указываем остальные технические параметры:

data = data – ссылаемся на наш датафрейм.
chains = 4 – количество цепочек MCMC [12]. Чем больше, тем надежнее результат и тем дороже его расчёт.
cores = 4 – количество ядер вашего процесса, которое будет использовано при расчете, чем больше – тем быстрее.
iter = 5000 – количество итераций, которые будут затрачены на модель.
warmup = 1000 – сколько итераций модель будет «разогреваться».
messages = FALSE, refresh = 0 – позволяют подавить излишние сообщения на этапе расчёта модели.

Теперь посмотрим на результаты модели

# Анализируем результат модели
precis(wide_prior)

             mean   sd  5.5% 94.5% n_eff Rhat4
a            2.58 1.03  0.97  4.20  2133     1
b_recruiters 0.10 0.11 -0.08  0.28  2148     1


# Экспонируем результаты
exp(coef(wide_prior))

           a b_recruiters 
   13.139963     1.108035

Мы начинаем прочтение сводки с Rhat4. Если этот показатель равен 1, значит с моделью все в порядке; если нет, то наша модель не сошлась. Одним из первых решений может быть попытка увеличения количества итераций. В нашем случае все в порядке.

Коэффициент b_recruiters низкий, всего лишь 10% после экспонирования, но всё же выше 6% из классической модели. При этом, если мы взглянем на колонки 5.5% и 94.5%, то увидим 89% интервал распределения коэффициента, который покрывает 0, что в байесовском подходе можно интерпретировать как незначимость.

А теперь настала пора раскрыть всю мощь метода – задать определенные априорные представления. Мы знаем из условия, что вклад одного рекрутера в приёмы варьируется от 50 до 75%. Для выражения этого бета-коэффициента мы также будем использовать нормальное распределение. Среднее между 50 и 75% равняется 62.5%. Поскольку мы хотим выразить процент, то это будет натуральный логарифм от 1.625. Стандартное отклонение, где три 3 сигмы (99%), должны покрывать диапазон от 50 до 75% примерно равняется 0.03.

# Модель с конкретным представление о бета-коэффициенте
narrow_prior <- ulam(
  alist(
    hires ~ dpois(mu),
    log(mu) <- a + b_recruiters * recruiters,
    
    a ~ dnorm(0, 3),
    b_recruiters ~ dnorm(log(1.625),0.03)
  ),
  data = data,
  chains = 4, cores = 4, 
  iter = 5000, warmup = 1000,
  messages = FALSE, refresh = 0
)

Прежде всего полезно оценить, насколько наш код действительно отражает априорное представление.

# Получаем выборки из априорного распределения
prior <- extract.prior(narrow_prior)

# Визуализируем гистограмму
library(ggplot2)

prior |> as.data.frame() |> 
  ggplot(aes(x = exp(b_recruiters))) +
  geom_histogram(fill = 'gray') +
  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()) + geom_vline(xintercept  = 1.625, colour = 'red', size = 2)

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

# Анализируем результат модели
precis(narrow_prior)

              mean   sd  5.5% 94.5% n_eff Rhat4
a            -0.70 0.29 -1.15 -0.24  2782     1
b_recruiters  0.46 0.03  0.41  0.51  2873     1


# Экспонируем результаты

exp(coef(narrow_prior))
           a b_recruiters 
   0.4978934    1.5863304

Наша модель нашла компромисс между априорным знанием и реальными данными, и вклад одного рекрутера оказался равен 58%. Интерсепт в данном случае не имеет логической интерпретации, так как утверждает, что при количестве рекрутеров равном нулю, количество приёмов снижается на 50%.

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

result <- as.data.frame(precis(narrow_prior))

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

Как видим, все наши коэффициенты значимы.

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

# Диагностика модели
traceplot(narrow_prior, window = c(1000, 5000))

Четыре цепочки MCMC отмечены разным цветом, перекрывают друг друга, какой-либо паттерн не наблюдается

Считаем прогноз

Мы построили модель, но изначальный вопрос, обращенный к нам, звучал так: какова вероятность, что подразделение сможет нанимать 50 и более сотрудников еженедельно. И вновь, байесовская модель позволяет нам ответить на этот вопрос. Сделаем прогноз на основе модели, при условии, что в штате работает 10 рекрутеров, что соответствует последней доступной дате в наших данных.

# Строим прогноз на основе модели
pr <- link(narrow_prior, data.frame(recruiters = 10))

# Созадим датасет
result = data.frame(
  quantity = round(pr),
  group = ifelse(round(pr) > 50, "Больше 50 приёмов", "Меньше 50 приёмов")
)

# Посчитаем проценты
percentage_above_50 <- length(pr[pr > 50]) / length(pr) * 100
percentage_label <- sprintf("Вероятность приёмов\n>=50: %.1f%%", percentage_above_50)


# Визуализируем
ggplot(result, aes(x = quantity, fill = group)) +
  geom_histogram(binwidth  = 1) +
  scale_fill_manual(values = c("Больше 50 приёмов" = "red", "Меньше 50 приёмов" = "gray")) +
  labs(x = "Кол-во приёмов",
       y = "Частота",
       fill = "Группа") +
  annotate("text", x = Inf, y = Inf, label = percentage_label, 
           hjust = 1.1, vjust = 1.5, size = 5, color = "black") +
  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())

Мы получили вероятностный ответ на основе модели, который с учётом имеющегося априорного знания, доступных данных за прошедшие 4 недели и основываясь на том, что в штате работает 10 рекрутеров, показывает, что вероятность приёма более 50 человек еженедельно составляет 53%.

При этом мы можем посчитать второй сценарий: что будет, если в штате будет 11 рекрутеров? В этом случае мы получим 100% вероятность выполнения плана. Это может показаться не сразу очевидным, если вы привыкли работать с линейными моделями, но здесь стоит ещё раз подчеркнуть, что природа наших коэффициентов основана не на сложении, а на умножении, то есть каждый рекрутер даёт рост числа приёмов на значительные 58%.

Заключение

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

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

Ссылки