R语言医学分析实战：心脏术后复发预测模型构建全流程

简介：一套面向医学数据分析与临床研究人员的R语言实战教程，聚焦心脏病术后复发风险预测这一典型场景，帮助读者从数据导入、探索性分析到建模评估完整走通流程。教程覆盖tidyverse数据处理、缺失值处理与类别不平衡方案，并讲解逻辑回归、随机森林、XGBoost等算法在Caret框架下的交叉验证与性能对比，结合ROC曲线和AUC值评估模型优劣，同时介绍变量重要性与模型解释方法。压缩包共906个文件，约41.52MB，以R脚本、RMarkdown源文件、HTML报告、CSV数据及PNG图表为主，另有少量PDF、XLSX及RDS对象文件，便于直接复现代码与查看结果。资源已吸引585人学习，适合具备基础统计知识或有一定R使用经验、希望用真实医学数据练习预测建模的读者。通过该教程可掌握完整的医学预测建模思路，获得可运行的脚本、数据预处理模板与可视化输出，为后续临床决策支持分析打下基础。

1. 心脏病术后复发预测：R医学分析怎么把“经验”变成“概率”

心内科查房时最常听到的对话是“这个病人术后半年大概率会复发，盯紧点”，但你再追问一句“大概率是多大概率”，答案就变成了“经验感觉”。这种事在临床数据里其实是可以算的。用 R 语言做医学分析，把心脏术后患者的病历、检验、手术记录和随访结果整理成结构化数据，训练一个复发预测模型，最终输出的是一个 0 到 1 之间的概率，而不是一句“感觉危险”。这篇笔记就讲一条完整可落地的路径：数据怎么整理、特征怎么构造、模型怎么选、结果怎么评估，以及那些让模型翻车的隐藏坑。适合正在做临床课题、想发预测模型类论文，或者科室内部想建术后预警机制的从业者，从零跑通到能答辩。

2. 术后复发预测的数据基础：先弄清手头有什么，再谈建模

2.1 复发预测建模到底需要哪些数据字段

先说结论：不是数据越多越好，是跟“复发”有因果链的字段才值得进模型。心脏术后复发的预测，常见的数据来源是三个系统——HIS 病历系统、检验科的 LIS 系统、随访门诊或电话随访记录。把这些数据导出来之后，字段大致分四类。

这里最容易犯的错是把“术后第 7 天指标”和“术前特征”混在一个模型里预测“是否复发”。后面避坑章会单独展开，现在你只需要记住一个原则：建模时要明确预测的时间原点，比如“从手术日开始的 12 个月内是否复发”，那所有特征必须是时间原点之前可获得的信息。

2.2 用 dplyr 完成清洗和衍生变量：一组可以直接改用的 R 代码

数据导出来通常是 CSV 或者 Excel，Excel 先另存为 CSV 再进 R，避免中文编码问题。下面这组代码是数据清洗的骨架，我一般会按项目改字段名后直接跑。

library(tidyverse) df <- read_csv("cardio_postop.csv", na = c("", "NA", "unknown")) df <- df %>% mutate( age_years = as.numeric(age), bmi = weight_kg / (height_cm / 100)^2, egfr = ifelse(creatinine > 0, 186 * (creatinine / 88.4)^(-1.154) * (age_years)^(-0.203) * ifelse(sex == "female", 0.742, 1), NA_real_), hypertension = ifelse(hypertension == "是" | hypertension == "1", 1, 0), diabetes = ifelse(diabetes == "是" | diabetes == "1", 1, 0), lvef_group = cut(lvef, breaks = c(0, 40, 50, 100), labels = c("重度", "中度", "正常")) ) %>% filter(!is.na(lvef), !is.na(relapse)) write_csv(df, "cardio_clean.csv")

代码逻辑分三段：第一部分读入数据并统一缺失值表示；第二部分用 mutate 做衍生变量，其中 BMI 和 eGFR 是连续性变量，高血压和糖尿病转成二分类 0/1，LVEF 按临床常用阈值切成三档，方便后续看非线性关系；第三部分过滤掉核心变量缺失的行并写出清洗后数据。

需要留意两个参数设置。一是 eGFR 公式，不同科室习惯用 CKD-EPI 或 MDRD，上面写的是简化 MDRD 公式，单位与国内检验结果对应，实际使用时换成你科室的惯例公式即可。二是 lvef_group 的切点不要照抄，40% 和 50% 是心衰指南里常用界值，如果你的样本量小，也可以把 LVEF 保留为连续变量，后面建模时用限制性立方样条处理非线性。

