lvzhaopeng的数字花园

R语言tidymodels那些事 （十）tidymodels 与贝叶斯模型比较（深度解读）

10分钟阅读

缘由

艺多不压身

碎碎念

我们继续飞行在云层和雾气之间，飞行在一个无光无影的空洞世界中。 ——《风沙星辰》

Ref

11 Comparing Models with Resampling | Tidy Modeling with R

解读

本文用 Ames 数据集，演示如何用 tidymodels 做预处理/交叉验证/模型比较，并用 tidyposterior + rstanarm 做贝叶斯模型比较，从而得到更稳健且对决策有实际意义的结论。

library(tidymodels)
library(modeldata)
library(tidyposterior)
library(rstanarm)   # 用于后验抽样（注意：安装 rstanarm 可能需要额外配置）
library(ggplot2)
library(ggrepel)
library(forcats)
library(dplyr)
library(tidyr)
library(purrr)
library(broom)
 
tidymodels_prefer()
 
data(ames)
ames <- mutate(ames, Sale_Price = log10(Sale_Price))
 
set.seed(502)
ames_split <- initial_split(ames, prop = 0.80, strata = Sale_Price)
ames_train <- training(ames_split)
ames_test  <- testing(ames_split)

1. 三套预处理（recipes）

我们实现三种预处理方案：basic、interact（加入交互）和 splines（在经纬度上加自然样条）。

basic_rec <- 
  recipe(Sale_Price ~ Neighborhood + Gr_Liv_Area + Year_Built + Bldg_Type + 
           Latitude + Longitude, data = ames_train) %>%
  step_log(Gr_Liv_Area, base = 10) %>% 
  step_other(Neighborhood, threshold = 0.01) %>% 
  step_dummy(all_nominal_predictors())
 
interaction_rec <- 
  basic_rec %>% 
  step_interact( ~ Gr_Liv_Area:starts_with("Bldg_Type_") ) 
 
spline_rec <- 
  interaction_rec %>% 
  step_ns(Latitude, Longitude, deg_free = 50)
 
preproc <- list(basic = basic_rec, 
                interact = interaction_rec, 
                splines = spline_rec)

2. 线性回归工作流集合与 10 折交叉验证

lm_models <- workflow_set(preproc, list(lm = linear_reg()), cross = FALSE)
 
set.seed(1001)
ames_folds <- vfold_cv(ames_train, v = 10)  # 10 折 CV
 
keep_pred <- control_resamples(save_pred = TRUE, save_workflow = TRUE)
 
lm_models <- 
  lm_models %>% 
  workflow_map("fit_resamples", 
               seed = 1101, verbose = TRUE,
               resamples = ames_folds, control = keep_pred)

3. 随机森林（作为比较）

rf_model <-
  rand_forest(trees = 1000) %>%
  set_engine("ranger") %>%
  set_mode("regression")
 
rf_wflow <-
  workflow() %>%
  add_formula(
    Sale_Price ~ Neighborhood + Gr_Liv_Area + Year_Built +
      Bldg_Type + Latitude + Longitude) %>%
  add_model(rf_model)
 
set.seed(1003)
rf_res <- rf_wflow %>% 
  fit_resamples(resamples = ames_folds, control = keep_pred)
 
four_models <- 
  as_workflow_set(random_forest = rf_res) %>% 
  bind_rows(lm_models)

4. 收集每一折的 R²（用于后续比较）

rsq_indiv_estimates <- 
  collect_metrics(four_models, summarize = FALSE) %>% 
  filter(.metric == "rsq")
 
# 宽格式：每一行 (id) = 一次 resample，列为不同模型的 R^2
rsq_wider <- 
  rsq_indiv_estimates %>% 
  select(wflow_id, .estimate, id) %>% 
  pivot_wider(id_cols = "id", names_from = "wflow_id", values_from = ".estimate")

说明：rsq_indiv_estimates 保留了 每一折（每个 id）上所有模型的 R²，而不是折平均。这样可以做配对比较（paired），因为在同一折上评估各模型，它们的表现是可直接比较的。

5. 频率学检验（相关性、配对 t 检验 / 回归）

# 相关性（说明：同一折上各模型的 R^2 通常高度相关）
corrr::correlate(rsq_wider %>% select(-id), quiet = TRUE)
 
# 可视化：每个 resample(id) 上的模型 R^2 路径图
rsq_indiv_estimates %>% 
  mutate(wflow_id = reorder(wflow_id, .estimate)) %>% 
  ggplot(aes(x = wflow_id, y = .estimate, group = id, color = id)) + 
  geom_line(alpha = .5, linewidth = 1.25) +
  theme_minimal() + theme(legend.position = "none")
 
# 配对相关检验（示例：basic_lm vs splines_lm）
rsq_wider %>% with(cor.test(basic_lm, splines_lm)) %>% broom::tidy()
 
# 配对差值分析（平均差值是否为 0）
compare_lm <- rsq_wider %>% mutate(difference = splines_lm - basic_lm)
 
# 1) 用 lm 测试差值均值
lm(difference ~ 1, data = compare_lm) %>% tidy(conf.int = TRUE) %>% 
  select(estimate, p.value, starts_with("conf"))
 
# 2) 或者用配对 t 检验
rsq_wider %>% with(t.test(splines_lm, basic_lm, paired = TRUE)) %>% tidy() %>% 
  select(estimate, p.value, starts_with("conf"))

解读（频率学）：你会得到一个点估计（均值差 ≈ 0.00913）和一个 p 值（非常小），说明在样本上 splines_lm 的均值 R² 显著高于 basic_lm。但这里仍然要问：这个差值在实际上重要吗？ —— 这正是我们要用贝叶斯方法深入判断的地方。

