Cox构建预测模型（3）：如何用R语言做cox先单后多回归分析！（附全套代码）...

上一篇为大家介绍了预测模型中基线表格的绘制方式，那么在预测模型构建过程中，通过Cox回归来筛选预测因子也是十分关键的。在影响因素研究中，大家较为熟知的筛选方式是先单后多，但是在预测模型文章中，更为常用的筛选方法是逐步回归法。

逐步回归，是通过逐步将自变量输入模型，如果模型具统计学意义，并将其纳入在回归模型中。同时移出不具有统计学意义的变量。最终得到一个自动拟合的回归模型。说到这里大家可能会有点理解了，逐步回归法的目的是构建一个拟合佳的模型，该变量在回归过程中是否具有统计学意义无关，这与先单后多有着本质的区别。

本文将继续以本文将以复现一篇seer公共数据库文章（Q1 IF=5.2)为例，为大家介绍用R语言进行Cox回归筛选变量的方法，同时，介绍一种更加便捷快速完成统计分析的“神器”——风暴统计。

主要内容包括： 一、文献解读 二、利用R语言复现 三、利用在线网站复现   四、小结

一、文献解读

案例文献是沈阳医学院公共卫生学院学者基于SEER数据库的一项回顾性研究，旨在建立一个列线图来预测老年恶性骨肿瘤（MBT）患者的总生存期（OS）。

1. 摘要

背景：恶性骨肿瘤（MBT）是老年患者死亡的原因之一。我们研究的目的是建立一个列线图来预测老年MBT患者的总生存期（OS）。

方法：从SEER数据库下载了2004年至2018年所有老年MBT患者的临床病理数据。他们被随机分配到训练集（70%）和验证集（30%）。采用单因素和多因素Cox回归分析确定老年MBT患者的独立危险因素。基于这些危险因素构建列线图，以预测老年MBT患者的1年，3年和5年OS。然后，利用一致性指数（C指数）、校准曲线和受试者工作曲线下面积（AUC）来评价预测模型的准确性和判别力。决策曲线分析（DCA）用于评估列线图的临床潜在应用价值。根据列线图上的分数，将患者分为高风险组和低风险组。Kaplan-Meier（K-M）曲线用于测试两名患者之间的生存差异。

结果：从SEER数据库下载了2004年至2018年所有老年MBT患者的临床病理数据。他们被随机分配到训练集（70%）和验证集（30%）。采用单因素和多因素Cox回归分析确定老年MBT患者的独立危险因素。基于这些危险因素构建列线图，以预测老年MBT患者的1年，3年和5年OS。然后，利用一致性指数（C指数）、校准曲线和受试者工作曲线下面积（AUC）来评价预测模型的准确性和判别力。决策曲线分析（DCA）用于评估列线图的临床潜在应用价值。根据列线图上的分数，将患者分为高风险组和低风险组。Kaplan-Meier（K-M）曲线用于测试两名患者之间的生存差异。

结论：我们建立了一个新的列线图来预测老年MBT患者的1年，3年，5年的OS。该预测模型可以帮助医生和患者制定治疗计划和后续策略。、

2. 数据介绍

文献共纳入1641名2004-2018诊断的老年MBT患者的临床病理数据。暴露因素包括年龄、种族、性别、诊断年份、组织学类型、分级、分期、原发位置、TNM分期、肿瘤大小、是否化疗、是否放疗及手术方式。

3.研究结果

筛选预测因子只用到训练集数据，通过单因素Cox回归与多因素Cox回归筛选预测变量。这里作者并未使用先单后多或者是逐步回归法筛选预测变量，可能从临床实际考虑的更多。但在实操过程中推荐大家优先考虑逐步回归法进行筛选，这样更有利于构建一个拟合效果好的预测模型。

二、利用R语言复现

下面我们利用从SEERStat中提取的相关数据进行复现。根据文献中的纳入排除标准，提取涉及的相关数据，最终共纳入1，574名患者（原文献1，641）。介于SEERStat数据库会有更新，因此提取的样本量与原文会有所出入，这里请大家多关注统计方法的运用！本次用到的是R版本是4.3.1。

开始回归分析前需要做好前期工作：①导入数据集；②按照7:3的比例将数据集拆分为训练集train和验证集test；③完成均衡性分析。接着再开展今天的回归分析。(详见下方链接）

如何利用R语言拆分数据集并做均衡性检验

1. 安装加载R包并导入数据

这里绘制均衡性表格，主要用到"autoReg"R包，请注意模型变量的筛选是根据训练集进行的，在进行数据拆分后，需要用训练集来进行Cox回归。

install.packages("autoReg")
library("autoReg")

2. Cox回归分析

autoReg包是一款功能强大的R包，可以一步到位实现批量单因素，多因素以及逐步回归法。可以通过调整“threshold= ”来改变限制P值的限制。

#cox回归模型构建
coxmod<-coxph(Surv(train$months,train$status)~age + Race + Sex + Year.of.diagnosis 
              + Histologic.type + Primary.Site + Grade + Stage + t_stage + n_stage
              + m_stage + Surgery + Radiation + Chemotherapy + Tumor_size,
data=train)
summary(coxmod)
#单因素＋P<0.05纳入多因素＋逐步回归后退法
ft3<-autoReg(coxmod,uni=TRUE,threshold=0.05, final= TRUE) 
myft(ft3)

代码解读：autoReg函数可以自动根据设定的阈值控制变量进入多因素回归模型，如果不限制阈值，全部变量进入多因素回归，可将阈值设置为1，如threshold=1；加上“final= TRUE”表示增加逐步向后回归的结果。

R语言结果如下：

三、利用在线网站复现

如果没有代码基础，或者希望通过更便捷的方式完成统计分析，推荐使用这个智能在线统计分析平台——风暴统计。一键成三线表，便捷又快速。

网址：www.medsta.cn（电脑端浏览器打开） 或medsta.cn（medical statistics缩写）

1.进入网站分析模块

电脑端打开风暴统计平台——“风暴智能统计”模块，点击“临床预测模型（最新）”，进入“cox预测模型”页面。

2.完成前期的准备工作

导入数据并根据需要做好整理转换，将数据集拆分为训练集与验证集，做好均衡性分析后，点击进入“Cox构建预测模型”。（详见下方链接）

如何利用R语言拆分数据集并做均衡性检验

3.Cox回归筛选预测变量

点击进入“Cox预测模型”模块，将回归因变量与回归自变量分别选入，此外时间节点的设置与研究数据相对应。

这里要特别注意回归方法的选择，如果仅使用先单后多进行筛选，则逐步回归法选择否，如果选择开展逐步回归法，可以通过多因素回归P值进行阈值的设定，下方的多因素回归结果即逐步回归结果，可以分别查看单因素、多因素（逐步）或先单后多的结果，同样支持导出excel或word三线表结果。

Word版结果如下，为了更加清晰的展示HR值与95%置信区间，下表隐去了β值与SE值，实际上网站给出的统计结果更加全面：

四、小结

通过对比R语言结果与风暴统计结果可以发现，两者的结果完全一致，并且风暴统计给出的统计量更加的全面丰富，且将P值与HR值分为两列，更加美观，如果您在实际中需要用到Cox回归筛选变量，不妨来风暴统计平台试一试！