2.3 复发标签的定义：没想清楚这个，后面全白做

标签是预测模型的终点，也是所有前期工作的地基。心脏术后“复发”在不同研究里定义差异很大：有的定义为心绞痛复发，有的定义为再入院，有的定义为再次血运重建，有的定义为心源性死亡。你在建模前必须做的不是跑代码，而是和科室对清楚“复发”的操作性定义。

我建议的标签设计是三列：follow_up_months是实际随访时长，relapse是随访期内是否发生目标事件，time_to_event是手术日到事件发生或末次随访的月份数。只保存relapse一个字段是不够的，因为随访时间短的未复发病人，很可能是“还没来得及复发”而不是“安全”。建模阶段如果你的随访时间差异大，常见做法是设定一个固定时间窗，比如 12 个月内复发记为 1，未满 12 个月就失访的样本直接剔除，或用生存分析处理删失。后面避坑章会再强调这一点，但数据准备阶段就要按这个思路去整理，否则后期很被动。

3. 从逻辑回归到随机森林：建模主流程与参数选择

3.1 第一步先跑逻辑回归：这个基线决定了你的模型下限

临床预测模型领域有个共识：先做逻辑回归，把基线立住，再考虑要不要上机器学习。逻辑回归的优势不是预测精度，而是可解释性和稳定性。医生接受一个模型的理由往往是“我能看懂变量怎么影响结果”，而不是“准确率 95%”。

下面这段代码用清洗后的数据做训练集和测试集拆分，然后训练逻辑回归模型。

set.seed(2024) train_id <- sample(seq_len(nrow(df)), size = floor(0.7 * nrow(df))) train <- df[train_id, ] test <- df[-train_id, ] fit_glm <- glm(relapse ~ age_years + bmi + egfr + lvef + hypertension + diabetes + surgery_duration + cpb_time, family = binomial(link = "logit"), data = train) summary(fit_glm)

训练集和测试集按 7:3 拆分，set.seed 固定随机种子是为了结果可复现。glm函数的family参数指定二分类逻辑回归，link = "logit"是默认的 logit 连接。输出里重点看每个变量的回归系数、p 值和置信区间：p 值小于 0.05 的变量进入后续精简模型，但不要机械地只按 p 值筛特征——有些变量临床上重要但统计不显著，比如高血压史，建议同时参考系数方向和临床意义。

跑完summary后我一般做两件事：第一件是用vif()检查多重共线性，第二件是看 AIC 而不是看准确率。逻辑回归的预测结果是一个概率，把它转成 0/1 标签需要阈值，而这个阈值的选择会影响准确率，基线阶段先不纠结阈值，用概率值本身去评估模型。

3.2 升级到随机森林与 XGBoost：什么时候值得换算法

逻辑回归只能捕捉线性关系和少量交互项，如果变量之间存在复杂的非线性交互，比如“高龄 + 低 LVEF + 手术时间长”的组合效应，树模型就有优势。随机森林和 XGBoost 在医学预测里的定位是“更高上限但更低可解释性”。样本量几百例的医学数据集，树模型的优势并不一定发挥得出来，很多时候逻辑回归和随机森林的 AUC 差距不大。

library(randomForest) set.seed(2024) fit_rf <- randomForest( factor(relapse) ~ age_years + bmi + egfr + lvef + hypertension + diabetes + surgery_duration + cpb_time, data = train, ntree = 500, mtry = 3, importance = TRUE ) pred_rf <- predict(fit_rf, newdata = test, type = "prob")[, 2]

ntree设 500，够用的同时不会太慢。mtry是每个节点随机抽样的变量个数，分类问题的默认值是特征数的平方根，如果总共有 8 个变量，mtry = 3 左右即可。importance = TRUE会计算每个变量的重要性分数，后面可以用来做特征筛选的可视化。注意randomForest对缺失值不友好，数据里不能有 NA，这也是前面清洗时强调过滤缺失行的原因。

如果你的数据量到了几千例，或者逻辑回归明显不够用，再考虑 XGBoost：

library(xgboost) train_x <- train %>% select(age_years, bmi, egfr, lvef, hypertension, diabetes, surgery_duration, cpb_time) dtrain <- xgb.DMatrix(data = as.matrix(train_x), label = train$relapse) params <- list( objective = "binary:logistic", eta = 0.05, max_depth = 3, subsample = 0.8, colsample_bytree = 0.8, eval_metric = "auc" ) fit_xgb <- xgb.train(params, dtrain, nrounds = 300, verbosity = 0)

