简介：这份资源是面向医学统计与临床研究方向的R语言实战教程，围绕心脏病术后复发预测这一具体课题，帮助具备基础R语法或统计背景的读者完成从数据到模型的完整分析链路。压缩包共906个文件，约41.52MB，以91个R脚本、20个Rmd报告、24个HTML文档和31张PNG图表为核心，辅以CSV、RDS、XLSX等数据文件及大量R包依赖目录，覆盖数据导入、清洗、建模与结果呈现各环节。教程内容涉及tidyverse数据整理、缺失值与不平衡数据处理、逻辑回归与随机森林、xgboost等模型对比，以及交叉验证、混淆矩阵、ROC曲线和AUC评估，并借助varImp等函数讨论特征重要性与模型可解释性。目前已有585人学习下载，适合希望把R语言应用于临床预测建模、需要可复现案例与排错参考的医学科研人员和数据分析学习者。
1. 心脏病术后复发预测：一份 R 语言医学分析教程能落地到什么程度
心脏病术后复发预测，听起来像只有三甲医院统计科才碰的课题，但真拿到数据你会发现，卡住大多数人的不是医学知识，而是 R 语言里那些琐碎到让人翻车的细节：变量类型不对、缺失值没处理干净、正负样本比例悬殊到模型直接摆烂。这份「R医学分析-心脏病术后复发预测教程」就是冲着这些具体问题来的，它把数据导入、清洗、建模、评估、解释串成一条完整链路，用的都是tidyverse、caret、randomForest、xgboost这些实战里真会装的包。适合谁？R 语言入门到中级、手头有医学随访数据、想把复发预测跑通并写进论文或报告的人。如果你还在纠结read.csv和read_csv的区别，或者glm()报错看不懂，这份教程能帮你把流程立起来。
2. 数据导入与清洗：从原始随访表到可建模数据框
2.1 为什么先啃readr和data.table而不是 base R
医学随访数据动辄几千行、上百列，用read.csv读进来慢不说，字符串列默认转 factor 的旧行为经常埋雷。教程里推荐readr::read_csv()，它读得快、不自动转 factor、列类型推断更保守，配合problems()能直接定位解析失败的行。如果数据量再大一点，比如几万行以上的多中心数据，data.table::fread()是更稳的选择，内存占用和速度都明显优于 base R。我一般会先fread探一下列数和分隔符，再用read_csv做正式导入，因为read_csv的col_types参数可以显式指定每列类型，避免后续建模时因子水平对不上。
library(readr) library(dplyr) # 显式指定列类型，避免自动推断把 ID 读成 numeric col_spec <- cols( patient_id = col_character(), age = col_double(), sex = col_factor(levels = c("Male", "Female")), hypertension = col_factor(levels = c("No", "Yes")), diabetes = col_factor(levels = c("No", "Yes")), cholesterol = col_double(), smoking = col_factor(levels = c("No", "Yes")), family_history = col_factor(levels = c("No", "Yes")), recurrence = col_factor(levels = c("No", "Yes")) ) df <- read_csv("heart_recurrence.csv", col_types = col_spec) glimpse(df)这段代码的关键在col_spec：patient_id必须是字符型，否则前导零会丢；二分类变量显式设成 factor 并指定 levels 顺序，后面confusionMatrix()才不会把正类搞反。glimpse()用来快速确认每列类型和缺失情况，比str()输出更紧凑。参数上，col_factor的levels顺序很重要，caret默认把第一个 level 当负类，如果你把 "Yes" 放前面，召回率算出来就是反的。
2.2 缺失值处理：is.na()只是起点，别急着na.omit()
教程里提到is.na()、complete.cases()、na.omit()，但直接na.omit()在医学数据里是血泪经验级别的坑——删着删着样本少一半，而且删掉的往往是有合并症的重症患者，模型学出来的结论直接偏掉。常见做法是先做缺失模式可视化，用naniar包看每个变量的缺失比例和共现关系，再决定是删除、插补还是保留。对于缺失比例低于 5% 的连续变量，中位数插补够用；分类变量可以单独设一个 "Unknown" 水平；缺失比例超过 20% 的变量，要么找临床意义解释，要么直接排除并在论文里说明。
library(naniar) library(mice) # 缺失模式可视化 vis_miss(df) # 连续变量中位数插补，分类变量加 Unknown 水平 df_imputed <- df %>% mutate( cholesterol = ifelse(is.na(cholesterol), median(cholesterol, na.rm = TRUE), cholesterol), smoking = forcats::fct_explicit_na(smoking, na_level = "Unknown") ) # 多变量插补示例（适合缺失较多的场景） imp <- mice(df, m = 5, method = "pmm", seed = 123) df_mice <- complete(imp, 1)vis_miss()输出的是缺失热图，能一眼看出缺失是随机的还是成片的。mice的m = 5表示生成 5 个插补数据集，method = "pmm"是预测均值匹配，适合连续变量。complete(imp, 1)取第一个插补集，正式分析时应该对每个插补集分别建模再合并结果，但教程阶段先用一个集跑通流程没问题。注意fct_explicit_na会把 NA 变成一个显式水平，后面建模时这个水平也会参与，如果临床认为 "Unknown" 没有意义，就在建模前用step_unknown或手动排除。
2.3 不平衡数据的处理边界
心脏病术后复发率通常不高，正负样本 1:5 甚至 1:10 很常见。教程提到过采样、欠采样和 SMOTE，这里要讲清楚边界：欠采样会丢信息，样本量本来就少的时候别用；过采样简单复制会让模型过拟合；SMOTE 在smotefamily或themis包里都有实现，但它对分类变量的处理需要额外注意，默认的 KNN 插值可能在因子变量上生成无意义的组合。我一般会先跑一版原始数据看基线，再用themis::step_smote()在recipe流程里做，这样交叉验证时重采样只在训练折内进行，避免数据泄漏。
library(themis) library(recipes) rec <- recipe(recurrence ~ ., data = df_imputed) %>% step_normalize(all_numeric()) %>% step_smote(recurrence, over_ratio = 0.5) # 在 caret 或 tidymodels 流程里应用over_ratio = 0.5表示少数类扩到多数类的 50%，不是 1:1，这样比全量过采样稳一些。step_normalize放在 SMOTE 前面，因为 SMOTE 基于距离计算，量纲不统一会出问题。这一步的坑在于：如果你在完整数据集上先 SMOTE 再交叉验证，验证集里混入了合成样本，AUC 会虚高，论文送审被质疑数据泄漏就麻烦了。
3. 建模与调参：glm、randomForest、xgboost怎么选怎么跑
3.1 逻辑回归作为基线：系数就是可解释性
医学领域逻辑回归仍然是首选基线，因为它的系数可以直接换算成 OR 值，临床医生看得懂。glm(recurrence ~ ., family = binomial())跑起来简单，但要注意：因子变量的参考水平、连续变量的线性假设、多重共线性。教程里用cor()看相关性，但cor()只对连续变量有效，分类变量要用vcd::assocstats()或caret::findCorrelation()。我一般会先跑单因素逻辑回归筛变量，再把 P 值小于 0.1 的放进多因素模型，最后用step()做逐步回归，虽然逐步回归有争议，但在变量数不多的时候作为探索够用。
library(broom) # 单因素筛选 univars <- names(df_imputed)[!names(df_imputed) %in% c("patient_id", "recurrence")] uni_results <- lapply(univars, function(v) { f <- as.formula(paste("recurrence ~", v)) tidy(glm(f, data = df_imputed, family = binomial())) }) names(uni_results) <- univars # 多因素模型 multi_model <- glm(recurrence ~ age + sex + hypertension + diabetes + cholesterol + smoking + family_history, data = df_imputed, family = binomial()) summary(multi_model) exp(coef(multi_model)) # OR 值tidy()把模型结果转成数据框，方便批量提取 P 值和系数。exp(coef())得到 OR 值，大于 1 表示风险增加，小于 1 表示保护因素。注意glm默认把 factor 的第一个 level 当参考，所以sexMale的系数是相对于 Female 的。如果某个变量在单因素里显著、多因素里不显著，别急着删，可能是共线性或交互作用，试试加交互项或做分层分析。
3.2 随机森林与 XGBoost：调参不是玄学，有顺序
随机森林和 XGBoost 在医学预测里用得多，但调参顺序错了就是浪费时间。randomForest先调mtry（每棵树分裂时随机选的变量数），再调ntree（树的数量），nodesize对不平衡数据影响大，可以适当调小。XGBoost 参数多，但核心就几个：max_depth、eta、nrounds、subsample、colsample_bytree。教程里用caret::train()做交叉验证，trainControl设method = "cv"、number = 5、classProbs = TRUE、summaryFunction = twoClassSummary，这样输出的是 AUC 而不是准确率，对不平衡数据更合理。
library(caret) library(randomForest) library(xgboost) ctrl <- trainControl( method = "cv", number = 5, classProbs = TRUE, summaryFunction = twoClassSummary, savePredictions = "final" ) # 随机森林 rf_grid <- expand.grid(mtry = c(2, 3, 4, 5)) rf_model <- train(recurrence ~ ., data = df_imputed, method = "rf", trControl = ctrl, tuneGrid = rf_grid, metric = "ROC", ntree = 500) # XGBoost xgb_grid <- expand.grid( nrounds = c(100, 200), max_depth = c(3, 5), eta = c(0.05, 0.1), gamma = 0, colsample_bytree = 0.8, min_child_weight = 1, subsample = 0.8 ) xgb_model <- train(recurrence ~ ., data = df_imputed, method = "xgbTree", trControl = ctrl, tuneGrid = xgb_grid, metric = "ROC")twoClassSummary会输出 ROC、敏感度、特异度，metric = "ROC"让调参以 AUC 为目标。savePredictions = "final"保存最终模型的预测结果，后面画 ROC 曲线直接用。随机森林的mtry一般从sqrt(变量数)附近试，XGBoost 的eta小一点、nrounds大一点通常更稳，但训练时间会拉长。注意train()默认会做预处理，如果已经手动标准化过，加preProcess = NULL避免重复。
3.3 交叉验证的坑：数据泄漏和分层
交叉验证最怕数据泄漏。比如你先在完整数据上做了 SMOTE，再交给train()做 CV，验证折里混了合成样本，AUC 虚高。正确做法是把重采样嵌进recipe或trainControl的sampling参数里。另外，createFolds()默认随机分折，不平衡数据要加list = FALSE并手动分层，或者直接用trainControl的classProbs配合twoClassSummary，caret会自动做分层。如果样本量很小（比如少于 200），5 折 CV 每折验证集才 40 个样本，AUC 波动会很大，这时候用留一法或重复 5 折更稳，但计算量成倍增加。
# 分层交叉验证 set.seed(123) folds <- createFolds(df_imputed$recurrence, k = 5, list = TRUE, returnTrain = FALSE) ctrl_strat <- trainControl( method = "cv", number = 5, index = lapply(folds, function(x) setdiff(1:nrow(df_imputed), x)), classProbs = TRUE, summaryFunction = twoClassSummary )index参数手动指定每折的训练集索引，setdiff把验证折排除掉，这样保证每折的正负比例和整体一致。如果不想手动写，caret在method = "cv"时对 factor 类型的因变量默认就是分层抽样，但显式写出来更放心。
4. 模型评估与解释：混淆矩阵、ROC 和特征重要性
4.1confusionMatrix()的正类方向别搞反
caret::confusionMatrix()输出一堆指标，但正类方向搞反是新手最常见的翻车点。confusionMatrix(data = pred, reference = truth, positive = "Yes")里的positive参数必须显式指定，否则 R 按字母顺序把 "No" 当正类，敏感度和特异度直接对调。教程里提到准确率、精确率、召回率、F1，这些在confusionMatrix的输出里都有，但医学场景下召回率（敏感度）通常比精确率重要，因为漏诊一个复发患者的代价比误诊大。AUC 用pROC包算，roc(response, predictor)然后auc()，多模型比较用roc.test()做 DeLong 检验。
library(pROC) # 假设 pred_prob 是模型输出的正类概率 roc_obj <- roc(df_imputed$recurrence, pred_prob, levels = c("No", "Yes")) plot(roc_obj, print.auc = TRUE, main = "ROC Curve") auc(roc_obj) # 混淆矩阵 pred_class <- ifelse(pred_prob > 0.5, "Yes", "No") confusionMatrix(factor(pred_class, levels = c("No", "Yes")), df_imputed$recurrence, positive = "Yes")roc()的levels参数指定负类在前、正类在后，print.auc = TRUE直接在图上标 AUC。阈值 0.5 不是固定的，医学场景可以根据约登指数找最佳截断点，pROC::coords(roc_obj, "best", best.method = "youden")能直接算出来。注意confusionMatrix的data和reference都必须是 factor 且 levels 一致，否则报错。
4.2 特征重要性：varImp()和 SHAP 的取舍
随机森林和 XGBoost 的特征重要性用caret::varImp()就能出，但它给的是基于节点不纯度或增益的排名，对分类变量水平多的会有偏。更稳的做法是用iml或shapviz包算 SHAP 值，能看出每个特征对单个预测的贡献方向。教程里提到varImp()，作为快速筛选够用，但如果要写进论文，SHAP 图比重要性条形图更有说服力。逻辑回归直接看系数和 OR 值，配合forestmodel包画森林图，临床医生接受度最高。
library(iml) # 用 iml 算 SHAP 值 predictor <- Predictor$new(rf_model, data = df_imputed[, -1], y = df_imputed$recurrence) shap <- Shapley$new(predictor, x.interest = df_imputed[1, -1]) plot(shap) # 逻辑回归森林图 library(forestmodel) forest_model(multi_model)Predictor$new()包装模型和数据，Shapley$new()对单个样本算 SHAP，plot()出图。forest_model()直接接受glm对象，输出 OR 值和置信区间的森林图，比手动ggplot省事。注意 SHAP 计算量随特征数指数增长，特征多的时候用shapviz的采样版本或fastshap近似。
5. 避坑与排查：五个真实翻车记录
5.1 现象：confusionMatrix报错 "levels of data and reference must be the same"
原因：模型预测的 factor levels 和真实标签的 levels 顺序或内容不一致，常见于predict()输出丢了 levels 或 SMOTE 后新生成的样本标签类型变了。解决：统一用factor(x, levels = c("No", "Yes"))强制转换，或者在trainControl里设classProbs = TRUE后直接用概率阈值分类，不走predict()的默认分类。
5.2 现象：随机森林varImp显示某个变量重要性为负
原因：randomForest包对分类变量的重要性计算在某些情况下会出负值，尤其是变量水平多但样本少的时候。解决：换ranger包跑随机森林，或者用party::cforest()的条件推断树，重要性更稳。如果只是筛选变量，负值直接当 0 处理，别硬解释。
5.3 现象：XGBoost 训练集 AUC 0.99，验证集 AUC 0.6
原因：典型过拟合，max_depth太深或nrounds太多，模型把训练集噪声也学了。解决：先降max_depth到 3 或 4，加subsample = 0.7、colsample_bytree = 0.7，再用xgb.cv()看早停轮数，early_stopping_rounds = 20能自动停在验证误差最低点。
5.4 现象：mice插补后模型结果和完整案例差很多
原因：插补模型本身有偏，或者插补时用了包含结局变量的信息导致泄漏。解决：插补时把结局变量排除在预测矩阵外，用predictorMatrix参数设结局列全为 0。插补后做敏感性分析，比较完整案例、中位数插补、多重插补三种结果，如果差异大就在论文里讨论。
5.5 现象：step_smote报错 "All columns must be numeric"
原因：themis::step_smote()默认要求所有预测变量是 numeric，因子变量需要先step_dummy()转成哑变量。解决：在recipe里先step_dummy(all_nominal(), -all_outcomes())，再step_smote()。注意哑变量化后变量数增加，SMOTE 的 KNN 距离计算会变慢，样本大时考虑step_adasyn()或step_rose()替代。
6. 把模型塞进临床工作流：一个可复现的预测脚本长什么样
跑通单个模型只是起点，真正落地要的是一个能重复执行的脚本：读新数据、套用训练好的预处理参数、输出每个患者的复发概率和风险分层。我一般会把recipe和模型一起存成.rds，新数据来了直接bake()和predict()，避免重新拟合预处理步骤导致结果对不上。下面这个脚本模板我用了很多次，改改路径和变量名就能套。
library(tidymodels) library(readr) # 保存训练好的 workflow final_wf <- workflow() %>% add_recipe(rec) %>% add_model(linear_reg() %>% set_engine("glm")) # 这里以逻辑回归为例 # 假设 final_wf 已经 fit 过 # saveRDS(final_wf, "heart_recurrence_model.rds") # 新数据预测 new_data <- read_csv("new_patients.csv") model <- readRDS("heart_recurrence_model.rds") preds <- predict(model, new_data, type = "prob") %>% bind_cols(new_data %>% select(patient_id)) %>% mutate( risk_level = case_when( .pred_Yes < 0.2 ~ "Low", .pred_Yes < 0.5 ~ "Medium", TRUE ~ "High" ) ) write_csv(preds, "recurrence_predictions.csv")这个脚本的关键在predict(model, new_data, type = "prob")，type = "prob"输出正类概率而不是硬分类，方便做风险分层。case_when的阈值 0.2 和 0.5 是根据约登指数和临床共识定的，不同数据集要重新校准。write_csv输出结果给临床端，patient_id保留用于回溯。注意新数据的列名和类型必须和训练时一致，read_csv的col_types最好也固定下来，否则因子水平对不上会报错。
验证模型稳定性我习惯做两件事：一是用bootstrap重采样 1000 次算 AUC 的置信区间，二是换一个时间段的随访数据做外部验证。内部验证 AUC 0.85、外部验证掉到 0.7 是常态，别慌，说明模型有过拟合，回头检查变量筛选和调参过程。如果外部验证 AUC 低于 0.65，基本要考虑重新设计特征或换模型了。
从那以后我每次跑医学预测模型，都强制走一遍「完整案例 → 多重插补 → SMOTE → 交叉验证 → 外部验证」的流程，少一步都不敢往论文里写。希望帮到你。
本文还有配套的精品资源，点击获取
