R: Caret包: Brier分数

Question

我想使用caret包中的train()函数执行逻辑回归。我的模型看起来像这样：

model <- train(Y ~.,
  data = train_data,
  family = "binomial",
  method = "glmnet")

使用得到的模型，我想要做出预测：

pred <- predict(model, newdata = test_data, s = "lambda.min", type = "prob")

现在，我想评估模型预测与实际测试数据相比有多好。为此，我知道如何接收ROC和AUC。然而，我也有兴趣收到布瑞尔的分数。Brier分数的公式与MSE几乎相同。我面临的问题是，预测中的type参数只允许"prob“(或我不感兴趣的”类“)，这给出了一个预测为1的概率(例如0.64)，以及为零的互补概率(例如0.37)。然而，对于Brier分数，我需要为包含两者信息的每个预测提供一个概率估计(例如，高于0.5的值表示1，低于0.5的值表示0)。我还没有找到在caret包中接收Brier分数的任何解决方案。我知道使用cv.glmnet包时，predict函数允许参数"response“，这将解决我的问题。然而，出于个人喜好，我想继续使用caret包。谢谢你的帮助！

Answer 1

如果我们使用维基对brier score的定义：

Brier分数最常见的公式是

​

​

其中f_t是预测的概率，o_t是(0或1)的实际结果，N是预测实例的数量。

在R中，如果你的标签是一个因子，那么逻辑回归将总是相对于第二个水平进行预测，这意味着你只需要计算概率和0/1。例如：

library(caret)
idx = sample(nrow(iris),100)
data = iris
data$Species = factor(ifelse(data$Species=="versicolor","v","o"))
levels(data$Species)
[1] "o" "v"

在本例中，o为0，v为1。

train_data = data[idx,]
test_data = data[-idx,]

model <- train(Species ~.,data = train_data,family = "binomial",method = "glmnet")

pred <- predict(model, newdata = test_data)

所以我们可以看到这个类的概率：

head(pred)
          o          v
1 0.8367885 0.16321154
2 0.7970508 0.20294924
3 0.6383656 0.36163437
4 0.9510763 0.04892370
5 0.9370721 0.06292789

计算分数的步骤：

f_t = pred[,2]
o_t = as.numeric(test_data$Species)-1
mean((f_t - o_t)^2)
[1] 0.32

Answer 2

我使用Brier分数在caret中调优我的模型以进行二进制分类。我确保“积极”类是第二个类，这是您将响应标记为"0:1“时的默认值。然后，我基于caret自己的汇总函数套件创建了这个主汇总函数，以返回我想要查看的所有指标：

BigSummary <- function (data, lev = NULL, model = NULL) {
  pr_auc <- try(MLmetrics::PRAUC(data[, lev[2]],
                                 ifelse(data$obs == lev[2], 1, 0)),
                silent = TRUE)
  brscore <- try(mean((data[, lev[2]] - ifelse(data$obs == lev[2], 1, 0)) ^ 2),
               silent = TRUE)
  rocObject <- try(pROC::roc(ifelse(data$obs == lev[2], 1, 0), data[, lev[2]],
                             direction = "<", quiet = TRUE), silent = TRUE)
  if (inherits(pr_auc, "try-error")) pr_auc <- NA
  if (inherits(brscore, "try-error")) brscore <- NA
  rocAUC <- if (inherits(rocObject, "try-error")) {
    NA
  } else {
    rocObject$auc
  }
  tmp <- unlist(e1071::classAgreement(table(data$obs,
                                            data$pred)))[c("diag", "kappa")]
  out <- c(Acc = tmp[[1]],
           Kappa = tmp[[2]],
           AUCROC = rocAUC,
           AUCPR = pr_auc,
           Brier = brscore,
           Precision = caret:::precision.default(data = data$pred,
                                                 reference = data$obs,
                                                 relevant = lev[2]),
           Recall = caret:::recall.default(data = data$pred,
                                           reference = data$obs,
                                           relevant = lev[2]),
           F = caret:::F_meas.default(data = data$pred, reference = data$obs,
                                      relevant = lev[2]))
  out
}

现在我可以简单地在trainControl中传递summaryFunction = BigSummary，然后在train调用中传递metric = "Brier", maximize = FALSE。