XGBoost 的参数比随机森林敏感，eta是学习率，设大了容易过拟合，设小了训练慢，0.05 是起步值。max_depth控制树深度，医学小样本数据深度到 3 到 5 就够了，再深就开始记住噪声。subsample和colsample_bytree是行采样和列采样比例，都在 0.8 附近是常规操作。nrounds是树的数量，这个值配合早停法使用更好。

3.3 模型对比不是只看准确率，要建一张统一评估表

很多教程到这里直接把随机森林和逻辑回归的准确率一比就结论了，这在小样本、类别不平衡的医学数据上很容易误判。复发事件通常发生率在 10% 到 30%，如果一个模型把所有样本都预测为“不复发”，准确率也有 70% 到 90%，看着很高，实际上没有任何临床价值。

正确的做法是同时记录 AUC、灵敏度、特异度、阳性预测值、阴性预测值，并且固定同一个阈值去比较。下面这段代码把三个模型放在同一套评估框架里。

library(pROC) pred_glm <- predict(fit_glm, newdata = test, type = "response") pred_rf <- predict(fit_rf, newdata = test, type = "prob")[, 2] roc_glm <- roc(test$relapse, pred_glm) roc_rf <- roc(test$relapse, pred_rf) # 以 Youden index 为阈值切分 thr_glm <- coords(roc_glm, "best", ret = "threshold")$threshold thr_rf <- coords(roc_rf, "best", ret = "threshold")$threshold glm_tbl <- data.frame( model = "glm", threshold = thr_glm, auc = as.numeric(auc(roc_glm)), sensitivity = coords(roc_glm, thr_glm, ret = "sensitivity")$sensitivity, specificity = coords(roc_glm, thr_glm, ret = "specificity")$specificity ) rf_tbl <- data.frame( model = "randomForest", threshold = thr_rf, auc = as.numeric(auc(roc_rf)), sensitivity = coords(roc_rf, thr_rf, ret = "sensitivity")$sensitivity, specificity = coords(roc_rf, thr_rf, ret = "specificity")$specificity ) bind_rows(glm_tbl, rf_tbl)

coords(roc, "best", ret = "threshold")是 pROC 包自动选择最优阈值的方式，默认用 Youden index，即灵敏度加特异度减一的最大值。这里有个细节：“最优阈值”只是统计意义上的最优，临床上可能更看重灵敏度，宁可误报也不漏报，那阈值就应该往低调。对比表的意义在于让你看清模型之间真实的差距，而不是靠一个准确率拍板。

4. 评估一个复发预测模型：AUC、校准度和 Brier 分数

4.1 区分度：ROC 曲线与 AUC 的正确读法

区分度回答的问题是：随机抽一个复发病人和一个未复发病人，模型给前者更高风险概率的概率有多大。AUC 是 0.5 到 1 之间的数，0.5 等于瞎猜，0.7 到 0.8 是医学预测模型的常见水平，0.9 以上要么运气好要么有问题，先怀疑标签泄漏，别高兴太早。

画 ROC 曲线的代码很简单，难的是解读。

plot(roc_glm, col = "blue", main = "ROC Curve") plot(roc_rf, col = "red", add = TRUE) legend("bottomright", legend = c("GLM", "RandomForest"), col = c("blue", "red"), lwd = 2) roc.test(roc_glm, roc_rf) # DeLong 检验比较两个模型

roc.test是 pROC 包提供的 Delong 检验，用来判断两个 AUC 的差异是否显著。小样本数据里 AUC 差异 0.03 到 0.05 就很可能不显著，这时选逻辑回归更合理，因为简单、可解释、好发表。比较模型时也要报告 AUC 的置信区间，ci.auc(roc_glm)能得到 95% 区间，如果两个模型的置信区间大幅度重叠，那它们的区分能力本质上差不多。

4.2 校准度：模型说 70% 就真的是 70% 吗

AUC 高只能说明“排序排得好”，不能说明“概率值本身可信”。一个模型把所有预测概率都抬高到 0.8 到 0.9 之间，只要排序对，AUC 照样好看，但每个概率值都偏离真实风险。校准度就是检验这个的：把预测概率分箱，对比每个箱内的平均预测概率和实际复发率。

