临床预测模型R语言全流程：从Lasso筛选到列线图与DCA实战-北京尧图网络

1. 为什么你跑不通那些“公开代码”先搞清预测模型的全貌再动手经常有临床的朋友带着同一类问题来找我“师兄我从网上下了一份列线图的R代码数据格式也对好了可一运行就是报错卡在glm那一步到底哪错了”问得多了我就发现问题几乎都不是出在最后那几行画图代码上而是出在一个大家普遍忽略的前提上——临床预测模型的R语言代码从来不是“复制粘贴就能跑”的东西。它背后是一整条流水线数据清洗、变量筛选、模型训练、内外部验证、可视化呈现。你下载的那份代码只是流水线末端的“画图车间”前面的数据如果没按它预期的格式准备好后面自然全线崩溃。这篇帖子我想聊的不是某一份具体代码的讲解而是把这条流水线完整地捋一遍从一堆原始病历数据出发到最终产出森林图、列线图、校准曲线和DCA曲线每一步的R语言实战代码怎么写、为什么这么写、哪些地方容易踩坑。我会尽量用一套可以跑通的最小示例来串联涉及glm逻辑回归建模、lasso筛选变量、用rms包做列线图、pROC和ggplot2做ROC与校准曲线以及rmda或自写函数做DCA。末尾再讲几个我实际项目中踩过、也帮别人调过的经典报错。这篇内容适合谁两类人。一类是刚开始做临床预测模型的硕博生和临床医生你的数据已经有了但面对RStudio不知道先敲哪行命令另一类是已经能跑通基础代码但想知道“为什么这里要这样处理变量”“为什么列线图分数跟logit变换对不上”的进阶用户。如果你是后者可以直接跳到第3章以后的建模和验证部分。2. 数据预处理拿到原始病历后前两个小时的活全在这里这章节负责从入院记录、检验报告等原始数据转换成模型能够直接使用的结构化分析集。很多时候这部分被当成“体力活”但临床预测模型的结论可靠性恰恰取决于此。2.1 数据导入的第一步统一变量名与数据格式R里最痛苦的从来不是统计建模而是数据进来的时候各列的名字和你代码里写的不一致。医院导出的数据常见格式是Excel或CSV但列名可能是中文、带空格甚至同一个指标在不同批次的导出里叫法都不一样。我建议的第一步不是马上跑任何分析而是先做一次变量名的标准化映射。用一个简单的示例来说明假设你要做的是“急性胰腺炎患者发生重症的预测模型”原始数据里可能有这么几列library(readxl) rawdata - read_excel(your_raw_data.xlsx, sheet 1) # 看一眼真实列名 names(rawdata)通常你会看到类似患者编号、性别、年龄岁、入院24h内血钙这样带单位的列名。在R里带中文、带括号的列名虽然能用反引号访问但后续写公式时极其容易出错。所以第一步是清洗列名并做类型转换# 统一列名变成模型数据的标准格式 mydata - rawdata colnames(mydata) - c(id, age, sex, bmi, calcium, crp, apache2, outcome) # 结局变量转成因子0非重症1重症 mydata$outcome - factor(mydata$outcome, levels c(0, 1), labels c(no, yes))这里的核心习惯是先建一张“变量字典”把原始字段和模型变量名一一对应并确认每个变量的类型是数值型还是因子型。因子型变量在后续glm中会自动生成哑变量如果忘记转因子R会把它当连续变量代入这里的疏漏容易造成结果的明显偏差。这个过程虽然琐碎但如果没有在数据读取阶段就统一好后续每画一张图可能都要回头修。2.2 缺失值的处理不能无脑删除临床数据里缺失值几乎是必然存在的。很多教程会告诉你“用na.omit()删掉所有含缺失的行”这在变量少、缺失少的时候勉强可行但实际中如果直接删除常常会损掉大量有效样本——而且数据缺失机制往往不是完全随机的删掉后会导致选择偏倚。我推荐的做法是两步。第一步先看一下每个变量的缺失比例# 查看缺失率 missing_pct - sapply(mydata, function(x) sum(is.na(x)) / nrow(mydata)) sort(missing_pct, decreasing TRUE)缺失率在5%以下的连续变量可以直接用中位数或均值填补影响很小。缺失率在5%~20%的变量建议用mice包做多重插补尤其当它是某个临床评分如APACHE II评分的分项时。缺失率超过40%的变量除非有极强的临床依据否则建议直接在建模阶段剔除。第二步对于关键变量明确区分“真缺失”和“临床未测”。比如血钙值缺失可能是因为医生判断没有必要测这种“未测”本身可能隐含疾病严重度信息。可以通过增加一个“是否检测”的指示变量来捕捉这种信息# 为钙离子列增加一个缺失指示变量 mydata$calcium_miss - ifelse(is.na(mydata$calcium), 1, 0) mydata$calcium[is.na(mydata$calcium)] - median(mydata$calcium, na.rm TRUE)这个操作在很多预测建模实践中都被证明能小幅提升模型表现因为“医生决定测不测某项指标”往往和病情有关。这不是让你把所有缺失都这样处理而是提醒你缺失值处理本身就是建模策略的一部分不应该无脑交给na.omit()。2.3 连续变量要不要分箱临床预测模型里经常涉及连续变量分箱。比如年龄有些人喜欢分成“≤60”和“60”或按四分位数分组。R里的切分很简单library(dplyr) mydata - mydata %% mutate(age_group cut(age, breaks c(0, 60, 75, 120), right FALSE, labels c(60, 61-75, 75)))但这里我要多说一句自己的体会。列线图Nomogram里连续变量并不一定要分箱。分箱会损失信息量而且分箱点选取不当会让模型拟合产生较大偏差。如果你做的预测模型最终目标是临床快速使用分箱能提升易用性但如果追求预测精度建议先用限制性立方样条RCS在rms包里查看非线性趋势再决定是原样进入模型还是分箱。怎么快速判断变量是否线性可以用rcs()做一下试探性建模library(rms) dd - datadist(mydata) options(datadist dd) fit_rcs - lrm(outcome ~ rcs(age, 4) calcium apache2, data mydata) anova(fit_rcs)看age的非线性检验P值。如果P值不显著说明把它当线性项处理没有大问题如果显著再考虑分箱或保留样条项。这一步虽然多花了五分钟但能让你在论文审稿时少一个被质疑的点审稿人通常会问“连续变量分箱的截断值依据是什么”如果你能说“我先用RCS检验了非线性发现阈值大约在62岁所以按60/75分箱”这就比凭空切分有说服力得多。3. 变量筛选与核心建模你为什么需要lasso又为什么不能只靠lasso数据准备好了下面是重头戏筛选变量并拟合预测模型。临床预测模型的主流模型仍然是逻辑回归。虽然机器学习模型随机森林、XGBoost在很多数据集上精度更高但临床场景最看重的往往是可解释性——你需要知道每个变量权重对应的临床意义所以逻辑回归加列线图依然是论文和实际临床工具的主流形态。3.1 先跑一个全变量模型看看基线表现不要一上来就做筛选。先把你临床上认为相关的变量全部放进模型得到一个“全模型”作为后续所有筛选策略的参照。这个模型的AUC、校准表现就是你后续比较的基线。# 全变量逻辑回归 full_model - glm(outcome ~ age sex bmi calcium crp apache2, data mydata, family binomial()) summary(full_model)看结果的时候关注两件事。第一样本量是否够。经验法则是模型中每个参数至少需要10~15个阳性事件。如果你的结局是重症患者人数共80例那模型最多容纳的预测变量数大约就是5~8个80/108。如果全模型放了15个变量那大概率过拟合后面需要大力筛选。第二变量的系数方向和临床认知是否一致。如果血钙的系数显著为正即血钙越高重症风险越高和临床直觉矛盾要回去检查是不是编码反了或者存在共线性。3.2 用lasso筛变量为什么它是高维变量的默认选项很多临床预测模型教程会直接建议用逐步回归stepwise来筛选变量然后在括号里提一句lasso。但近些年统计学家对逐步回归的态度已经比较谨慎——它基于P值反复进出变量会导致系数估计偏倚而且标准误偏小。相比之下lasso最小绝对收缩和选择算子通过给系数施加L1惩罚把不重要的变量系数压缩到0从而同时实现变量选择和系数收缩在高维相关变量存在时比逐步回归稳定得多。在R里做lasso最常用的包是glmnet。逻辑回归时设定family binomiallibrary(glmnet) # glmnet要求输入为矩阵且不能有缺失 x - as.matrix(mydata[, c(age, sex, bmi, calcium, crp, apache2)]) y - as.numeric(mydata$outcome yes) # 这里用alpha1代表lasso set.seed(2024) cv_fit - cv.glmnet(x, y, family binomial, alpha 1, nfolds 10) # 画出交叉验证误差曲线找到最优lambda plot(cv_fit)cv.glmnet会做10折交叉验证自动选出两个关键的lambda值lambda.min交叉验证误差最小的lambda和lambda.1se在最小误差一个标准误范围内的最大lambda模型更简洁。实际使用中我通常取lambda.1se对应的变量集因为它在精度损失很小的情况下变量数少一半临床使用更方便。# 提取最优lambda下的系数 coef_min - coef(cv_fit, s lambda.1se) # 查看哪些变量被保留 coef_min - as.matrix(coef_min) selected_vars - rownames(coef_min)[coef_min[, 1] ! 0] selected_vars一个小提醒glmnet要求输入数据完整所以要在前面缺失值填补完成后才能运行。如果你的数据里还有NA这里会直接报错。3.3 用lasso筛完再用临床判断兜底lasso选出来的变量集可以当作一个强有力的起点但不建议完全放弃临床判断直接照单全收。原因是lasso只在乎预测精度它可能选入某个实验室指标但对临床来说采集成本高或漏掉某个已经被大量文献证明的关键变量因为与另一个变量存在强相关导致它在惩罚下被压缩为0。我的一般操作是把lasso选出的变量作为候选集再和临床团队过一遍如果某个在临床上公认重要的变量没有被选入我会手动保留它然后用选定的变量重新拟合逻辑回归模型得到不带惩罚项的普通glm用于后续列线图制作。因为列线图本质上展示的是原始逻辑回归系数如果你直接用glmnet的系数做列线图会涉及尺度变换的问题比较麻烦。# 用lasso筛出的变量建立传统逻辑回归 lasso_selected - c(age, calcium, apache2) # 示例 final_model - glm(outcome ~ age calcium apache2, data mydata, family binomial()) summary(final_model)注意这一步不要再用逐步回归去筛了。lasso已经起到筛选作用你这一步只是把模型形式换成无惩罚的glm以便后续画列线图和计算更常规的统计量。如果再跑一次stepwise相当于二次筛选增加了过拟合风险。3.4 训练集与验证集划分不该忽略set.seed如果你想用一个相对独立的数据集来评估模型表现需要把数据拆成训练集和测试集。R里的标准操作是用caret或rsample包做分层抽样library(caret) set.seed(123) train_index - createDataPartition(mydata$outcome, p 0.7, list FALSE) train_data - mydata[train_index, ] test_data - mydata[-train_index, ]这里有两个关键点。第一为什么一定要set.seed()因为数据划分是随机的如果不设种子每次运行代码得到的训练集都不一样模型结果就无法复现。投稿时审稿人会要求你提供可复现的完整代码种子就是复现的一部分。第二为什么要分层抽样createDataPartition而不是简单的随机抽样sample()当结局比例较低时比如重症比例只有15%简单随机抽样很容易让测试集里一个阳性都没有或训练集里阳性占比过低。分层抽样保证划分后每一组里阳性比例与总体大致一致。3.5 临床预测模型中的样本量估算一个常被忽略的前提写完模型后建议顺手算一下你当前样本量是否支持这个模型。有一个便捷的做法是用pmsampsize包library(pmsampsize) pmsampsize(type b, # binary 结局 rsq 0.3, # 预估模型R方拿不准就填0.2左右 parameters 6, # 预计纳入的候选变量数 prevalence 0.15) # 结局事件比例它会输出两个重要的量所需的最小样本量以及所需的最小阳性事件数。如果你的实际样本量远低于这个值那就得考虑减少预测变量数量或者明确说明这是探索性分析。好多人吭哧吭哧建模画图最后审稿人说“样本量不足”再从第一步重来那是相当难受的。这个包的使用只需要一行代码能帮你提前发现这种结构性硬伤。4. 模型好不好看全看可视化森林图、列线图、校准曲线、DCA一网打尽临床预测模型的展示环节几乎决定了一篇论文给人的第一印象。很多人的模型统计量做得很扎实但图做得不够直观结果论文说服力打折。这一节我把四个最常用的图逐个讲清楚提供可以直接改参数运行的代码。4.1 森林图用forestploter展示单因素与多因素结果森林图的场景有两个一是展示单因素分析中每个变量和结局的关系未调整的OR值二是展示多因素模型中变量的调整后OR值。很多论文会两张图并排左边是单因素结果、右边是多因素结果。forestploter是近几年维护比较活跃的做森林图的包语法比forestplot更直观。先构造一个数据集包含变量名、OR、置信区间、P值这些列然后直接绘图library(forestploter) # 组装一个示例数据框 plot_data - data.frame( Variable c(Age, Sex (Male vs Female), BMI, Calcium, APACHE II), Label c(Age, Sex, BMI, Calcium, APACHE II), OR c(1.03, 1.42, 1.12, 0.82, 1.18), Low c(1.01, 0.88, 0.98, 0.71, 1.10), High c(1.05, 2.28, 1.28, 0.95, 1.27), P c(0.001, 0.145, 0.085, 0.012, 0.001) ) # 在变量名前创建一列空白占位以便放森林图CI区间 plot_data$ - paste(rep( , nrow(plot_data)), collapse ) # 指定绘图 forest_plot - forest(plot_data[, c(Variable, , OR, P)], est plot_data$OR, lower plot_data$Low, upper plot_data$High, ci_column 2, # 森林图放在第2列 ref_line 1, xlab OR (95% CI), ticks_at c(0.5, 1, 1.5, 2))有一点要注意ci_column指定的是包含空格的占位列这一列在视觉上只是给森林图留出空间。如果在实际项目中想显示OR的数值列需要把OR列放在est参数里同时再复制一列用于文本展示。具体的列布局可以自己微调但核心逻辑就是文本列空格列数据标注列。如果你用的是普通的forestplot包语法会不一样所以我通常建议锁定一个包用熟。4.2 列线图rms包里的经典操作列线图是临床预测模型中最经典的可视化输出。它把复杂的逻辑回归公式转换成一个带刻度的计算图读者按变量取值向上去线、对应顶部得到总分总分向下对应预测概率。R里做列线图一般用rms包而不是基础glm因为rms提供了一系列配套函数lrm、nomogram、calibrate、val.surv等专门为回归建模和可视化设计。建模代码library(rms) # 先声明数据分布 dd - datadist(mydata) options(datadist dd) # 用lrm拟合逻辑回归等效于glm binomial nom_model - lrm(outcome ~ age calcium apache2, data mydata, x TRUE, y TRUE) # 绘制列线图 nom_plot - nomogram(nom_model, fun plogis, fun.at c(0.05, 0.2, 0.5, 0.8, 0.95), lp FALSE, funlabel Risk of Severe AP) plot(nom_plot)这里有一个容易被忽视的细节lrm拟合时一定要设x TRUE, y TRUE否则后续做calibrate校准曲线时会报错“does not have x or y”。如果你已经用普通glm拟合好了模型想不重跑直接用rms画列线图也不是不行但很麻烦——nomogram要求模型对象是lrm。所以我一般直接从建模就切换到rms体系。还有一点fun plogis的含义是把线性预测值logit转换为概率。列线图最底部那一行的数值就是你的结局预测概率。这个转换和lp FALSE配合表示直接显示概率刻度而不是线性预测值。4.3 校准曲线预测概率和实际发生率是否一致校准曲线评价的是模型预测的概率与实际观测到的结局频率是否一致。通俗讲如果模型预测某一组病人重症概率是30%那这群人里实际发生重症的比例是否也在30%左右。校准曲线如果贴近对角线说明校准良好。R里两种常见做法。第一种是用rms自带的calibrate基于Bootstrap重采样cal - calibrate(nom_model, method boot, B 1000) plot(cal)第二种是ggplot2手工绘制把预测概率分箱比如按十分位然后计算每个分箱内的实际阳性率再画散点或折线。很多论文要求展示的是第二种。代码思路如下library(ggplot2) # 计算每个样本的预测概率 mydata$pred_prob - predict(nom_model, type fitted) # 按预测概率十分位分组 mydata$pred_group - cut(mydata$pred_prob, breaks quantile(mydata$pred_prob, probs seq(0, 1, 0.1)), include.lowest TRUE) # 计算每组的预测均值与实际发生率 cal_data - mydata %% group_by(pred_group) %% summarise(pred_mean mean(pred_prob), obs_rate mean(outcome yes), n n()) ggplot(cal_data, aes(x pred_mean, y obs_rate)) geom_point(size 2) geom_line() geom_abline(slope 1, intercept 0, linetype dashed) labs(x Predicted Probability, y Observed Rate) coord_equal()这个代码生成的校准图比较直观但不少审稿人更信任calibrate自带的Bootstrap曲线因为它在图中自带了一条平滑的Bootstrap校正曲线能直观看到模型的过拟合程度。我的做法是两种都生成论文正文用calibrate的版本补充材料放分箱校准图。4.4 受试者工作特征曲线ROC与曲线下面积AUCROC曲线用于评价模型的区分度即患病的和没患病的模型能不能分清。pROC包是做这个的标准选择。library(pROC) # 预测概率 pred_prob_test - predict(final_model, newdata test_data, type response) # 计算ROC roc_curve - roc(test_data$outcome, pred_prob_test, levels c(no, yes), direction ) auc_value - auc(roc_curve) ci_auc - ci.auc(roc_curve) # 输出AUC的95%置信区间 # 画图 plot(roc_curve, col #E64B35, lwd 2, main ROC Curve of Severe AP Prediction Model) legend(bottomright, legend paste0(AUC , round(auc_value, 3), (95% CI: , round(ci_auc[1], 3), -, round(ci_auc[3], 3), )), bty n)这里有个容易踩的坑auc()会默认按字母顺序把no当成阳性组而no是阴性所以计算结果方向是反的AUC会变成1-AUC。解决方式就是在roc()里明确指定levels c(no, yes)并设置direction 告诉R“高水平代表yes组”。这个方向问题如果不调整出来的AUC可能小于0.5很多人以为模型极差其实是方向反了。如果需要同时比较多个模型的AUC是否有显著差异用pROC::roc.test()roc_model1 - roc(test_data$outcome, pred_prob_model1, levels c(no, yes)) roc_model2 - roc(test_data$outcome, pred_prob_model2, levels c(no, yes)) roc.test(roc_model1, roc_model2, method delong)DeLong检验是临床论文里比较两个AUC最常用的方法P值小于0.05说明两个模型区分度有统计学差异。4.5 决策曲线分析DCA判断模型是否有临床净获益DCA一度是临床预测模型论文的“加分项”现在几乎是“必选项”了。它通过权衡“真阳性带来的获益”和“假阳性带来的伤害”展示在不同阈值概率下使用这个模型相对于“所有人都治疗”和“所有人都不治疗”两种策略的净获益。R里最方便的是rmda包library(rmda) # 用原始数据构建DCA dca_model - decision_curve(outcome ~ age calcium apache2, data mydata, family binomial(link logit), thresholds seq(0, 0.9, by 0.01), confidence.level 0.95) plot_decision_curve(dca_model, curve.names Full model, cost.benefit.axis TRUE, standardize FALSE, col #E64B35)thresholds参数指的就是医生决定采取干预时所用的概率阈值比如“预测概率超过20%就按重症处置”20%就是阈值概率。DCA会在这个阈值的全区间上计算净获益。如果模型曲线在大部分阈值区间都高于“Treat All”和“Treat None”两条参考线说明模型有临床应用价值。有人问DCA要不要在测试集上做。稳妥的做法是在测试集上做或者在全数据集上Bootstrap来做。rmda默认支持Bootstrap置信区间但跑起来稍慢。如果数据量不大直接用全数据DCA也能说明问题。5. 模型评价C指数、NRI、IDI以及那三个“反向暴击”的坑做完可视化不等于工作结束。一个可靠的临床预测模型必须经过严谨区分度、校准度和临床实用性三层评价外加内部验证。这一节我不只讲怎么算还会专门列出我自己踩过、也替别人调过的三个典型坑。5.1 C指数与AUC同一件事的两种叫法逻辑回归模型的C指数concordance index在二分类结局中和AUC是等价的都表示随机抽取一对患者模型把阳性患者的风险排得比阴性患者高的概率。rms包里的lrm对象会直接输出C指数nom_model输出结果的C那一行就是对训练集样本的C指数等价于训练集AUC。但如果你用glm建模可以用somers2()函数来计算library(Hmisc) somers2(predict(final_model, type response), as.numeric(mydata$outcome yes))somers2()返回的C即为C指数Dxy是Somers D两者关系是C 0.5 Dxy/2。5.2 内部验证为什么训练集AUC虚高Bootstrap告诉你真相很多新手只看训练集AUC觉得0.85就很好了。但模型在训练集上的表现通常虚高因为模型“记住”了一部分噪声。如果你在测试集上跌到0.75并不能完全说明模型差有时只是测试集样本量小、波动大。更稳的内部验证方式是Bootstrap重采样反复从全量数据中有放回地抽样每次拟合模型再验证计算乐观度并校正。rms里的validate函数能直接输出校正后的C指数set.seed(2024) val - validate(nom_model, method boot, B 1000) val输出结果中重点关注Dxy这一行的optimism和corrected列。corrected列对应的C指数就是校正后的C 0.5 corrected Dxy/2。如果校正后的C指数比训练集C指数低0.05以上说明过拟合程度可能会引起审稿人关注你需要在讨论里说明。5.3 新增变量到底有没有用NRI与IDI临床预测模型的进阶问题经常是“我在旧模型基础上加了一个新生物标志物它值不值得纳入”只看AUC提升往往不敏感——AUC可能只提升了0.01但临床净重分类改善却很明显。这时需要NRINet Reclassification Improvement和IDIIntegrated Discrimination Improvement。R里可以用nricens包比较嵌套模型library(nricens) # 标准方法基于样本事件/非事件两组的重分类表 nb - nricens(mdl.std list(glm(outcome ~ age calcium, data mydata, family binomial())), mdl.new list(glm(outcome ~ age calcium crp, data mydata, family binomial())), updown category, cut c(0.2, 0.4))cut是临床分类的阈值一般结合临床实践设定。如果你不确定阈值怎么选可以用updown diff这种基于连续预测值变化的方法不依赖分类阈值的NRI。做完这个之后描述结果时需要小心NRI的点估计范围可能很大不要只报一个值还要给置信区间。如果写论文我会把NRI的95%CI也一并报告方法是用Bootstrap抽样1000次重复估计NRI取2.5%和97.5%分位数作为置信区间。5.4 三个经典“反向暴击”坑这一节写三个我自己在真实项目中反复遇到过的问题全部是“代码能跑但结果对不上”的类型。第一个坑哑变量编码不一致。比如glm里性别被自动编码成0/1但当你用predict函数时如果新数据的因子水平顺序和训练数据不一致预测结果会完全错掉。解决办法是所有建模和预测都用同一个数据框切片不要在中间过程手动把因子转成数字。第二个坑lrm和glm的系数看起来差距很大。这通常不是bug而是两者对截距项或者编码方式的处理略有差异。lrm默认使用tol1e-12的奇异值处理在完全共线性变量存在时它能强行给出一个系数而glm会直接报错。如果在glm里遇到coefficients not defined because of singularities这类警告第一步不是删代码而是去查变量之间的相关性找出完全共线的那个列。第三个坑校准曲线在低风险人群里系统性偏移。很多临床数据里低风险人群占大头如果模型的预测概率整体偏低校准曲线前端会明显偏离对角线。这时不要第一时间认为是模型错了先检查你的事件率。比如数据里重症发生率是10%那模型预测的最高分段也可能只有30%—40%而你画图时xlim和ylim仍默认从0到1就会显得校准曲线“压扁”在对角线下方。解决办法是在两个坐标轴上统一用模型预测概率的实际范围。这个属于画图细节但在审稿人眼里就是校准不佳容易引起误解。6. 一份可运行的整套“最小工作流”代码从glm到列线图到外部验证前面每一章都是散点式的你可能需要一个能把所有环节串起来的“脚手架”。下面我给出一份尽量精简但能端到端跑通的最小工作流把前面的关键步骤压缩成一个流程。代码可以直接跑通拿自己的数据替换后基本就能出一套初稿图表。为了不无限拉长这里去掉了缺失填补的细节假设cleaned_data已经是一个无缺失、列名规范的数据框。library(rms) library(pROC) library(glmnet) library(rmda) library(forestploter) # 1. 准备RCS环境切分数据 set.seed(2024) cleaned_data - cleaned_data[sample(1:nrow(cleaned_data)), ] train_idx - 1:floor(0.7 * nrow(cleaned_data)) train_data - cleaned_data[train_idx, ] test_data - cleaned_data[-train_idx, ] # 2. 变量筛选lasso方案 x_train - as.matrix(train_data[, c(age, sex, bmi, calcium, crp, apache2)]) y_train - as.numeric(train_data$outcome yes) set.seed(123) cv_fit - cv.glmnet(x_train, y_train, family binomial, alpha 1) coef_lasso - coef(cv_fit, s lambda.1se) selected - rownames(coef_lasso)[which(coef_lasso[, 1] ! 0)] selected - setdiff(selected, (Intercept)) selected # 3. 用lrm拟合最终模型 dd - datadist(train_data) options(datadist dd) formula_text - paste(outcome ~, paste(selected, collapse )) final_lrm - lrm(as.formula(formula_text), data train_data, x TRUE, y TRUE) # 4. 列线图 nom - nomogram(final_lrm, fun plogis, fun.at c(0.05, 0.2, 0.5, 0.8, 0.95), funlabel Risk) plot(nom) # 5. 校准曲线Bootstrap内验证 cal - calibrate(final_lrm, method boot, B 500) plot(cal) # 6. 测试集ROC test_pred_prob - predict(final_lrm, newdata test_data, type fitted) roc_test - roc(test_data$outcome, test_pred_prob, levels c(no, yes), direction ) plot(roc_test) text(0.3, 0.3, paste0(AUC, round(auc(roc_test), 3))) # 7. DCA dca_res - decision_curve(as.formula(formula_text), data train_data, family binomial(link logit), thresholds seq(0, 0.9, by 0.01)) plot_decision_curve(dca_res, curve.names Model, standardize FALSE)这份最小工作流大约能应对70%的二分类型临床预测模型初稿需求。如果你做的是生存分析比如Cox回归做预后模型那就把lrm换成cph把nomogram里的fun换成Surv函数校准曲线用val.survROC要改用时间依赖ROCtimeROC包或survivalROC包框架整体一样但细节函数要替换。7. 写在外面的话代码能跑通只是起点最后说一点代码之外的东西。我自己帮人调过很多“预测模型”代码最终发现最花时间的不是报错本身而是使用者并不清楚这个报错在告诉他什么。比如glm报NA系数说明有共线性calibrate报错说明建模的时候忘了xTRUE, yTRUEDCA画出来曲线全部低于Treat All线大部分时候不是模型差而是阈值范围没设对或结局方向反了。R语言的报错信息确实不友好但如果你把每一段报错都当作数据告诉你的线索会慢慢练出一种排查直觉。我自己的习惯是每拿到新数据先跑一个很小的全流程建模ROC校准确认整体逻辑畅通然后再回去精修变量筛选、做插补、调图表细节。不要一上来就追求把列线图画得精美万分——如果前面的数据有误、模型有过拟合画得越漂亮后面返工越痛苦。希望这套从数据清洗、lasso筛选、逻辑回归建模到森林图、列线图、校准曲线、ROC和DCA的流程能帮你省下最开始的那些弯路。你如果在自己数据上跑出问题报错信息直接发过来我们可以接着往下排。

