ROC与R中线性回归的交叉验证

Question

一个由两部分组成的问题：我试图弄清楚：
(1) 如何使用 lm() 生成线性回归的 ROC 曲线（如果正确的话？？），以及 (2) 如何通过 k 折交叉验证来实现，这样我就可以得到平均 ROC 曲线（和 AUC）。

如果结果是连续变量，它必须转换为二进制变量，对吗？通常我会使用glm(..., family = 'binomial') 来拟合逻辑回归模型，但这是最合适的方法吗？ （看来我只是在拟合不同的模型。）

我想从cvAUC 包的rdrr.io website 中得到类似下面的图（红线是平均 ROC 曲线，虚线是 k 倍 ROC 曲线），但我不知道如何到达那里我的数据。

以data(USArrests) 为例：

library(dplyr)
library(pROC)
data(USArrests)

# create train and test sets
set.seed(2021)
dat <- mutate(USArrests, index=1:nrow(USArrests))
train.dat <- sample_frac(dat, 0.5) # splits `dat` in half
test.dat <- subset(dat, !dat$index %in% train.dat$index) # uses other half to test

# trying to build predictions with lm()
fit <- lm(Murder ~ Assault, data = train.dat)
predicted <- predict(fit, test.dat, type = "response")

# roc curve
roc(test.dat$Murder ~ predicted, plot = TRUE, print.auc = TRUE) # AUC = 1.000

上面的代码得到了结果，但是给出了警告：

警告信息： 在 roc.default(response, m[[predictors]], ...) 中： “响应”有两个以上的层次。考虑明确设置“级别”或改用“multiclass.roc”

我不知道该怎么做。它还得到了 AUC = 1.000——这种方法是否错误，为什么？

此外，它仅适用于一个训练/测试集。我不确定如何使用 k 折集进行训练。我想我必须以某种方式将它与caret::train() 结合起来。我尝试使用来自ROC curve from training data in caret 的随机森林模型的 ROC 解决方案，但它不适用于我的代码。

例子：

library(caret)
library(MLeval)

train_control <- trainControl(method = "cv", number = 10, savePredictions = TRUE)
rfFit <- train(Murder ~ Assault, data = USArrests, trControl = train_control, method = "lm")

rfFit$pred$mtry # NULL
res <- MLeval::evalm(rfFit) # error with error message below

MLeval：机器学习模型评估
输入：插入符号训练函数对象
不平均概率。
第 1 组类型：cv
[.data.frame(preds, c(G1, G2, "obs")) 中的错误： 选择了未定义的列

Answer 1

如果将其切换为 0/1 变量，则可以像这样进行交叉验证：

USArrests <- USArrests %>% 
  mutate(Murder01 = as.numeric(Murder > mean(Murder, na.rm=TRUE)))

# create train and test sets
set.seed(2021)
cvfun <- function(split, ...){
  mod <- glm(Murder01 ~ Assault, data=analysis(split), family=binomial)
  fit <- predict(mod, newdata=assessment(split), type="response")
  data.frame(fit = fit, y = model.response(model.frame(formula(mod), data=assessment(split))))
}
library(rsample)
library(purrr)
library(tidyverse)
cv_out <- vfold_cv(USArrests, v=10, repeats = 5) %>% 
    mutate(fit = map(splits, cvfun)) %>% 
    unnest(fit) %>% 
    group_by(id) %>% 
    summarise(auc = roc(y, fit, plot=FALSE)$auc[1])

cv_out
# # A tibble: 5 x 2
#   id        auc
# * <chr>   <dbl>
# 1 Repeat1 0.936
# 2 Repeat2 0.928
# 3 Repeat3 0.937
# 4 Repeat4 0.918
# 5 Repeat5 0.942

也就是说，我不确定这是否比在线性模型上使用 R 平方或 MSE 更好。而且，我对本教程中的算法实际上在做一些统计上有意义的事情并不是非常有信心。我肯定是错的，并且会听从具有更多专业知识的人，但我看不出它有多大意义，而且在这种情况下它肯定不会产生有意义的东西。您认为只有完美预测才会出现 AUC 为 1。

此外，我不确定这些数字有什么证明价值。通常，您会希望使用这种分析来调整模型规范——通常是通过找到接近最优的超参数值。您可以想象使用不同的模型规范来执行此操作。例如，您可以评估具有 Assault 中的二次多项式的模型与线性模型的相对预测能力，如下所示。

cvfun2 <- function(split, ...){
  mod <- glm(Murder01 ~ poly(Assault, 2),  data=analysis(split), family=binomial)
  fit <- predict(mod, newdata=assessment(split), type="response")
  data.frame(fit = fit, y = model.response(model.frame(formula(mod), data=assessment(split))))
}

cv_out2 <- vfold_cv(USArrests, v=10, repeats = 5) %>% 
    mutate(fit = map(splits, cvfun2)) %>% 
    unnest(fit) %>% 
    group_by(id) %>% 
    summarise(auc = roc(y, fit, plot=FALSE)$auc[1])

mean(cv_out2$auc)
# [1] 0.9123994
mean(cv_out$auc)
# [1] 0.9320451

---

编辑 - 制作 ROC 图

cv_out_plot <- vfold_cv(USArrests, v=10, repeats = 5) %>% 
  mutate(fit = map(splits, cvfun)) %>% 
  unnest(fit) %>% 
  group_by(id) %>% 
  summarise(sens = roc(y, fit, plot=FALSE)$sensitivities, 
         spec = roc(y, fit, plot=FALSE)$specificities, 
         obs = 1:length(sens))
ave <- cv_out_plot %>% 
  ungroup %>% 
  group_by(obs) %>% 
  summarise(sens = mean(sens), 
            spec = mean(spec), 
            id = "Average")

cv_out_plot <- bind_rows(cv_out_plot, ave) %>% 
  mutate(col = factor(ifelse(id == "Average", "Average", "Individual"), 
                      levels=c("Individual", "Average")))



ggplot(cv_out_plot , aes(x=1-sens, y=spec, group=id, colour=col)) + 
  geom_line(aes(size=col, alpha=col)) + 
  scale_colour_manual(values=c("black", "red")) + 
  scale_size_manual(values=c(.5,1.25)) + 
  scale_alpha_manual(values=c(.3, 1)) + 
  theme_classic() + 
  theme(legend.position=c(.75, .15)) + 
  labs(x="1-Sensitivity", y="Specificity", colour="", alpha="", size="")