library(tidyverse) cal_tbl <- tibble( pred = pred_glm, actual = test$relapse ) %>% mutate(bin = cut(pred, breaks = seq(0, 1, by = 0.1), include.lowest = TRUE, right = FALSE)) %>% group_by(bin) %>% summarise( mean_pred = mean(pred), mean_actual = mean(actual), n = n() ) ggplot(cal_tbl, aes(x = mean_pred, y = mean_actual)) + geom_point(size = 3) + geom_abline(slope = 1, intercept = 0, linetype = "dashed") + coord_equal() + labs(x = "预测复发概率", y = "实际复发率")

按 0.1 的间隔分箱后，如果点大多数落在对角线上，说明校准良好。样本量小的时候分箱会不均匀，某些箱可能只有几个样本，这时分箱宽度可以放宽到 0.2，或者用rms::calibrate()做平滑校准曲线。校准不好常见两个原因：一是模型遗漏了重要的预测变量，二是训练集和测试集的疾病谱不同。

4.3 Brier 分数：区分度和校准度的“合订本”

Brier 分数是预测概率与实际结果之间均方误差，一个数字同时惩罚排序错误和概率偏差。公式很直观：对每个样本算(预测概率 - 实际结局)²，再取平均。实际结局是 0 或 1，所以一个完美模型的 Brier 是 0，AUC 很高但校准稀烂的模型，Brier 会暴露得很彻底。

brier_glm <- mean((test$relapse - pred_glm)^2) brier_rf <- mean((test$relapse - pred_rf)^2)

在无信息模型的基准值上校准理解更有用：如果全部预测为事件发生率 p，那么基准 Brier 是p * (1 - p)。你的模型 Brier 要明显低于这个基准才有增量价值。比如事件发生率为 0.2，基准 Brier 就是 0.16，模型 Brier 如果在 0.1 以下，算表现不错。

5. 复发预测最容易踩的 5 个坑：现象、原因、解决

5.1 坑一：术后变量混进特征集，出现标签泄漏

现象：模型 AUC 高达 0.95，科室同事觉得不可思议，推敲后也觉得数据没问题，但换一批病人就崩到 0.6。

原因：你把术后第七天的并发症、术后用药情况这些“发生后才能知道”的变量，当成了预测因子。复发和术后并发症共享同一个病因链，模型实际上是学到了“术后已经出问题的人更容易复发”。这是数据前置窗口的问题。

解决：建模前写一份变量清单，逐项标注“在预测时间原点是否可得”。只保留术前和术中变量，或者如果你明确预测的是“术后 30 天内的早期复发”，才允许纳入术后早期可测的指标。这个规则必须在特征筛选之前定死，而不是事后补救。

5.2 坑二：随访时间不足，未复发被当成安全

现象：模型的阴性预测值高得离谱，但临床反馈说“你们说安全的病人，半年后还是回来了”。

原因：随访 1 个月和随访 24 个月的病人都标记为“未复发”，前者很可能只是还没到复发时间。这种数据叫“删失”，直接用 logistic 回归处理会把未复发组污染成“安全样本”。

解决：最彻底的办法是改用生存分析，用 Cox 比例风险模型或随机生存森林，把time_to_event作为 Y 变量而不是简单的 0/1。退一步的做法是只保留随访时间达到结局定义时间窗的样本，比如研究 12 个月内复发，就只保留随访满 12 个月或有事件发生的病人。代价是样本量损失，但模型可信度大幅提升。

5.3 坑三：复发率太低，准确率和 AUC 双双失真

现象：复发率只有 8%，模型把所有病人预测为不复发，准确率 92%，AUC 也可能到 0.75 以上。

原因：类别不平衡让模型倾向预测多数类。准确率在这个场景毫无意义，AUC 也会因为多数类样本主导排序而产生乐观估计。

解决：评估指标换成灵敏度、特异度、阳性预测值、阴性预测值、Brier 分数。模型层面可以尝试对少数类加权，逻辑回归里设置weights参数，随机森林里设sampsize平衡两类样本。但要注意，做了平衡处理后的预测概率不再是真实风险，模型校准会受到影响。

5.4 坑四：连续变量直接进逻辑回归，线性假定被无视

现象：年龄和肌酐放进逻辑回归后系数不显著，但你用分组箱图看，年龄和复发率明显呈 U 型关系。

原因：逻辑回归假设连续变量与 logit 概率之间是线性关系。心功能指标往往不是线性影响复发风险的，LVEF 低是高风险，但 LVEF 从 55% 升到 70%，风险不一定继续下降。

