具有 R 中线性回归交叉验证的 ROC
ROC with cross-validation for linear regression in R
一个由两部分组成的问题：我想弄清楚：
(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
上面的代码得到了结果，但是给出了警告：
Warning message:
In roc.default(response, m[[predictors]], ...) :
'response' has more than two levels. Consider setting 'levels' explicitly or using 'multiclass.roc' instead
根据它的建议，我不知道该怎么办。它还得到了 AUC = 1.000——这种方法是否错误，为什么？
此外，它只适用于一组 train/test。我不确定如何使用 k 折集进行训练。我想我必须以某种方式将它与 caret::train() 结合起来。我尝试使用来自 的随机森林模型的 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: Machine Learning Model Evaluation
Input: caret train function object
Not averaging probs.
Group 1 type: cv
Error in [.data.frame(preds, c(G1, G2, "obs")) :
undefined columns selected
如果将其切换为 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="")
一个由两部分组成的问题：我想弄清楚：
(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
上面的代码得到了结果，但是给出了警告：
Warning message: In roc.default(response, m[[predictors]], ...) : 'response' has more than two levels. Consider setting 'levels' explicitly or using 'multiclass.roc' instead
根据它的建议，我不知道该怎么办。它还得到了 AUC = 1.000——这种方法是否错误，为什么？
此外，它只适用于一组 train/test。我不确定如何使用 k 折集进行训练。我想我必须以某种方式将它与 caret::train() 结合起来。我尝试使用来自
示例：
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: Machine Learning Model Evaluation
Input: caret train function object
Not averaging probs.
Group 1 type: cv
Error in[.data.frame(preds, c(G1, G2, "obs")) : undefined columns selected
如果将其切换为 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="")