我一直在使用Rstudio中的caret包中的gbm来寻找故障发生的概率。
我使用了Youden的J统计量来找到最佳分类的阈值，即0.63。现在我该如何使用这个阈值？我认为最好的方法是将这个阈值以某种方式纳入caret中的gbm模型中，以获得更准确的预测，然后在训练数据上重新运行模型？目前它默认使用0.5的阈值，我找不到明显的方法来更新阈值。
另一种方法是，阈值只是用来将测试数据的预测结果划分到正确的类别中？这似乎更直接，但如果假设应该根据新的阈值更新概率，那么我如何在ROC_AUC图上反映这种变化呢？
任何帮助将不胜感激。谢谢
编辑：我正在处理的完整代码如下：
library(datasets)library(caret)library(MLeval)library(dplyr)data(iris)data <- as.data.frame(iris)# create classdata$class <- ifelse(data$Species == "setosa", "yes", "no")# split into train and testtrain <- data %>% sample_frac(.70)test <- data %>% sample_frac(.30)# Set up control function for trainingctrl <- trainControl(method = "cv", number = 5, returnResamp = 'none', summaryFunction = twoClassSummary, classProbs = T, savePredictions = T, verboseIter = F)# Set up trainng grid - this is based on a hyper-parameter tune that was recently donegbmGrid <- expand.grid(interaction.depth = 10, n.trees = 20000, shrinkage = 0.01, n.minobsinnode = 4) # Build a standard classifier using a gradient boosted machineset.seed(5627)gbm_iris <- train(class ~ ., data = train, method = "gbm", metric = "ROC", tuneGrid = gbmGrid, verbose = FALSE, trControl = ctrl)# Calcuate best thresholdscaret::thresholder(gbm_iris, threshold = seq(.01,0.99, by = 0.01), final = TRUE, statistics = "all")pred <- predict(gbm_iris, newdata = test, type = "prob")roc <- evalm(data.frame(pred, test$class))回答：
你的代码中有几个问题。我将使用mlbench包中的PimaIndiansDiabetes数据集，因为它比iris数据集更适合。
首先，用于将数据拆分为训练集和测试集的代码：
train <- data %>% sample_frac(.70)test <- data %>% sample_frac(.30)不适合使用，因为训练集中出现的一些行也会出现在测试集中。
此外，避免使用函数名作为对象名，这将在长远来看为你节省很多麻烦。
data(iris)data <- as.data.frame(iris) #不好的对象名举个例子：
library(caret)library(ModelMetrics)library(dplyr)library(mlbench)data(PimaIndiansDiabetes, package = "mlbench")创建训练集和测试集，你可以使用基础R的sample来抽样行，或者使用caret::createDataPartition。createDataPartition更可取，因为它试图保留响应的分布。
set.seed(123)ind <- createDataPartition(PimaIndiansDiabetes$diabetes, 0.7)tr <- PimaIndiansDiabetes[ind$Resample1,]ts <- PimaIndiansDiabetes[-ind$Resample1,]这样，训练集中的行就不会出现在测试集中。
让我们创建模型：
ctrl <- trainControl(method = "cv", number = 5, returnResamp = 'none', summaryFunction = twoClassSummary, classProbs = T, savePredictions = T, verboseIter = F)gbmGrid <- expand.grid(interaction.depth = 10, n.trees = 200, shrinkage = 0.01, n.minobsinnode = 4) set.seed(5627)gbm_pima <- train(diabetes ~ ., data = tr, method = "gbm", #使用xgboost metric = "ROC", tuneGrid = gbmGrid, verbose = FALSE, trControl = ctrl)为thresholder创建一个概率向量
probs <- seq(.1, 0.9, by = 0.02)ths <- thresholder(gbm_pima, threshold = probs, final = TRUE, statistics = "all")head(ths)Sensitivity Specificity Pos Pred Value Neg Pred Value Precision Recall F1 Prevalence Detection Rate Detection Prevalence1 200 10 0.01 4 0.10 1.000 0.02222222 0.6562315 1.0000000 0.6562315 1.000 0.7924209 0.6510595 0.6510595 0.99220782 200 10 0.01 4 0.12 1.000 0.05213675 0.6633439 1.0000000 0.6633439 1.000 0.7975413 0.6510595 0.6510595 0.98178403 200 10 0.01 4 0.14 0.992 0.05954416 0.6633932 0.8666667 0.6633932 0.992 0.7949393 0.6510595 0.6458647 0.97399184 200 10 0.01 4 0.16 0.984 0.07435897 0.6654277 0.7936508 0.6654277 0.984 0.7936383 0.6510595 0.6406699 0.96360225 200 10 0.01 4 0.18 0.984 0.14188034 0.6821550 0.8750000 0.6821550 0.984 0.8053941 0.6510595 0.6406699 0.94012306 200 10 0.01 4 0.20 0.980 0.17179487 0.6886786 0.8833333 0.6886786 0.980 0.8086204 0.6510595 0.6380725 0.9271018 Balanced Accuracy Accuracy Kappa J Dist1 0.5111111 0.6588517 0.02833828 0.02222222 0.97777782 0.5260684 0.6692755 0.06586592 0.05213675 0.94786323 0.5257721 0.6666781 0.06435166 0.05154416 0.94063574 0.5291795 0.6666781 0.07134190 0.05835897 0.92602505 0.5629402 0.6901572 0.15350721 0.12588034 0.85853086 0.5758974 0.6979836 0.18460584 0.15179487 0.8288729根据你 preferred metric提取阈值概率
ths %>% mutate(prob = probs) %>% filter(J == max(J)) %>% pull(prob) -> thresh_probthresh_prob0.74在测试数据上进行预测
pred <- predict(gbm_pima, newdata = ts, type = "prob")根据测试集中的响应创建一个数值响应（0或1），因为这对于ModelMetrics包中的函数是必需的
real <- as.numeric(factor(ts$diabetes))-1ModelMetrics::sensitivity(real, pred$pos, cutoff = thresh_prob)0.2238806 #根据这一点，很明显在这些测试数据上选择的阈值不是最优的ModelMetrics::specificity(real, pred$pos, cutoff = thresh_prob)0.956ModelMetrics::kappa(real, pred$pos, cutoff = thresh_prob)0.2144026 #根据这一点，很明显在这些测试数据上选择的阈值不是最优的ModelMetrics::mcc(real, pred$pos, cutoff = thresh_prob)0.2776309 #根据这一点，很明显在这些测试数据上选择的阈值不是最优的ModelMetrics::auc(real, pred$pos)0.8047463 #AUC良好，但mcc和kappa低表明阈值选择不佳Auc是所有阈值上的度量，因此不需要指定截止阈值。
由于只使用了一个训练/测试拆分，性能评估将存在偏差。最好使用嵌套重采样，以便在多个训练/测试拆分上进行评估。这里是一种执行嵌套重采样的方法。
编辑：回答评论中的问题。
要创建roc曲线，你不需要在所有阈值上计算敏感性和特异性，你可以只使用指定的包来完成这项任务。结果可能会更加可信。
我更喜欢使用pROC包：
library(pROC)roc.obj <- roc(real, pred$pos)plot(roc.obj, print.thres = "best")
图中的最佳阈值是测试数据上特异性+敏感性最高的阈值。显然，这个阈值（0.289）比基于交叉验证预测获得的阈值（0.74）低得多。这就是我说如果你在交叉验证预测上调整阈值，并将由此获得的性能作为阈值成功的指标，将会存在相当大的乐观偏见的原因。
在上面的例子中，不调整阈值会导致在测试集上表现更好。这可能普遍适用于Pima Indians数据集，或者这可能是一个不幸的训练/测试拆分情况。因此，最好使用嵌套重采样来验证这种情况。
