【二分类因变量机器学习】R语言07-XGBoost机器学习模型（标准化代码）

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

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

# 加载必要的包
if (!require(pacman)) install.packages("pacman")
pacman::p_load(readxl, writexl, xgboost, tidymodels, DALEXtra, auditor, vivo, pROC, 
               ggplot2, officer, flextable, gridExtra, showtext)

# 设置全局中文字体支持
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"
}

# 1. 数据准备
data <- read_excel("示例数据.xlsx")

# 将分类变量转换为因子
data$`肥胖程度` <- as.factor(data$`肥胖程度`)
data$`教育水平` <- as.factor(data$`教育水平`)
data$`血型` <- as.factor(data$`血型`)
data$`指标8` <- as.factor(data$`指标8`)
data$`结局` <- as.factor(data$`结局`)

# 分层抽样划分训练集和测试集
set.seed(123)
train_index <- sample(1:nrow(data), nrow(data)*0.7)
Train <- data[train_index, ]
Test <- data[-train_index, ]

# 保存数据集
write_xlsx(Train, "Results-xgb/Train.xlsx")
write_xlsx(Test, "Results-xgb/Test.xlsx")
write_xlsx(data, "Results-xgb/All.xlsx")

# 2. XGBoost模型训练

# 4. 模型评估
# ROC曲线
roc_result <- roc(Test$结局, xg_pred_prob_value)

# 绘制ROC曲线
jpeg("Results-xgb/ROC_curve.jpg", width = 800, height = 600)
plot(roc_result, main = "XGBoost模型ROC曲线", 
     print.auc = TRUE, print.auc.x = 0.6, print.auc.y = 0.4,
     auc.polygon = TRUE, max.auc.polygon = TRUE, grid = TRUE)
dev.off()

pdf("Results-xgb/ROC_curve.pdf", width = 8, height = 6)
plot(roc_result, main = "XGBoost模型ROC曲线", 
     print.auc = TRUE, print.auc.x = 0.6, print.auc.y = 0.4,
     auc.polygon = TRUE, max.auc.polygon = TRUE, grid = TRUE)
dev.off()

# 计算AUC值
auc_value <- auc(roc_result)

# 混淆矩阵
conf_matrix <- confusionMatrix(xg_pred_class$.pred_class, Test$结局)

# 计算Brier得分
brier_score <- mean((as.numeric(xg_pred_prob_value) - as.numeric(Test$结局))^2)

# 5. 模型解释

pdf("Results-xgb/shap_plot.pdf", width = 8, height = 6)
plot(xgb_shap, show_boxplots = FALSE)
dev.off()

# 模型解释绘制ROC曲线
model_eval <- model_evaluation(xgb_exp)

jpeg("Results-xgb/model_evaluation.jpg", width = 800, height = 600)
plot(model_eval)
dev.off()

pdf("Results-xgb/model_evaluation.pdf", width = 8, height = 6)
plot(model_eval)
dev.off()

# 模型解释部分依赖图（PDP）
xgb_profiles <- model_profile(xgb_exp)

jpeg("Results-xgb/partial_dependence.jpg", width = 800, height = 600)
plot(xgb_profiles)
dev.off()

pdf("Results-xgb/partial_dependence.pdf", width = 8, height = 6)
plot(xgb_profiles)
dev.off()

# 模型解释变量重要性
xgb_vi <- model_parts(xgb_exp)

jpeg("Results-xgb/variable_importance.jpg", width = 800, height = 600)
plot(xgb_vi)
dev.off()

pdf("Results-xgb/variable_importance.pdf", width = 8, height = 6)
plot(xgb_vi)
dev.off()

# 模型解释CP图
xgb_cp <- predict_profile(xgb_exp, new_observation = Test[1, ])

jpeg("Results-xgb/ceteris_paribus.jpg", width = 800, height = 600)
plot(xgb_cp)
dev.off()

pdf("Results-xgb/ceteris_paribus.pdf", width = 8, height = 6)
plot(xgb_cp)
dev.off()

# 6. 保存结果
# 保存预测结果
results_table <- data.frame(
  实际值 = Test$结局,
  预测概率 = xg_pred_prob_value,
  预测类别 = xg_pred_class$.pred_class
)
write_xlsx(results_table, "Results-xgb/predictions.xlsx")

# 保存模型性能指标
performance_table <- data.frame(
  指标 = c("准确率", "灵敏度", "特异度", "AUC", "Brier得分"),
  值 = c(round(conf_matrix$overall[1], 4), 
        round(conf_matrix$byClass[1], 4), 
        round(conf_matrix$byClass[2], 4),
        round(auc_value, 4),
        round(brier_score, 4))
)
write_xlsx(performance_table, "Results-xgb/performance_metrics.xlsx")

# 保存变量重要性数据
vi_data <- as.data.frame(xgb_vi)
write_xlsx(vi_data, "Results-xgb/variable_importance_data.xlsx")

# 7. 创建Word报告
doc <- read_docx()
doc <- body_add_par(doc, "XGBoost模型分析报告", style = "heading 1")
doc <- body_add_par(doc, paste0("分析日期: ", Sys.Date()), style = "Normal")

doc <- body_add_par(doc, "数据概览", style = "heading 2")
doc <- body_add_par(doc, paste0("总样本量: ", nrow(data)), style = "Normal")
doc <- body_add_par(doc, paste0("训练集样本量: ", nrow(Train)), style = "Normal")
doc <- body_add_par(doc, paste0("测试集样本量: ", nrow(Test)), style = "Normal")

doc <- body_add_par(doc, "模型性能指标", style = "heading 2")
doc <- body_add_flextable(doc, flextable(performance_table))

doc <- body_add_par(doc, "ROC曲线", style = "heading 2")
doc <- body_add_img(doc, "Results-xgb/ROC_curve.jpg", width = 6, height = 5)

doc <- body_add_par(doc, "变量重要性", style = "heading 2")
doc <- body_add_img(doc, "Results-xgb/variable_importance.jpg", width = 6, height = 5)

doc <- body_add_par(doc, "模型解释", style = "heading 2")
doc <- body_add_par(doc, "Breakdown分析", style = "heading 3")
doc <- body_add_img(doc, "Results-xgb/breakdown_plot.jpg", width = 6, height = 5)

doc <- body_add_par(doc, "SHAP值分析", style = "heading 3")
doc <- body_add_img(doc, "Results-xgb/shap_plot.jpg", width = 6, height = 5)

doc <- body_add_par(doc, "部分依赖图", style = "heading 3")
doc <- body_add_img(doc, "Results-xgb/partial_dependence.jpg", width = 6, height = 5)

doc <- body_add_par(doc, "Ceteris Paribus分析", style = "heading 3")
doc <- body_add_img(doc, "Results-xgb/ceteris_paribus.jpg", width = 6, height = 5)

doc <- body_add_par(doc, "结论", style = "heading 2")
doc <- body_add_par(doc, paste0(
"XGBoost模型在测试集上表现", 
  ifelse(auc_value > 0.7, "良好", "一般"),
"，AUC值为", round(auc_value, 3), 
"。模型准确率为", round(conf_matrix$overall[1], 3), "。"
), style = "Normal")

print(doc, target = "Results-xgb/XGBoost分析报告.docx")

# 8. 保存工作空间
save.image("Results-xgb/XGBoost分析.RData")

cat("分析完成！所有结果已保存到 Results-xgb 文件夹中。")