R语言tidymodels那些事（十六） 模型预测的真实性

缘由

艺多不压身

碎碎念

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

Ref

19 When Should You Trust Your Predictions? | Tidy Modeling with R

解读

本文通过四个部分展示 tidymodels 的建模与诊断流程，重点说明 为什么要引入模糊预测、特征工程、未来数据预测与可解释性分析。

---

一、二分类模拟 + 等概率区间

模拟数据

我们首先模拟一个带噪声的二分类数据集。这里故意设置特征间有相关性，并在生成标签时引入误差，从而让模型面临“不确定”的情况。

library(tidymodels)
tidymodels_prefer()
 
 
simulate_two_classes <- 
  function (n, error = 0.1, eqn = quote(-1 - 2 * x - 0.2 * x^2 + 2 * y^2))  {
    sigma <- matrix(c(1, 0.7, 0.7, 1), nrow = 2, ncol = 2)
    dat <- MASS::mvrnorm(n = n, mu = c(0, 0), Sigma = sigma)
    colnames(dat) <- c("x", "y")
    cls <- paste0("class_", 1:2)
    dat <- 
      as_tibble(dat) %>% 
      mutate(
        linear_pred = !!eqn,
        linear_pred = linear_pred + rnorm(n, sd = error),
        prob = binomial()$linkinv(linear_pred),
        class = ifelse(prob > runif(n), cls[1], cls[2]),
        class = factor(class, levels = cls)
      )
    dplyr::select(dat, x, y, class)
  }
 
set.seed(1901)
training_set <- simulate_two_classes(200)
testing_set  <- simulate_two_classes(50)

为什么这样做：

• 模拟数据比真实数据更可控，可以演示“分类困难”场景。
• 加入噪声后，模型预测不会百分百准确，这正是现实问题的本质。

---

拟合逻辑回归模型

使用 tidymodels 拟合一个带二次项的逻辑回归模型。

two_class_mod <- 
  logistic_reg() %>% 
  set_engine("stan", seed = 1902) %>% 
  fit(class ~ . + I(x^2)+ I(y^2), data = training_set)
 
print(two_class_mod, digits = 3)

优缺点：

• 逻辑回归直观、可解释，但对非线性边界的刻画有限。
• 加入多项式项能提高拟合能力，但也容易过拟合。

---

常规预测 vs 模糊预测

先做常规预测：

test_pred <- augment(two_class_mod, testing_set)
predict(two_class_mod, testing_set)

然后引入 probably::make_two_class_pred()，为预测加上 模糊区间 (Equivocal Zone)：

library(probably)
 
training_set$class=factor(training_set$class,level=c("class_1","class_2"))
lvls <- levels(training_set$class)
 
test_pred <- 
  test_pred %>% 
  mutate(.pred_with_eqz = make_two_class_pred(.pred_class_1, lvls, buffer = 0.15))

为什么要模糊预测：

• 常规预测：即使模型置信度很低，也会强行给出一个类别 → 误导风险高。
• 模糊预测：引入“不可报告”区间，当概率落在模糊区时提示“结果不确定” → 更符合真实业务决策。
• 缺点是：有时会牺牲预测覆盖率。

---

对比结果

test_pred %>% conf_mat(class, .pred_class)
test_pred %>% conf_mat(class, .pred_with_eqz)

对比发现：

• 常规预测准确率较低，因为强行分类。
• 模糊预测虽然减少了报告样本，但整体准确率更高。

我们还可以系统评估不同 buffer 大小的影响：

eq_zone_results <- function(buffer) {
  test_pred <- 
    test_pred %>% 
    mutate(.pred_with_eqz = make_two_class_pred(.pred_class_1, lvls, buffer = buffer))
  acc <- test_pred %>% accuracy(class, .pred_with_eqz)
  rep_rate <- reportable_rate(test_pred$.pred_with_eqz)
  tibble(accuracy = acc$.estimate, reportable = rep_rate, buffer = buffer)
}
 
map(seq(0, .1, length.out = 40), eq_zone_results) %>% 
  list_rbind() %>% 
  pivot_longer(c(-buffer), names_to = "statistic", values_to = "value") %>% 
  ggplot(aes(x = buffer, y = value, lty = statistic)) + 
  geom_step(linewidth = 1.2, alpha = 0.8) + 
  labs(y = NULL, lty = NULL)

下面输出每个分类的预测结果

test_pred <- 
	test_pred %>% 
	bind_cols(
		predict(two_class_mod, testing_set, type = "pred_int", std_error = TRUE)
	)
#输出
 .pred_class .pred_class_1 .pred_class_2      x      y class  
   <fct>               <dbl>         <dbl>  <dbl>  <dbl> <fct>  
 1 class_2           0.0256          0.974  1.12  -0.176 class_2
 2 class_1           0.555           0.445 -0.126 -0.582 class_2
 3 class_2           0.00620         0.994  1.92   0.615 class_2
 4 class_2           0.472           0.528 -0.400  0.252 class_2
 5 class_2           0.163           0.837  1.30   1.09  class_1

