Сравнение результатов выпускников школы с использованием байесовского подхода с помощью R и brms

Сравнение результатов выпускников школы с использованием байесовского подхода в R с помощью brms

ANOVA — Байесовский стиль

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

Существует ли какая-либо реальная разница в результатах учащихся в зависимости от сектора школы?

В штате Виктория, Австралия, правительство штата проводит ежегодное исследование для определения результатов после окончания школы (On Track Survey). Этот набор данных (доступный по лицензии CC BY 4.0) доступен на протяжении нескольких лет, и мы рассмотрим только самый последний доступный на момент написания этой статьи, 2021 год.

В этой статье мы попытаемся ответить на эти вопросы, применяя методы байесовского анализа с использованием пакета R brms.

• Оптимизация поездки по Европе генетические алгоритмы и Google Maps API решают задачу коммивояжера
• Как преуспеть в курсе по науке о данных Полное руководство
• Создание инструмента отслеживания обновлений проекта в Excel с использованием VBA

Загрузка библиотек и набора данных

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

library(tidyverse) #Мета-пакет Tidyverselibrary(brms) #Пакет байесовского моделированияlibrary(tidybayes) #Функции помощи и графические элементы для байесовских моделейlibrary(readxl)df <- read_excel("DestinationData2022.xlsx",     sheet = "SCHOOL PUBLICATION TABLE 2022",     skip = 3)colnames(df) <- c('vcaa_code',                   'school_name',                   'sector',                   'locality',                   'total_completed_year12',                   'on_track_consenters',                   'respondants',                   'bachelors',                   'deferred',                   'tafe',                   'apprentice_trainee',                   'employed',                   'looking_for_work',                   'other') df <- drop_na(df)

Вот пример того, как выглядит наш набор данных.

В наборе данных 14 полей.

• VCAA Code — Административный идентификатор для каждого кода
• Название школы
• Сектор — Г = государственные, К = католические и Н = независимые
• Местоположение или пригород
• Всего учеников, закончивших 12-й класс
• Согласовавших участие в опросе
• Ответивших
• Остальные столбцы представляют проценты отклика для каждого результата

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