解决：先对每个连续变量做非线性检查，最简单的办法是画loess平滑曲线看趋势，出现“倒 U”或“U”型就考虑用限制性立方样条（rms::rcs()）或直接转成有序分类变量。转分类的代价是损失连续信息，样本量够的时候样条是更好的选择。

5.5 坑五：缺失数据直接删行或均值填充，信息浪费还引入偏倚

现象：建模时为了省事删掉了所有带缺失值的样本，样本量从 800 掉到 300，模型 AUC 下降，而且新数据预测效果不稳定。

原因：不同变量的缺失模式可能和结局有关，比如病情重的患者往往检验项目做得更全，直接删除这部分样本等同于删除了高危人群。均值填充则压缩了变量方差，削弱了它和其他变量的相关结构。

解决：少于 5% 缺失的变量用中位数填充可以接受，但更稳妥的做法是多重插补，用mice包生成多个数据集，每个数据集各建一次模型，最后合并系数估计。如果缺失集中在少数变量且缺失率过高，优先考虑在模型中剔除变量而不是插补。这些选择都要在方法学部分写清楚，审稿人和临床医生都会追问。

6. 把模型用起来：列线图、决策曲线和外部验证

6.1 列线图：把回归系数变成医生床边能用的工具

逻辑回归模型的最大落地优势是能转成列线图（nomogram）：每个变量按取值对应一个分数，所有分数相加得到总分，再映射到复发概率。这个工具不需要电脑，一张纸就能在床边用。

library(rms) ddist <- datadist(train) options(datadist = "ddist") fit_lrm <- lrm(relapse ~ age_years + bmi + egfr + lvef + hypertension + diabetes + surgery_duration + cpb_time, data = train) nom <- nomogram(fit_lrm, fun = plogis, lp = FALSE, funlabel = "复发概率") plot(nom)

lrm函数和glm的区别是它基于最大似然估计且专为回归模型可视化设计，datadist是 rms 包画列线图必须的前置设置。fun = plogis告诉模型把线性预测值转成概率。列线图发布时最好附上受试者工作特征曲线和校准曲线的截图，这是预测模型类论文的标配三件套。

6.2 用决策曲线分析选阈值，而不是用 0.5 一刀切

所有预测模型最终都要面对同一个问题：预测概率超过多少算高危？有人习惯用 0.5，但这个值对复发率只有 10% 的疾病明显不合理。决策曲线分析（DCA）回答的是“在不同阈值下，使用这个模型对比不使用它，患者的净获益是多少”。

library(rmda) test$pred_glm <- pred_glm dc <- decision_curve(relapse ~ pred_glm, data = test, study.design = "case-control", policy = "optimum", fitted.risk = TRUE) plot_decision_curve(dc)

DCA 图的横轴是阈概率，纵轴是净获益率。模型曲线高于“全部治疗”和“全部不治疗”两条参考线时，说明在这个阈值下使用模型是合算的。你想要的临床场景是“宁可多干预也不漏掉复发”，就选曲线较高区间偏左侧的阈值；想控制过度治疗，选偏右侧的阈值。这个步骤虽然不直接影响模型训练，但它决定了模型如何真正进入决策流程。

6.3 外部验证：别让模型的成绩单只在一家医院有效

内部测试集做得再好，也回答不了一个问题：换个医院、换个地区的人群，模型还灵不灵。外部验证是预测模型证据链里最后一环，也是最常被跳过的一环。做法很简单：拿到另一家医院的数据，用已训练好的模型直接预测，再重新计算 AUC、校准曲线和 Brier 分数，看性能掉了多少。

ext <- read_csv("external_cohort.csv") ext$pred_external <- predict(fit_glm, newdata = ext, type = "response") roc_ext <- roc(ext$relapse, ext$pred_external) auc(roc_ext) brier_ext <- mean((ext$relapse - ext$pred_external)^2)

如果外部验证的 AUC 比内部验证明显下降，不要急着调参数补救，先检查变量定义的一致性——比如 LVEF 是超声测的还是核素测的，肌酐是入院时还是术前一天的，这些细微差别在单中心数据里不存在，跨中心一定会有。预测模型的终点不是模型训完，而是它能在你拿不到原始数据的地方依然站得住。这些年我自己的教训是：建模花一周，数据清洗花两周，外部验证花一个月。后面两步没有捷径可以省，希望帮到你。

本文还有配套的精品资源，点击获取