---

二、Chicago 数据预测

数据准备

真实数据，目标是预测芝加哥地铁的日均客流。
我们选取 ridership（客流量）、date 和各个车站数据。拆分训练集（去掉最后 14 天）和测试集（最后 14 天）。

data(Chicago)
 
Chicago <- Chicago %>% select(ridership, date, one_of(stations))
n <- nrow(Chicago)
 
Chicago_train <- Chicago %>% slice(1:(n - 14))
Chicago_test  <- Chicago %>% slice((n - 13):n)

特征工程 + PLS 降维

使用 recipe 做：

• 日期特征；
• 节假日变量；
• 虚拟变量转换；
• 去掉常量列；
• 对车站客流标准化 + PLS 降维（减少维度，提高效率）。

base_recipe <-
  recipe(ridership ~ ., data = Chicago_train) %>%
  step_date(date) %>%
  step_holiday(date, keep_original_cols = FALSE) %>%
  step_dummy(all_nominal()) %>%
  step_zv(all_predictors()) %>%
  step_normalize(!!!stations)%>%
  step_pls(!!!stations, num_comp = 10, outcome = vars(ridership))

优缺点：

• 优点：有效压缩高维输入、缓解多重共线性。
• 缺点：PLS 牺牲了解释性（不能直观解释每个车站的贡献）。

---

建模与预测

lm_spec <-
  linear_reg() %>%
  set_engine("lm") 
 
lm_wflow <-
  workflow() %>%
  add_recipe(base_recipe) %>%
  add_model(lm_spec)
 
set.seed(1902)
lm_fit <- fit(lm_wflow, data = Chicago_train)
 
res_test <-
  predict(lm_fit, Chicago_test) %>%
  bind_cols(
    predict(lm_fit, Chicago_test, type = "pred_int"),
    Chicago_test
  )
 
res_test %>% select(date, ridership, starts_with(".pred"))
res_test %>% rmse(ridership, .pred)

• RMSE 较低，说明模型预测合理。
• 提供置信区间，能更好支持风险控制。

---

三、未来数据预测：Chicago_2020

我们人为创建 Chicago_2020，日期设为 2020 年 6 月，保持结构一致。这样可以模拟模型在未来的预测表现。

Chicago_2020 <- Chicago_test %>% 
  mutate(date = seq(as.Date("2020-06-01"), by = "day", length.out = 14))
 
res_2020 <-
  predict(lm_fit, Chicago_2020) %>%
  bind_cols(
    predict(lm_fit, Chicago_2020, type = "pred_int"),
    Chicago_2020
  ) 
 
res_2020 %>% select(date, contains(".pred"))
res_2020 %>% rmse(ridership, .pred)

• 优点：在业务中常用于“未来场景模拟”。
• 缺点：若新数据分布与训练集差异过大，预测可能失真。

---

四、模型可解释性诊断：applicable

最后用 applicable 包检验新数据是否与训练数据分布接近。
通过 PCA 投影后计算 马氏距离，数值越大，表示新样本“离群”，预测越可能不可靠。

library(applicable)
pca_stat <- apd_pca(~ ., data = Chicago_train %>% select(one_of(stations)), 
                    threshold = 0.99)
pca_stat
 
#x 轴显示了距离的值，y 轴展示了分布的百分位数。例如，训练集样本中有一半的距离小于 3.7。autoplot(pca_stat, distance) + labs(x = "distance")
 
score(pca_stat, Chicago_test) %>% select(starts_with("distance")) #训练集中的得分
score(pca_stat, Chicago_2020) %>% select(starts_with("distance")) #新的预测中的得分 用于衡量结果的可信性

为什么要做：

• 优点：帮助判断模型是否在“熟悉的区域”预测。
• 缺点：如果阈值过严，可能把合理的新数据误判为异常。

常见问题

1. 如果未来真的出现了“异常趋势”（超出了训练数据的分布，比如疫情期间工作日客流骤降），我们现在做的特征工程和模型是不是就失效了？

未来的模式可能跟过去不同，那么这些特征仍然被模型用，却解释力大大降低。你可以把“异常趋势”理解为 分布外 (out-of-distribution, OOD) 问题，常见应对方式有：

• 异常点监控
  • 在预测时，先用 applicable::apd_pca() 这类工具检测输入是否偏离训练分布。
  • 一旦发现异常，就不要盲目给出预测，而是提示“预测不可靠”

• 鲁棒方法
  • 使用对异常点更不敏感的模型（如基于树的集成模型、分位数回归），它们对极端值的适应性更强。

• 预测区间 + 不确定性估计
  • 不仅给出点预测，还要输出预测区间。
  • 如果异常点来了，区间会快速变宽，提示“预测不确定”。