临床预测模型R语言全流程：从Lasso筛选到列线图与DCA实战

相关推荐

Vue大文件上传组件实战：切片上传、断点续传与Web Worker哈希计算

Element UI树形拖拽三级菜单拖不动？根因是children字段缺失

从电机控制到车规芯片平台开发：FOC、CAN与AUTOSAR的技术迁移路线

更多资讯

分析赚钱方式对人的影响与物化

搜索引擎算法更新应对指南：SEO型网站稳固自然流量的核心打法

事件驱动架构从理论到生产落地-架构师必知必会的EDA实战指南

书霸AI文献综述：把散乱资料变成研究地图

从分层架构到整洁架构-软件架构设计思想的演进与抉择

QT手写MQTT客户端：协议详解与工程实践指南

今日更新

能源系统对接碳交易的路径：从碳核算到碳资产运营

智能硬件数据驱动运营：订单日记与预测性维护实践

逆变器下垂控制与虚拟同步机对比：Simulink仿真建模与参数整定

本周精选

超人会飞不算本事：系统稳定依赖清晰规则与边界设计

超人VS蜘蛛侠：拆解超级IP的影响力与传播方法论

基于CNN的调制信号识别：MATLAB实现时频图分类实战

本月热门

基于YOLO的焊缝缺陷检测：从数据标注到工业部署全流程实践

分布式服务的决策记录

创业产品的首版范围

读完文章，想聊聊您的网站？