如何使用R来做Cox回归?来看详细教程!
回归家族三兄弟
Cox回归是个有魔性的统计学方法。介绍Cox回归分析之前我们先介绍一下回归家族的另外两个成员:大哥——线性回归,二哥——逻辑回归。
大哥线性回归是个意气汉子,办事情从来直来直去。二哥逻辑回归最喜欢干搬箱子的活, 把不同种类的箱子搬来搬去搬来搬去。三弟Cox回归最魔性,他可以看到一段时间以后将会发生些什么。
三兄弟着眼的都是“关系”,而他们擅长的事务大有不同:大哥最擅长解决“连续”的问题;二哥最擅长解决“分类”的问题;三弟最擅长解决“和时间有关系”的问题。
举个栗子。
隔壁张大爷有个苹果园,三个瓜娃子从小到大没少偷苹果吃,张大爷都睁一只眼闭一只眼。年复一年的,张大爷渐渐老了,打理一个大果园有些力不从心,果园的收成慢慢地不如从前,大爷无计可施,只好在一天找来了三兄弟。
张大爷慢慢在果园里踱来踱去,说“三个娃啊,你们都长大了,有些事我弄不明白只好来问你们了。我自己感觉啊,种出来的苹果每年口感都不太一样,有的年份更甜一些,有的年份比较一般。我心里大概有个谱,一年里晴的时候多吧,春天冷一些吧,苹果好像就甜一些。不知道你们能不能把这事给我讲明白?”
大哥马上跳了出来:“大爷,这个问题让我来!甜度与阳光和温差有直接关系,我们可以用线性分析来做。每天的日照时间长度和早晚温差都不同,只要有足够的以前年份的数据就可以找出日照和温差的影响,然后我们就可以根据今年的阳光和温差情况预测今年的苹果甜度了。”
回归家族大娃,手里高高举着牌子“线性回归:专门解决因变量是“连续”的问题。”
张大爷点头:“大娃你说的好。那知道了甜度之后,我还有个问题。有的年份苹果卖的价钱好,有的年份卖的价钱就一般。 你也知道,有的年份苹果甜,有的不甜;有的年份来收的人多,有时候来的人少;有的年份能卖到出口经销商那去,有时候不能。你们说,这些情况和苹果的价钱到底有没有联系呢?”
话音未落,二哥就跳来出来:“这个我懂。大爷其实您每年都对苹果能卖多少钱心里有个打算吧,可能有的年份能比您心理价位卖的还高一些,有的年份比您心理价位低一些。咱不看每年具体苹果卖多少钱,年景不同没法比,就看这个价钱是不是比您的心理价位低,好吧。如果实际卖价比您心理价位低了,我们就说卖的不好,如果实际卖价比您心理价位高,我们就说卖的好。
那大爷您这个问题涉及到的就是各种情况,比如您提到的甜不甜啊、有没有大经销商来采购啊,和您的苹果卖的好不好的关系。大哥的方法太直,不适合这个问题。咱们要把方法转化一下,比如用一个指数函数…” 二哥拉住大爷的袖子“哎大爷别走啊,我还要给你解释一下指数函数是啥呢…”
好容易解释清楚了指数函数是什么,二哥默默举起了牌子“逻辑回归:专门解决因变量是“分类”的问题。”
大爷袖着手走来走去,“娃们啊,我还有个问题。咱们这个苹果熟透了,就从树上落下来,可是落下来的很多都磕坏了。可苹果没熟透的时候如果摘下来,味道又不够。每年果园里摘苹果都要讲究,不能太早也不能太晚。咱们怎么能知道什么时候苹果都会掉下来,什么时候还没熟呢?”
三娃推开了两个哥哥,“大爷,这个问题我懂。时间越长,苹果落的越多。” 大爷说,“也不能这么说,风速,日晒,化肥都会影响苹果的成熟。” 三娃继续说:“风速,日晒,化肥都对苹果成熟有影响。当苹果完全成熟这一事件发生时,它从树上落下,那我们就知道,它从结果子到落下一共用了多久。风速、日晒和化肥对苹果落下这件事的影响就体现在这个用了多久落下来的时间上。”大爷点点头,“说的对啊,三娃。”
三娃举起牌子“Cox回归:专门解决各因素影响下,某事件在一段时间内发生概率的问题。”
三娃的本领
三娃的本领:不同因素,对某一事件在时间t上发生的概率的影响作用 。
在生存问题中,我们关注的是两件事:事件到底发生了没有(事件发生的概率),事件用多长时间发生的(生存时间)。这两件事是一起存在的,只有当事情发生了,才有生存时间可谈。
回到果园的问题上,从苹果树结果子开始,我们就天天盯着树上的苹果,一天两天,一个月两个月(事件一直没发生)。终于有一天,苹果从树上“哐”掉了下来(事件终于发生了),我们乐疯了好吗,终于可以记下从结果子到现在的时间了(生存时间)。这个生存时间和苹果落地事件是联系在一起的。
在整个过程中,风力、日晒和肥料可能会对苹果落地事件有影响作用,即可能对(苹果在果树上的)生存时间有影响作用。这就是Cox回归研究的内容。
总结一下,在这个果园模型中,苹果落地是事件,从结果到苹果落地的时间为生存时间,这是两个因变量;风力、日晒和肥料是可能的影响因素,是三个自变量。
那么我们可以假设,时间的改变为(∆t), 生存时间为T, 那么在t时刻单位时间内事件发生概率为:
同时我们知道,苹果落地这一事件(Y)和时间(t)以及三个因素日照(x1)、风力(x1)和肥料(x1)有关。假设存在一个固定的初始概率h0(t),代表基本状态下苹果落地的概率, 某时间下苹果落地这一事件发生的概率与各个因素之间的关系为:
对于一个固定的问题(同一个果园,同一批果树),如果我们拿出两个观察到的样本i和i´,它们各自的自变量值为{Xi1,Xi2…Xik}和{Xi´1,Xi´2…Xi´n},那么在时间t时,在这两个样本上发生该事件(苹果坠落)的概率是等比的,写为:
其中,
这一“等比性”是Cox回归的一大特点。
在R中进行Cox回归
今天我们来讲一讲如何在R中进行Cox分析。
为什么要用R呢?
因为R中Cox函数的表示方法都与Cox分析的数学表达式很类似,能帮助我们直观的看到Cox回归中的各项参数。另外,R中可以根据需求更改Cox回归中的自变量(即x)。(当然我一定不会说最重要的是R里久经考验的安装包==)
好啦,我们来看看R中要怎么做Cox回归分析。当然,现在我们要进入到正经的工作状态,不用苹果举例子,回到我们的临床数据上。
首先,我们要下载安装并在系统中载入 “survival”和“survminer”包。
survival包是R中经典的分析包之一,主攻就是生存分析,里面有很多与Cox回归相关的命令。“survminer”是与它配合使用的分析包, 如果小伙伴们想要根据生存数据作图的话,survminer是个好选择。
图1. 下载安装及载入“survival”包
在“survival”包中有很多数据,我们可以调出数据列表:
图2. 代码及数据列表结果
我们选用肺癌病人数据,并查看数据前几行。查看前几行是我们接触新数据到第一步,它可以帮助我们直观的了解数据结构。这一癌症数据是survival包中自带的默认数据。我们可以看到,数据格式是每一个病人数据为一个横行,采集的病人数据有inst、time等,在数据采集的过程中有空白项(NA)的出现。
图3. 代码及癌症病人数据前几行直观图
在R中的数据各项指标都已经数字化了,age年龄下的数字好理解,但是inst,status中的1,2都分别代表了什么呢?我们可以查看R中的帮助文件来了解数据详情:
图4. 调取数据帮助文件代码及数据解释
对于 “survival”包中的癌症患者数据,我们想要研究的是,随着时间的改变,性别(sex)与年龄(age)对患者死亡情况的影响。 在这个临床问题中,性别与年龄对应了果园问题中的日照、风力和肥料,它们都是影响因素,是自变量;患者的治愈情况与苹果坠落情况对应,都是我们要研究的具体事件,是因变量。
在R中构建的函数格式与数学公式极为相似。我们首先把问题简化为病人死亡事件和年龄(age)、性别(sex)之间的关系。简化模型是为了方便介绍代码的构成以及如何阅读回归结果,小伙伴们还是要在回归模型中带入全部的影响因子啊。回归分析的结果我们着重要看什么呢?要看各影响因素的影响是不是显著的(P值是不是小于0.05或0.1),以及各因素的影响因子值(coef)是多少。
如果在Cox模型中自变量为age,sex,因变量为status(status=2代表病人死亡),R中代码和回归结果如下:
图5. 简化版Cox回归经典代码与回归结果
根据这份肺癌病人数据,我们可以做更全面的分析,理论结构和简化版完全一致,只是多出了两个自变量: ph.karno和wt.loss。我们在R中的代码和回归结果如下:
图6. 完整版Cox回归代码与回归结果
这样Cox回归分析就做完了吗?当然不是啦,小伙伴们还记得上学时候老师经常说考试做完卷子要干嘛?回头检查啊!不要扔下笔就跑了啊。
做Cox回归分析很重要的的一步就是检验模型。之前我们说过,Cox回归一个很重要的特点是这个回归模型具有等比性。因此,我们得到一个Cox回归模型后,首先要检验的就是它的“等比性”。
图7. Cox回归模型检验及结果
我们发现,检验结果中ph.karno具有显著性(p<0.05),说明ph.karno这一自变量不符合“等比性”要求,这一回归模型不准确。 当某一因素不符合要求时,我们的处理方法是:将此因素分层分析。
我们采用这种处理方法的原因是,当这一因素不满足“可比性”时,它在数据中不均一。不均一是什么意思呢?通俗的讲就是不均匀,分层了,比如水就是均一溶剂,水油混合物(想想家里梳妆台上的水油卸妆液)是不均一的分层溶剂,在统计学上是一个道理。而且从分析结果上看,它在数据中的差别具有统计学显著性。在这种情况下,已经不能把此因素作为整体分析(水和油不能混在一起,要分开),必须进行分层。我们对ph.karno分层分析如下:
图8. 分层Cox回归模型代码与结果
我们也要检验分层后的Cox模型中其他自变量是否已经符合“可比性”。
图9. 分层Cox回归模型检验代码及结果
可以看出,分层后的Cox回归模型各影响因素都通过了可比性检验,这一分层方法是正确的。那么分层后的Cox回归模型该如何解读呢?
图10. 分层Cox回归模型综合结果
首先,分层之后因为数据缺失损失了14个数据,样本总量为214,病人死亡事件152起。其次,在分层后Cox回归分析采纳的三个自变量age,sex,wt.loss中我们发现性别(sex)与病人死亡事件有非常显著(P<0.001)的负面关系(性别因素的影响因素为负)。sex为定性因素,根据数据索引,sex=1为男性,sex=2为女性,同时 status=1时病人生存, status=2时病人死亡。
根据前文对Cox回归的描述,我们研究的是多因素影响下,在一段时间内某事件发生概率的问题,在本组数据中是肺癌病人的死亡率问题,我们得出的结论是,男性病人死亡率高于女性病人。
以上,就是Cox回归分析对生存数据的分析过程,由问题起源到数学原理再到编程应用,小伙伴们学会了嘛?有任何意见和建议,欢迎给我们留言哦。
相关阅读
关注医咖会,学习统计学~
有临床研究设计或统计学方面的难题?快加小咖个人微信(xys2016ykf),加入医咖会统计讨论群,和小伙伴们一起交流学习吧。我们诚邀各位小伙伴加入我们,一起创作有价值的内容,将知识共享给更多人!
点击左下角“阅读原文”,看看医咖会既往推送了哪些研究设计或统计学文章。