df_long <- df |>   mutate(across(5:14, ~ as.numeric(.x)), #преобразование всех числовых полей из символов в числа         across(8:14, ~ .x/100 * respondants), #расчет количества по пропорциям         across(8:14, ~ round(.x, 0)), #округление до целых чисел         respondants = bachelors + deferred + tafe + apprentice_trainee + employed + looking_for_work + other) |> #пересчет общего числа респондентов |>   filter(sector != 'A') |>  #Низкий объем   select(sector, school_name, 7:14) |>   pivot_longer(c(-1, -2, -3), names_to = 'outcome', values_to = 'proportion') |>   mutate(proportion = proportion/respondants)

Анализ исходных данных

Давайте кратко визуализируем и суммируем наши данные. В этом наборе данных указано 157 школ в государственном секторе, 57 школ в независимом секторе и 50 школ в католическом секторе.

df |>   mutate(sector = fct_infreq(sector)) |>   ggplot(aes(sector)) +    geom_bar(aes(fill = sector)) +     geom_label(aes(label = after_stat(count)), stat = 'count', vjust = 1.5) +    labs(x = 'Сектор', y = 'Количество', title = 'Количество школ по секторам', fill = 'Сектор') +    scale_fill_viridis_d(begin = 0.2, end = 0.8) +    theme_ggdist()

df_long |>   ggplot(aes(sector, proportion, fill = outcome)) +    geom_boxplot(alpha = 0.8) +    geom_point(position = position_dodge(width=0.75), alpha = 0.1, color = 'grey',  aes(group = outcome)) +    labs(x = 'Сектор', y = 'Пропорция', fill = 'Результат', title = 'Ящик с усами пропорций по секторам и результатам') +    scale_fill_viridis_d(begin = 0.2, end = 0.8) +    theme_ggdist()

Распределения рассказывают интересную историю. Результат “Бакалавр” показывает наибольшую изменчивость во всех секторах, при этом в независимых школах наибольшая медианная пропорция студентов выбирает образование на уровне бакалавра. Интересно отметить, что в государственных школах наибольшая медианная пропорция студентов переходит к трудоустройству после окончания школы. Все остальные результаты не слишком сильно отличаются – мы скоро вернемся к этому.

Статистическое моделирование – регрессия с бета-вероятностью

Мы сосредоточены на пропорциях студентов в школах и их результатах после окончания школы. В таких случаях лучше всего использовать бета-вероятность. brms – великолепный пакет от Бюркнера для разработки байесовских моделей. Нашей целью является моделирование распределения пропорций по результатам и секторам.

Бета-регрессия моделирует два параметра, μ и φ. μ представляет среднюю пропорцию, а φ является мерой точности (или дисперсии).

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

Мы просим модель создать отдельные бета-термы для каждой комбинации сектора и результата, с объединенным φ-термом, то есть моделировать разные средние пропорции с одинаковой дисперсией. Мы ожидаем, что 50% этих пропорций будут находиться между (-3, 1) на логит-шкале или (0.01, 0.73) в виде пропорций. Это достаточно широкий, но основанный на информации априорный. Глобальная оценка Phi является положительным числом, поэтому мы используем гамма-распределение, которое достаточно широко.

# Моделирование пропорции с взаимодействием сектора и результата с использованием бета-объединенного терма Phi m3a <- brm(  proportion ~ sector:outcome + 0,   family = Beta,  data = df_long |> mutate(proportion = proportion + 0.01), # Бета-регрессия не может принимать результаты, равные 0, поэтому мы здесь добавляем небольшое приращение  prior = c(prior(normal(-1, 2), class = 'b'),            prior(gamma(6, 1), class = 'phi')),  iter = 8000, warmup = 1000, cores = 4, chains = 4, seed = 246, save_pars = save_pars(all = T),  control = list(adapt_delta = 0.99, max_treedepth = 15)) |>   add_criterion(c('waic', 'loo'), moment_match = T)summary(m3a)

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

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

m3b <- brm(  bf(proportion ~ sector:outcome + 0,     phi ~ sector:outcome + 0),  family = Beta,  data = df_long |> mutate(proportion = proportion + 0.01),  prior = c(prior(normal(-1, 2), class = 'b')),  iter = 8000, warmup = 1000, cores = 4, chains = 4, seed = 246, save_pars = save_pars(all = T),  control = list(adapt_delta = 0.99)) |>   add_criterion(c('waic', 'loo'), moment_match = T)summary(m3b)

Диагностика и сравнение моделей

У нас есть две модели, и теперь мы сравним их точность прогнозирования на новых данных с помощью байесовской оценки LOO ожидаемой логарифмической плотности по точкам (elpd_loo).

comp <- loo_compare(m3a, m3b)print(comp, simplify = F)

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

alt_df <- df_long |>   select(sector, outcome, proportion) |>   rename(value = proportion) |>   mutate(y = 'y',          draw = 99) |>   select(sector, outcome, draw, value, y)sim_df <- expand_grid(sector = c('C', 'I', 'G'),            outcome = unique(df_long$outcome)) |>   add_predicted_draws(m3b, ndraws = 1200) |>   rename(value = .prediction) |>   ungroup() |>   mutate(y = 'y_rep',         draw = rep(seq(from = 1, to = 50, length.out = 50), times = 504)) |>   select(sector, outcome, draw, value, y) |>   bind_rows(alt_df)sim_df |>   ggplot(aes(value, group = draw)) +    geom_density(aes(color = y)) +    facet_grid(outcome ~ sector, scales = 'free_y') +    scale_color_manual(name = '',                        values = c('y' = "darkblue",                                  'y_rep' = "grey")) +    theme_ggdist() +    labs(y = 'Density', x = 'y', title = 'Распределение наблюдаемых и воспроизведенных пропорций по сектору и результату')

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

ANOVA — байесовский стиль

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

Ниже мы генерируем набор ожидаемых значений для каждой комбинации сектора и результата взаимодействия, а затем используем отличную функцию tidybayes::compare_levels(), которая выполняет тяжелую работу. Она вычисляет разницу в постериорных средних между секторами для каждого результата. Мы исключили результат “другое” для краткости.

new_df <- expand_grid(sector = c('I', 'G', 'C'),           outcome = c('apprentice_trainee', 'bachelors', 'deferred', 'employed', 'looking_for_work', 'tafe'))epred_draws(m3b, newdata = new_df) |>   compare_levels(.epred, by = sector, comparison = rlang::exprs(C - G, I - G, C - I)) |>     mutate(sector = fct_inorder(sector),         sector = fct_shift(sector, -1),          sector = fct_rev(sector))  |>   ggplot(aes(x = .epred, y = sector, fill = sector)) +      stat_halfeye() +      geom_vline(xintercept = 0, lty = 2) +       facet_wrap(~ outcome, scales = 'free_x') +      theme_ggdist() +      theme(legend.position = 'bottom') +      scale_fill_viridis_d(begin = 0.4, end = 0.8) +      labs(x = 'Пропорциональная разница', y = 'Сектор', title = 'Различия в постериорных средних по секторам и результатам', fill = 'Сектор')

В качестве альтернативы мы можем суммировать все эти распределения в удобной таблице для более простого понимания и 89% доверительным интервалом.

marg_gt <- epred_draws(m3b, newdata = new_df) |>   compare_levels(.epred, by = sector, comparison = rlang::exprs(C - G, I - G, C - I)) |>   median_qi(.width = .89) |>   mutate(across(where(is.numeric), ~round(.x, 3))) |>   select(-c(.point, .interval, .width)) |>   arrange(outcome, sector) |>   rename(diff_in_means = .epred,          Q5.5 = .lower,          Q94.5 = .upper) |>   group_by(outcome) |>   gt() |>   tab_header(title = 'Секторальные маргинальные распределения по результатам') |>   #tab_stubhead(label = 'Сравнение секторов') |>   fmt_percent() |>   gtExtras::gt_theme_pff()

Например, если мы хотим суммировать сравнение в официальном отчете, мы можем написать следующее.

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

В среднем, 42,5% (от 39,5% до 45,6%) учеников государственных школ, 53,2% (от 47,7% до 58,4%) учеников католических школ и 65% (от 60,1% до 69,7%) учеников частных школ начинают обучение в бакалавриате после окончания 12 класса.

Существует 89% вероятность того, что разница в поступлении на бакалавриат между католическими и государственными студентами составляет от 5,6% до 15,7% средним значением 10,7%. Кроме того, разница в поступлении на бакалавриат между частными и государственными студентами составляет от 17,8% до 27% средним значением 22,5%.

Эти различия существенны, и существует 100% вероятность, что эти различия не равны нулю.

Резюме и заключение

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

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

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