Байесовский фреймворк и HR задачи. Часть 3 - MCMC

Томас Байес сформулировал байесовский вывод ещё в 18 веке, но массовую популярность он стал приобретать только в последние десятилетия. В то время как фреймворк проверки нулевой гипотезы был предложен Роналдом Фишером в 20 веке спустя 200 лет после Байеса.

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

В прошлых статьях мы рассматривали аналитический способ расчёта [1] и сеточное приближение для получения байесовского вывода [2]. Эти способы хороши для простых моделей, но для сложных моделей, те в которых надо оценить сразу несколько параметров, они не подходят. Для таких моделей используют методы, требующие вычислительные способности современных компьютеров – это марковские цепи Монте-Карло (MCMC). Природа MCMC, на мой взгляд, достаточно сложна, нас же интересует суть – с помощью этого метода мы численно получаем большое количество выборок и можем аппроксимировать интересующее распределение, даже если оно очень сложное и его трудно рассчитать напрямую [3].

Для получения MCMC существуют специальные инструменты – это JAGS и STAN. Оба инструменты доступны в R. Их сравнения и описание выходят за пределы статьи. Сегодня мы будем практиковаться с JAGS – для этого его нужно скачать и установить на ваш компьютер [4], а в R поставить библиотеку rjags.

На примере простой задачи мы попробуем залезть под капот JAGS и будем писать модели самостоятельно, но в практике, я рекомендую, обращаться к высокоуровневым скриптам Джона Крушке, которые идут дополнением к его книге "Doing Bayesian Data Analysis" и сэкономят много вашего времени [5].

Представим, что мы ищем решение достаточно распространенной HR задачи – повысятся ли продажи сотрудников, после прохождения тренинга или нет?

Начнем с наших данных. Пусть у нас есть 25 сотрудников, которых мы обучили и, по которым, нам известен средний уровень продаж до и после обучения, выраженный в рублях. Мы оценили дельту в этих средних и хотим найти ответ на вопрос - какова разница в среднем и стандартном отклонении продаж до и после обучения?

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

Мы переходим в R, где и решаем эту задачу. Прежде всего сгенерируем наши данные на основе нормального распределения [6].

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

#Зерно для воспроизведения результата
set.seed(1800)

#Данные по продажам до обучения: среднее = 5000, 
#стандартное отклонение = 200
before <- rnorm(25, 5000, 200)

#Данные по продажам после обучения: среднее = 6000, 
#стандартное отклонение = 800
after <- rnorm(25, 5500, 900)

#Создаем датасет с данными и дельтой
data <- as.data.frame(cbind(before, after, diff = after - before))

Посмотрим на распределение нашей дельты, которую мы передали в переменную diff.

Медиана для средней дельты = 205 рублям, для стандартного отклонения 1054. Чуть менее половины наблюдений меньше нуля, что наталкивает на мысль, что в обучении смысла не было. Но мы продолжим наш анализ.

Стоит сразу отметить, что синтаксис rjags довольно сильно отличается от привычного R. При создании модели MCMC с помощью JAGS мы должны выполнить 4 шага:

Определение модели.
Компиляция модели.
Симуляция апостериорного распределения.
Проверка качества и интерпретация результата.

Модель задается в виде текста. В "блоке" likelihood мы определяем результат для каждого нашего сотрудника с помощью функции плотности вероятности нормального распределения. Обратите внимание, что в JAGS используется не привычное нам стандартное отклонение, а точность (precision), по сути, это обратная величина, её можно задать как 1 / s^2 или s^(-2).

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

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

#1. Определяем нашу модель
model <- "model{

    # Likelihood
    for(i in 1:length(Y)) {
        Y[i] ~ dnorm(m, s^(-2))
    }


    # Prior for m and s
    m ~ dnorm(1000, 500^(-2))
    s ~ dunif(0, 2000)
    
   
}"

После текстового определения модели переходим к её компиляции

#2. Компилируем модель
model_jags <- jags.model(

  #Устанавливаем связь с моделью
  textConnection(model), 
  
  #Передаем нашу дельту в качестве данных
  data = list(Y = data$diff), 
  
  #Настраиваем генератор случайных чисел
  inits = list(.RNG.name = "base::Wichmann-Hill", .RNG.seed = 1989))

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

#3. Симулируем апостериорное распределение c помощью MCMC
model_sim <- coda.samples(

  #Указываем скомпилированную модель
  model = model_jags,
  
  #Выбираем переменные, которые оцениваем
  variable.names = c("m", "s"), 
  
  #Кол-во итераций
  n.iter = 10000)

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

# Передадим результат модели в таблицу
result <- as.data.frame(
  cbind(model_sim[[1]], i = 1:length(model_sim[[1]]))
                       )

# Построим графики трассировки      
# Для среднего                                        
ggplot(result, aes(x = i, y = m)) + 
  geom_line(color='#096CA9')
  
# Для стандартного отклонения
ggplot(result, aes(x = i, y = s)) + 
  geom_line(color='#096CA9')

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

  
  #Построим графики для оцениваемых параметров m и s
  
  ggplot(result, aes(x = m)) + 
  geom_histogram(bins=28, color = '#7604A8', fill = '#096CA9') +
  scale_x_continuous(name = "m - среднее") +
  ylab("Количество")
  
  ggplot(result, aes(x = s)) + 
  geom_histogram(bins=28, color = '#7604A8', fill = '#096CA9') +
  scale_x_continuous(name = "s - стандартное отклонение") +
  ylab("Количество")

Медиана для среднего составляет 413 рублей, для стандартного отклонение 1094 рублей. Что является компромиссом между априорными данными и наблюдаемыми на 25 сотрудниках.

95% наблюдений для среднего (m) лежат от 20 до 827 рублей - довольно большой разбег. Я использовал обычные процентили для интервала, т.к. данные нормально распределены, но обычно, в байесовском фреймворке для этого высчитывают HDI - интервал наибольшей плотности, который эффективен на всех типах распределения. Можно воспользоваться функцией HPDPI() из пакет rethinking.

Мы видим на графике, что наше распределение для среднего перекрывает 0, оценим долю случаев ниже нуля - это всего лишь 2%.

(sum(result$m < 0) / length(result$m)) * 100

А что будет, если мы для сравнения подходов, проверим гипотезу с помощью фреймворка тестирования нулевой гипотезы? Давайте возьмем Т-критерий Стьюдента.

t.test(after, before, paired=T)

Получаем p-value = 0.1645, что однозначно говорит о том, что мы не можем отклонить нулевую гипотезу.

Но результат, полученный с помощью байесовского вывода намного более содержателен:

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

Ссылки:

Версия статьи на английском.