6. 贝叶斯方法（tidyposterior + rstanarm）：为什么用它？

  1. 交叉验证产生的每折结果并非独立（同一折上不同模型的结果相关），频率学上的 t 检验/置信区间常常忽略这种配对结构或不能充分建模它，导致对不确定性估计不够保守或不够准确。
  2. 频率学结果通常给出是否“显著”（p 值），但**无法直接提供“模型 A 比模型 B 更好的概率”**这类对决策更有价值的量化结论。
  3. 我们还关心“差异是否在实践上可忽略（practical equivalence）”，即 ROPE（Region of Practical Equivalence）。

**tidyposterior ：

  • 把每个 resample（fold）看作一个观测单元，建立层次贝叶斯模型：模型性能 = workflow 固有效应 + resample 的随机偏差（捕捉 resample 相关性），并对参数引入先验。
  • 用 rstanarm 做后验抽样，得到每个 workflow 的 R² 的后验分布（不是仅仅一个点估计）。
  • 基于后验样本可以直接计算两模型差值的后验分布、P(difference > 0)、可信区间和 ROPE 中的质量份额（pract_equiv 等）。

7. 用 tidyposterior 拟合并抽样（得到后验）

rsq_anova <-
  perf_mod(
    four_models,
    metric = "rsq",
    prior_intercept = rstanarm::student_t(df = 1),
    chains = 4,
    iter = 5000,
    seed = 1102
  )
 
# 抽取后验样本（tidy 格式，列名通常是 model / posterior）
model_post <- rsq_anova %>% tidy(seed = 1103)
glimpse(model_post)
# -> 例如： model (workflow id), posterior (每个后验样本的 R^2 数值)
 
 
autoplot(rsq_anova) +
	geom_text_repel(aes(label = workflow), nudge_x = 1/8, nudge_y = 1/100) +
	theme(legend.position = "none")

注意：perf_mod() 内部构建了一个层次模型（workflow-level 参数 + resample-level随机效应），这可以把相同 resample 上不同 workflow 的相关性纳入模型，从而得到更合理的不确定性估计（比简单 t 检验更稳健）。

8. 可视化：每个模型的后验分布（R²）

（默认 posterior 的范围通常在 0–1 之间，因为我们在建模 R²）

# 密度图 + 标注均值与 95% 区间
model_post %>% 
	mutate(model = forcats::fct_inorder(model)) %>%  mutate(model = forcats::fct_inorder(model)) %>%
	ggplot(aes(x = posterior)) + 
	geom_histogram(bins = 50, color = "white", fill = "blue", alpha = 0.4) + 
	facet_wrap(~ model, ncol = 1)

说明：横轴 0–1 表示 R² 可能取值范围；纵轴样本计数。后验分布的“位置”表示哪个模型的 R² 更高，分布的宽窄表示不确定性大小。

9. 模型间差值的后验（contrast_models）

rqs_diff <- contrast_models(
  rsq_anova,
  list_1 = "splines_lm",
  list_2 = "basic_lm",
  seed = 1104
)

这里 rqs_diff 就是 差值的后验分布样本池，每个样本 = 一次 MCMC 抽样下 splines_lm 与 basic_lm 的 R² 差值。

差值直方图可视化

rqs_diff %>% 
  as_tibble() %>% 
  ggplot(aes(x = difference)) + 
  geom_vline(xintercept = 0, lty = 2) + 
  geom_histogram(bins = 50, color = "white", fill = "red", alpha = 0.4)
  • 横坐标：差值 (splines_lm - basic_lm)。
  • 纵坐标：每个区间的后验样本数（count）。
  • 竖线 0：差值为零的参考线。

如果大多数后验质量在 >0 的区域，说明 splines_lm 的 R² 高于 basic_lm。

  1. summary(rqs_diff) %>% select(-starts_with("pract"))

这是标准的贝叶斯后验摘要：

contrastprobabilitymeanloweruppersize
splines_lm vs basic_lm1.000.009130.005100.01320

解释：

  • mean = 0.00913 → 差值的后验均值。
  • lower, upper → 95% 可信区间，大约 [0.005, 0.013]。
  • probability = 1.00 → 后验概率 P(diff > 0) ≈ 100%。
  • size = 0 → 因为你没指定 ROPE，所以这一列显示默认值。

👉 这是在问：差值是否大于零？差不多是多少？
结论：splines_lm 比 basic_lm 更好，平均高 0.009，几乎肯定 >0。

② summary(rqs_diff, size = 0.02) %>% select(contrast, starts_with("pract"))

这是基于 ROPE 的“实践等效”分析：

contrastpract_negpract_equivpract_pos
splines_lm vs basic_lm010

解释：

  • pract_neg = 0 → 没有证据支持 splines_lm 明显更差。
  • pract_pos = 0 → 没有证据支持 splines_lm 明显更好（大于 +0.02）。
  • pract_equiv = 1 → 100% 的后验质量落在 [-0.02, +0.02] 区间内。

这是在问：差值是否大到值得在意？
结论：虽然差值 >0 的概率是 100%，但提升幅度太小（<0.02），在实践上可以看作等效。

两者的差异（重点）

  • 不带 size → “统计角度” → 差值几乎肯定 >0。

  • 带 size → “实践角度” → 差值虽然显著，但幅度不足以超出 ±0.02 的等效区间。

  • 第一份 summary = splines_lm 确实更好（统计上显著）。

  • 第二份 summary = 但更好得不够多，可以忽略（实践上等效）。

关系图谱

反向链接