【二分类因变量机器学习】R语言11-九种模型五种评价批处理（标准化代码）

# 设置工作目录和清理环境
rm(list = ls())
if (!is.null(dev.list())) dev.off()
setwd("C:/Users/hyy/Desktop/")

# 创建结果文件夹
if (!dir.exists("Results-comp")) dir.create("Results-comp")

# 加载必要的包
if (!require(pacman)) install.packages("pacman")
pacman::p_load(mlr3verse, pROC, ggplot2, patchwork, ggsci, dcurves, probably, showtext,
               tidymodels, caret, writexl, openxlsx, readxl, dplyr, tidyr, ranger, kknn)

# 设置全局中文字体支持
font_add(family = "simhei", regular = "simhei.ttf") # 添加黑体
font_add(family = "simsun", regular = "simsun.ttc") # 添加宋体
showtext_auto(enable = TRUE) # 启用showtext

# 设置图形设备参数（适用于基础绘图）
if (.Platform$OS.type == "windows") {
  windowsFonts(SimHei = windowsFont("SimHei"))
  windowsFonts(SimSun = windowsFont("SimSun"))
  chinese_font <- "SimHei"
} else {
  chinese_font <- "sans"
}

# 一、数据加载与预处理
# 读取Excel数据
data_path <- "示例数据.xlsx"
data <- read_excel(data_path, sheet = "示例数据")

# 只保留结局和指标1-6列（避免因子水平不一致问题）
data <- data %>% dplyr::select(结局, 指标1, 指标2, 指标3, 指标4, 指标5, 指标6)

# 数据预处理：将结局转换为因子
data$结局 <- as.factor(data$结局)

# 查看数据结构
str(data)

# 创建mlr3任务
task <- as_task_classif(data, target = "结局", positive = "1")
set.seed(2025)
split <- partition(task, ratio = 0.7)

# 导出数据集至Excel并查看
write_xlsx(data, "C:/Users/hyy/Desktop/Results-comp/示例数据集.xlsx")

# 二、多模型训练与预测
# 定义8个学习器（确保概率预测）

# 三、单模型ROC绘制函数
plot_single_roc <- function(model_name, color) {
if (is.null(predictions[[model_name]])) return(NULL)

# 检查测试集是否包含两个类别
  test_truth <- task$truth(split$test)
if (length(unique(test_truth)) < 2) {
    message("测试集只包含一个类别，无法计算ROC: ", model_name)
    return(NULL)
  }

  roc_obj <- roc(
    response = as.numeric(test_truth) - 1,
    predictor = predictions[[model_name]]$prob
  )
  ggroc(roc_obj, legacy.axes = TRUE) +
    geom_abline(slope = 1, linetype = "dashed", alpha = 0.5) +
    labs(title = paste(model_name, "ROC曲线"),
         subtitle = paste0("AUC = ", round(roc_obj$auc, 3))) +
    scale_color_manual(values = color) +
    theme_minimal(base_size = 10) +
    theme(legend.position = "none")
}

# 四、多图排版与组合绘制
colors <- pal_npg("nrc")(8)
single_plots <- list()
for (name in names(learners)) {
  plot <- plot_single_roc(name, colors[match(name, names(learners))])
if (!is.null(plot)) {
    single_plots[[name]] <- plot
  }
}

if (length(single_plots) > 0) {
# 动态调整布局
  n_plots <- length(single_plots)
  n_col <- min(4, n_plots)
  n_row <- ceiling(n_plots / n_col)

  combined_plots <- wrap_plots(single_plots, ncol = n_col, nrow = n_row) +
    plot_annotation(tag_levels = 'A')

# 绘制所有模型在同一坐标系下的ROC曲线
  roc_data <- list()
for (name in names(predictions)) {
    test_truth <- task$truth(split$test)
    if (length(unique(test_truth)) >= 2) {
      roc_obj <- roc(as.numeric(test_truth) - 1, predictions[[name]]$prob)
      roc_data[[name]] <- data.frame(
        Model = name,
        FPR = 1 - roc_obj$specificities,
        TPR = roc_obj$sensitivities
      )
    }
  }

if (length(roc_data) > 0) {
    roc_combined <- do.call(rbind, roc_data)
    overlay_plot <- ggplot(roc_combined, aes(FPR, TPR, color = Model)) +
      geom_line(linewidth = 0.8) +
      geom_abline(slope = 1, linetype = "dashed", alpha = 0.5) +
      scale_color_manual(values = colors) +
      labs(x = "1 - Specificity", y = "Sensitivity", 
           title = "多模型ROC曲线对比") +
      theme_minimal(base_size = 12) +
      theme(legend.position = "right")
    
    # 整合所有可视化结果
    final_output <- (combined_plots / overlay_plot) + 
      plot_layout(heights = c(2, 1))
    
    # 输出高清图
    ggsave("C:/Users/hyy/Desktop/Results-comp/roc_combined.png", final_output, 
           width = 16, height = 12, dpi = 300)
  }
}

