贝叶斯比较
很多人在离开学校时都会思考未来。我们从小就被问到长大后想要成为什么，然后花了13年在前期教育中。在公共政策领域，政府、天主教和独立学校体系之间的差异被广泛讨论，特别是在非政府学校的公共资助以及资源分配方面。
不同学校部门的学生结果是否存在实质差异？
在澳大利亚的维多利亚州，州政府每年都会进行一项调查来了解高中毕业后的结果（On Track Survey）。这个数据集（根据 CC BY 4.0 许可证提供）跨足了多年，我们将仅关注此文章撰写时最新的数据，即 2021 年的数据。
本文将尝试回答这些问题，应用贝叶斯分析方法，使用 R 包 brms。
加载库和数据集
以下是加载所需库和数据集的代码。数据集以宽格式呈现，有许多合并的单元格。R 不太适应这种格式，因此我们需要做一些硬编码来重新标记向量并创建一个整洁的数据框。
library(tidyverse) #Tidyverse meta package
library(brms) #Bayesian Modelling Package
library(tidybayes) #Tidy Helper Functions and Visualisation Geoms for Bayesian Models
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 代码 — 每个代码的管理 ID
学校名称
部门 — G = 政府，C = 天主教，I = 独立
地区或郊区
完成 12 年级的总学生人数
同意追踪
应答者
其他列代表每种结果的应答者百分比
在我们的横断面研究中，我们对部门和结果的比例感兴趣，因此我们需要将这个宽数据框转换为长数据格式。
df_long <- df |>
mutate(across(5:14, ~ as.numeric(.x)), #convert all numeric fields captured as characters
across(8:14, ~ .x/100 * respondants), #calculate counts from proportions
across(8:14, ~ round(.x, 0)), #round to whole integers
respondants = bachelors + deferred + tafe + apprentice_trainee + employed + looking_for_work + other) |> #recalculate total respondents |>
filter(sector != 'A') |> #Low volume
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 = 'Sector', y = 'Count', title = 'Count of Schools by Sector', fill = 'Sector') +
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 = 'Sector', y = 'Proportion', fill = 'Outcome', title = 'Boxplot of Respondant Proportions by Sector and Outcome') +
scale_fill_viridis_d(begin = 0.2, end = 0.8) +
theme_ggdist()
这些分布呈现了一个有趣的故事。学士学位的结果在所有部门中的变异性最大，独立学校拥有最高的学生比例中位数，选择攻读本科教育。值得注意的是，政府学校在高中毕业后就业的学生中，具有最高的中位数比例。其他所有结果变化不大，我们将很快回顾这些结果。
统计建模 — 贝叶斯似然回归
我们关注学校毕业生的比例和他们在高中毕业后的结果。在这种情况下，贝叶斯似然是我们最好的选择。brms 是 Buerkner 开发的出色的包，用于构建贝叶斯模型。我们的目标是模拟结果和部门的比例分布。
贝叶斯回归模型包括两个参数，μ 和 φ。μ代表均值比例，φ代表精度（或方差）的度量。
我们的第一个模型如下所示，请注意，我们已经从部门和结果之间的交互开始，因为这是我们希望模型回答的问题，因此这是一个ANOVA模型。
我们要求模型为每个部门和结果的组合创建单独的 Beta 项，带有一个汇总的 φ 项，或者用相同的方差模拟不同比例的均值。我们期望这些比例中的 50% 位于 logit 比例尺上的 (-3, 1) 或比例上的 (0.01, 0.73)。这是足够宽泛但有信息的先验知识。全局 Phi 估计是一个正数，所以我们使用了足够宽泛的伽马分布。
# Modelling Proportion with Sector:Outcome Interaction term using Beta - Pooled Phi term
m3a <- brm(
proportion ~ sector:outcome + 0,
family = Beta,
data = df_long |> mutate(proportion = proportion + 0.01), # Beta regression can't take outcomes that are 0 so we fudge here by adding tiny increment
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% 的增量 — 这是因为 Beta 回归不能处理零值。我们尝试过使用零膨胀贝叶斯回归来模拟这一点，但需要更长时间来收敛。
同样，我们可以在没有汇总的 φ 的情况下建模，从上面观察到的分布情况来看，这在直觉上是有道理的，每个部门和结果的组合之间存在差异。模型定义如下。
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 = 'Distribution of Observed and Replicated Proportions by Sector and Outcome')
模型m3b，鉴于其非汇总的方差或phi项，能够更好地捕捉各部门和结果比例分布的变化。
贝叶斯风格的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 = 'Proportional Difference', y = 'Sector', title = 'Differences in Posterior Means by Sector and Outcome', fill = 'Sector')
或者，我们可以总结所有这些分布，制作一个整洁的表格，以便更容易解释，并包括一个 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 = 'Sector Marginal Distributions by Outcome') |>
#tab_stubhead(label = 'Sector Comparison') |>
fmt_percent() |>
gtExtras::gt_theme_pff()
例如，如果我们要在正式报告中总结一个比较，我们可以写下以下内容。
政府学校的学生开始本科学位的可能性较低，而他们的同行在 VCE 完成后的天主教学校和独立学校中更高。
平均而言，政府学校学生中有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%。
总结与结论
在本文中，我们演示了如何使用贝叶斯建模来模拟比例数据，采用Beta似然函数，然后执行贝叶斯ANOVA以探索各个部门之间比例结果的差异。
我们并没有试图创建这些差异的因果理解。可以想象，有许多因素会影响个别学生的表现，如社会经济背景、父母的教育程度，以及学校层面的影响、资源等。这是一个非常复杂的公共政策领域，往往会陷入零和的争论中。
✄-----------------------------------------------看到这里，说明你喜欢这篇文章，请点击「在看」或顺手「转发」「点赞」。
欢迎微信搜索「panchuangxx」，添加小编磐小小仙微信，每日朋友圈更新一篇高质量推文（无广告），为您提供更多精彩内容。
▼ ▼ 扫描二维码添加小编 ▼ ▼