# 五、DCA曲线分析
truth <- as.numeric(task$truth(split$test)) - 1
dca_data <- data.frame(truth = truth)
for (model in names(predictions)) {
  dca_data[[model]] <- predictions[[model]]$prob
}

# 确保truth是因子（dca要求）
dca_data$truth <- factor(dca_data$truth, levels = c(0, 1))

# 计算DCA曲线数据
dca_result <- dca(truth ~ .,
                  data = dca_data,
                  thresholds = seq(0.01, 0.99, 0.01)
)

# 导出DCA结果
write_xlsx(dca_result[["dca"]], "C:/Users/hyy/Desktop/Results-comp/dca_result.xlsx")

# 可视化DCA曲线
dcaresult <- dca_result[["dca"]] %>%
  filter(!label %in% c("Treat All", "Treat None")) %>%
  droplevels()

single_dca_plot <- ggplot(dcaresult, aes(x = threshold, y = net_benefit)) +
  geom_line(aes(color = label), linewidth = 1) +
  facet_wrap(~label, ncol = 4) +
  scale_color_npg() +
  labs(x = "Threshold Probability", y = "Net Benefit", 
       title = "单模型DCA曲线对比") +
  theme_bw(base_size = 12) +
  theme(legend.position = "none")

overlay_dca_plot <- ggplot(dca_result[["dca"]], aes(x = threshold, y = net_benefit)) +
  geom_line(aes(color = label), linewidth = 1, alpha = 0.8) +
  scale_color_npg() +
  labs(x = "Threshold Probability", y = "Net Benefit",
       title = "多模型DCA曲线叠加对比") +
  theme_bw(base_size = 14) +
  theme(legend.position = "right")

combined_plot <- (single_dca_plot / overlay_dca_plot) + 
  plot_layout(heights = c(2, 1)) +
  plot_annotation(tag_levels = 'A')

ggsave("C:/Users/hyy/Desktop/Results-comp/dca_combined.png", combined_plot, 
       width = 16, height = 12, dpi = 300, bg = "white")

# 六、校准曲线分析
# 批量训练与生成校准数据

# 生成分面数据
facet_data <- list()
for (name in names(cal_data)) {
if (!is.null(cal_data[[name]])) {
    facet_data[[name]] <- data.frame(
      Model = name,
      pred_bin = cal_data[[name]]$pred,
      true_bin = as.numeric(cal_data[[name]]$truth) - 1
    )
  }
}
facet_data <- do.call(rbind, facet_data)

# 分面校准曲线
single_plot <- ggplot(facet_data, aes(x = pred_bin, y = true_bin)) +
  geom_point(color = "#E64B35FF", size = 2) +
  geom_smooth(method = "loess", color = "#4DBBD5FF", linewidth = 1) +
  geom_abline(linetype = "dashed", color = "grey40") +
  facet_wrap(~Model, ncol = 4) +
  labs(x = "预测概率分箱均值", y = "实际正例比例", 
       title = "单模型校准曲线对比") +
  theme_bw(base_size = 12) +
  theme(panel.grid.minor = element_blank())

# 叠加校准曲线
overlay_plot <- ggplot(facet_data, aes(x = pred_bin, y = true_bin, color = Model)) +
  geom_point(size = 2, alpha = 0.8) +
  geom_smooth(method = "loess", linewidth = 1, alpha = 0.6, se = FALSE) +
  geom_abline(linetype = "dashed", color = "grey40") +
  scale_color_npg() +
  labs(x = "预测概率分箱均值", y = "实际正例比例", 
       title = "多模型校准曲线对比") +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom",
        panel.grid.minor = element_blank())

# 多图排版与组合绘制
final_plot <- single_plot / overlay_plot + 
  plot_layout(heights = c(2, 1)) +
  plot_annotation(tag_levels = 'A')

# 输出高清图
ggsave("C:/Users/hyy/Desktop/Results-comp/calibration_combined.png", final_plot, 
       width = 16, height = 12, dpi = 300, bg = "white")

# 七、混淆矩阵分析
# 重新定义学习器（确保分类预测）
learners_response <- list(
"Logistic"   = lrn("classif.log_reg", predict_type = "response"),
"决策树"     = lrn("classif.rpart", predict_type = "response"),
"随机森林"   = lrn("classif.ranger", predict_type = "response"),
"SVM"        = lrn("classif.svm", type = "C-classification"),
"神经网络"   = lrn("classif.nnet", predict_type = "response"),
"XGBoost"    = lrn("classif.xgboost", predict_type = "response"),
"朴素贝叶斯" = lrn("classif.naive_bayes", predict_type = "response"),
"KNN"        = lrn("classif.kknn", predict_type = "response", k = 2)  # 调整k值
)

# 批量训练与预测
predictions_response <- list()
for (model_name in names(learners_response)) {
  tryCatch({
    l <- learners_response[[model_name]]
    l$train(task, split$train)
    pred <- l$predict(task, split$test)
    predictions_response[[model_name]] <- data.frame(
      truth = pred$truth,
      response = pred$response,
      model = l$id
    )
  }, error = function(e) {
    message("响应预测失败: ", model_name, " - ", e$message)
  })
}

# 自定义混淆矩阵绘图函数
plot_conf_matrix <- function(pred_df) {
  cm <- table(pred_df$response, pred_df$truth)
  cm_df <- as.data.frame(cm)
  colnames(cm_df) <- c("Predicted", "Actual", "Freq")

  ggplot(cm_df, aes(x = Actual, y = Predicted, fill = Freq)) +
    geom_tile(color = "white", linewidth = 0.5) +
    geom_text(aes(label = Freq), color = "white", size = 5) +
    scale_fill_material("blue-grey", limits = c(0, max(cm_df$Freq))) +
    labs(title = unique(pred_df$model)) +
    theme_minimal(base_size = 12) +
    theme(panel.grid = element_blank(),
          axis.text.x = element_text(angle = 45, hjust = 1))
}

# 生成单模型图列表
single_plots <- list()
for (model_name in names(predictions_response)) {
if (!is.null(predictions_response[[model_name]])) {
    single_plots[[model_name]] <- plot_conf_matrix(predictions_response[[model_name]])
  }
}

# 组合为2x4布局
if (length(single_plots) > 0) {
  n_plots <- length(single_plots)
  n_col <- min(4, n_plots)
  n_row <- ceiling(n_plots / n_col)

  combined_plots <- wrap_plots(single_plots, ncol = n_col, nrow = n_row) +
    plot_annotation(tag_levels = 'A', 
                    title = "多模型混淆矩阵对比") +
    plot_layout(guides = "collect") & 
    theme(legend.position = "bottom")

# 输出组合图
  ggsave("C:/Users/hyy/Desktop/Results-comp/conf_matrix_combined.png", combined_plots, 
         width = 16, height = 9, dpi = 300)
}

# 八、评估指标批量计算
# 重新训练模型（确保概率预测）
trained_models <- list()
for (model_name in names(learners)) {
  tryCatch({
    l <- learners[[model_name]]
    l$train(task, split$train)
    trained_models[[model_name]] <- l
  }, error = function(e) {
    message("模型训练失败: ", model_name, " - ", e$message)
  })
}

# 生成评估报告
eval_metrics <- list()
for (model_name in names(trained_models)) {
  tryCatch({
    model <- trained_models[[model_name]]
    pred <- model$predict(task, split$test)
    
    # 转换为caret需要的因子格式
    pred_factor <- factor(pred$response, levels = c("1", "0"))
    truth_factor <- factor(pred$truth, levels = c("1", "0"))
    
    # 计算混淆矩阵及指标
    cm <- confusionMatrix(pred_factor, truth_factor, positive = "1")
    
    # 提取核心指标
    eval_metrics[[model_name]] <- data.frame(
      Model = model_name,
      Accuracy = ifelse(is.numeric(cm$overall["Accuracy"]), cm$overall["Accuracy"], NA),
      Sensitivity = ifelse(is.numeric(cm$byClass["Sensitivity"]), cm$byClass["Sensitivity"], NA),
      Specificity = ifelse(is.numeric(cm$byClass["Specificity"]), cm$byClass["Specificity"], NA),
      F1 = ifelse(is.numeric(cm$byClass["F1"]), cm$byClass["F1"], NA),
      PPV = ifelse(is.numeric(cm$byClass["Pos Pred Value"]), cm$byClass["Pos Pred Value"], NA),
      NPV = ifelse(is.numeric(cm$byClass["Neg Pred Value"]), cm$byClass["Neg Pred Value"], NA),
      AUC = tryCatch({
        as.numeric(mlr3measures::auc(truth = pred$truth, prob = pred$prob[, "1"], positive = "1"))
      }, error = function(e) NA)
    )
  }, error = function(e) {
    message("评估指标计算失败: ", model_name, " - ", e$message)
    eval_metrics[[model_name]] <- data.frame(
      Model = model_name,
      Accuracy = NA,
      Sensitivity = NA,
      Specificity = NA,
      F1 = NA,
      PPV = NA,
      NPV = NA,
      AUC = NA
    )
  })
}

eval_metrics <- bind_rows(eval_metrics)

# 格式优化
eval_metrics[, -1] <- round(eval_metrics[, -1], 3)

# 结果输出
print(eval_metrics)
write_xlsx(eval_metrics, "C:/Users/hyy/Desktop/Results-comp/混淆矩阵result.xlsx")

# 生成Excel报告
wb <- createWorkbook()
addWorksheet(wb, "模型评估")
writeData(wb, "模型评估", eval_metrics, startRow = 2, startCol = 2)
setColWidths(wb, "模型评估", cols = 2:9, widths = 15)

# 为数值列添加格式
if (nrow(eval_metrics) > 0) {
  addStyle(wb, "模型评估", 
           style = createStyle(numFmt = "0.000"), 
           rows = 3:(nrow(eval_metrics)+2), 
           cols = 3:9, gridExpand = TRUE)
}

saveWorkbook(wb, "C:/Users/hyy/Desktop/Results-comp/model_evaluation.xlsx", overwrite = TRUE)

# 九、临床影响曲线（CIC）分析
# 批量生成CIC数据
truth <- factor(task$truth(split$test))  # 保持因子类型
cic_list <- list()
for (model_name in names(learners)) {
  tryCatch({
    l <- learners[[model_name]]
    l$train(task, split$train)
    pred <- l$predict(task, split$test)
    
    # 构建DCA数据框
    dca_data <- data.frame(
      truth = factor(pred$truth),
      prob = pred$prob[, "1"],  # 提取正类概率
      model = l$id
    )
    
    # 生成CIC曲线
    cic_result <- dca(truth ~ prob,
                      data = dca_data,
                      thresholds = seq(0, 0.5, 0.01)
    )
    cic_list[[model_name]] <- cic_result$dca
  }, error = function(e) {
    message("CIC分析失败: ", model_name, " - ", e$message)
  })
}

# 单模型绘图模板
plot_cic_single <- function(data, model_name) {
  ggplot(data, aes(x = threshold)) +
    geom_line(aes(y = net_benefit), color = "#2E8B57", linewidth = 1) +
    labs(title = model_name, x = "风险阈值", y = "净获益") +
    theme_bw() +
    theme(panel.grid.minor = element_blank(),
          plot.title = element_text(hjust = 0.5))
}

# 生成单模型图列表
single_plots <- list()
for (model_name in names(cic_list)) {
if (!is.null(cic_list[[model_name]])) {
    single_plots[[model_name]] <- plot_cic_single(cic_list[[model_name]], model_name)
  }
}

# 布局组合图
if (length(single_plots) > 0) {
  n_plots <- length(single_plots)
  n_col <- min(4, n_plots)
  n_row <- ceiling(n_plots / n_col)

  combined_plot <- wrap_plots(single_plots, ncol = n_col, nrow = n_row) +
    plot_annotation(title = "多模型临床影响曲线对比", tag_levels = "A") +
    plot_layout(guides = "collect")

# 叠加对比图
  overlay_data <- bind_rows(cic_list, .id = "model")
  overlay_plot <- ggplot(overlay_data, aes(x = threshold, y = net_benefit, color = model)) +
    geom_line(linewidth = 1, alpha = 0.8) +
    scale_color_npg() +
    labs(x = "风险阈值", y = "净获益") +
    theme_bw() +
    theme(legend.position = "bottom")

# 最终输出布局
  final_plot <- combined_plot / overlay_plot + 
    plot_layout(heights = c(2, 1))

# 图形保存
  ggsave("C:/Users/hyy/Desktop/Results-comp/CIC_Combined.png", final_plot, 
         width = 16, height = 12, dpi = 300)
}

# 输出所有结果完成消息
message("分析完成！所有结果已保存到 Results-comp 文件夹中")