手把手教Stata做生存分析:K-M曲线绘制和Logrank检验

专题合集更多教程

对于生存数据,1958年,E. L. Kaplan 和 Paul Meier 两位教授介绍了一种全新的、解决随访期间右删失 (right censoring) 问题的生存分析方法,被称作Kaplan-Meier方法。这种方法精确地记录并利用每个个体发生终点事件的具体时间,在任何一个终点事件发生的时间点计算出一个新的、基于之前所有信息的总生存率 (Cumulative survival) 。

 

相比于之前使用的寿命表法(Life-table method),这种方法更加充分地利用了信息,给出更准确的统计量。同时,作为一种非参数估计方法,不要求总体的分布形式,因此非常适合生存分析时使用。Kaplan-Meier 曲线(简称K-M曲线) 还可以很直观地表现出两组或多组的生存率或死亡率,非常适合在文章中进行展示。因此,K-M曲线也成为了临床研究中最常用的方法之一。

 

Why Stata?

 

Stata软件在进行生存分析的过程中具有很强大的功能。无论是 K-M曲线,还是Cox回归分析,甚至是一些更加复杂的参数分析,Stata都可以轻松完成。

 

相比于SPSS,Stata的可重复性更强,图像更加美观;而相比于R语言,Stata的代码又更加简便、易懂、上手快,甚至可以完全使用窗口菜单完成,非常适合有科研追求的医生入门。

 

今天,我们就一起来学习一下生存分析中的第一步、也是最重要的步骤之一:K-M曲线的绘制和Logrank检验。

 

我们将使用Stata自带的一个模拟的药物临床试验的数据集进行所有的演示,请大家在Command对话框中输入webuse drugtr以调入这个数据集。

 

屏幕显示:

 

请注意:Stata已经将这个数据集设置成了生存数据的格式,导入数据集后,请大家首先输入stset, clear命令恢复成普通的数据格式。这样才是我们在临床研究中见到的数据结构,我们将在Step 2中学习如何将其转换成生存数据格式。

 

Step 1

数据集的初步观测

 

首先,我们来观察一下数据集。大家可以在command line中输入list 查看所有数据,也可以输入list in 5/10查看第5到第10个观测值,屏幕显示:

 

我们可以看到,在这个数据集中,共有4个变量:studytime, died, drug, age。大家可以在command line中输入“codebook+空格+变量名字”,查看变量的标签、范围、是否有缺失值等基本特征。例如我们输入codebook drug, 屏幕显示:

 

我们可以看到共有48个数据,不重复的值有两个,分别是0 (安慰剂,从第一行的标签"Drug type (0=placebo)"看出来) 或 1 (试验药),没有缺失值。

 

同理,我们可以查看其它4个变量: studytime (标签:Months to death or end of exp. ,代表了这位患者的随访时间),died (标签:1 if patient died), age (标签: Patient's age at start of exp.)。

 

在基本了解数据后,我们可以正式开始数据分析了。

 

Step 2

声明数据为生存分析数据 & 数据再观测

 

Stata要求我们他这是一个生存分析的数据集,因此,我们需要至少告诉Stata如何判断终点事件(Failure variable)、如何判断随访时间(Time variable)。

 

在窗口菜单下,我们可以通过Statistics > Survival analysis > Setup and utilities > Declare data to be survival-time data找到如下的对话框。

 

在Time variable一栏选择数据集中的 "studytime"变量,在failure event一栏选择数据集中的"died"变量,并在Failure values的框中输入1。若多个值代表终点事件的发生,我们也可以在Failure values一栏输入多个值,每个值之间用空格隔开。

 

你可能会注意到,在Main菜单界面上还有一栏"Multiple-record ID variable"。Stata默认数据集中的一行数据代表一位患者。而在某些研究中,一位患者可能会有多行数据, 这一栏便是告诉Stata用哪个变量区分不同患者,通常选择数据集中代表患者ID的变量。然而,由于这种情况在实际的临床分析中并不常见,我们不详细阐述。在我们这个例子中,因为每个患者只有一行数据,我们可以不指定ID 。 

 

点击OK后,屏幕上Results窗口内展示出如下画面。

 

如果你愿意使用代码,或者希望你的合作者能够重复你的结果,请使用上图第一行的代码:stset studytime, failure(died==1)。这是你在窗口菜单操作时Stata在背后执行的命令,我们也可以直接输入这行命令得到相同结果。

 

我们可以看到一共有48个患者,31人发生终点事件,共有744段person-time,随访时间最长的人是t=39。

 

使用stsum以及stdescribe命令,我们可以得到更多的关于这个数据集的信息。请注意:必须要在指定数据集为生存分析数据集之后 (换句话说,stset之后才能使用任何其他的st开始的命令)。

 

屏幕显示:

 

 

通过sts list命令可以展示Stata绘制K-M曲线所用的表格,我们可以看到每一个时间节点上的at-risk的人数、发生终点事件的人数、失访的人数,以及总的生存率。我们只展示前5行。

 

 

Step 3

K-M曲线的绘制

 

接下来,我们可以画K-M曲线。

 

在菜单栏中选择Statistics > Survival analysis > Graphs > Survivor and cumulative hazard functions,我们可以选择生存曲线(Function选项中Graph Kaplan-Meier survivor function) 或死亡曲线(Function选项中Graph Kaplan-Meier failure function),并且可以在下方Grouping variables处选择按照drug的数值分组绘图。

 

点击OK后,屏幕上会出现代码sts graph, by(drug)以及图像(下图)。

 

 

这是最基本、最简单的K-M曲线形式,在Stata中,我们还可以做很多改变,让图像变得更美观,表达更多的信息。

 

例如,在刚刚 Graph the survivor and cumulative hazard functions的窗口菜单中,利用第二个标签页(if/in)可以对数据库进行符合某种条件的筛选。例如我们只关注age<50岁的人。

 

点击OK后,屏幕上会出现代码sts graph if age<50, by(drug) 以及图像(下图)。

 

利用第三个标签页(At-risk table)勾选Show at-risk table beneath graph (截图省略),可以在x轴下方显示每个时间节点上at-risk的人数,便于读者更好了解K-M曲线在某一时间节点下降的原因。点击OK后,屏幕上会出现代码sts graph, by(drug) risktable以及图像(下图)。

 

如果把代码中的risktable改成atrisk,我们可以将每个时间点at risk的人数直接标记在曲线上。下图代码: sts graph, by(drug) atrisk

 

利用剩下的几个标签页,我们还可以增加置信区间(CI)、更改坐标轴标签(Y axis & X axis)、标签(Legend)、标题(Titles), 图像整体风格(Overall)等等,大家可以自己点击操作进行尝试。我们在此通过代码做一个示范,大家可以根据需要更改代码的内容,适应自己的研究需要 (下列图像选项都可以不用代码,而通过窗口菜单完成)。

 

这张图的代码是:

 

label define drug_label 0 "安慰剂" 1 "试验药"

label values drug drug_label

sts graph, by(drug) ci atrisk xlabel(0(5)40) ylabel(0(0.2)1) legend(cols(2)) xtitle("x轴的标签") ytitle("y轴的标签") title("这里是标题") subtitle("这里是副标题") caption("这里是注释") note("这里可以写数据来源")

 

解释:

 

代码第1~2 行:定义一个名为drug_label的标签(你可以随意取名字),让数值为0标记为“安慰剂”,数值1标记为“试验药”。用这个标签去把drug这个变量的数值标记成“安慰剂”或者“试验药”。

 

请注意:代码第1、2行也可以通过Data > Variables Manager来通过窗口菜单操作: 

 

代码第3行:ci代表置信区间;如果大家对于xlabel, legend, title等命令有疑问,可以参考我们之前推送过的文章:如何用Stata作漂亮的图?来看超详细教程!

 

请注意:第3行的代码需要写在同一行,不要手动换行!如需换行,请在需要换行的地方加上“///”,如下图:

 

 

Step 4

检验组间差别

 

K-M曲线完成后,我们可能用肉眼看出了两组有显著差别。但是,我们还需要用统计学方法检测一下(很多时候,当效应值很小时,即使肉眼看不出来,也可能会有显著差别,反之亦然)。我们可以通过输入sts test加上分组变量的名字(比如这个例子中就是sts test drug)来进行检测。

 

Stata默认使用Log-rank test进行检测,结果如下:

 

在这个检验中,P<0.0001,因此我们拒绝零假设(两个药物组没有显著性差别)。因此得出结论,试验药能够显著提高生存率。

 

在下一期中,我们将继续学习如何利用Stata进行Cox回归的分析。敬请期待!

 

附录:do file

 

最后附上本次的do file,请大家把代码拷贝至Do-file Editor中(下图左4图标"Do-file Editor") 

 

并用光标选中希望运行的代码(下图我选中了前3行代码)

 

点击Do(上图右上角的图标"Do")

 

############Do-file开始############

clear all //清除系统现存的数据

set more off //让Stata一口气输出一个命令的所有结果

 

webuse drugtr //调用数据, 请注意需要联网

stset, clear //取消这个数据集默认的生存数据格式,便于我们后面演示

list in 1/10 //查看第1-10个观测值

codebook studytime died drug age //查看这4个变量的基本信息

stset studytime, failure(died==1) time0(t0) scale(1) //告诉Stata这是生存数据

stsum 

stdescribe

sts list

 

sts graph, by(drug) //最基本的样式

sts graph if age<50, by(drug)  //限制age<50

sts graph, by(drug) risktable //加上risk table

sts graph, by(drug) atrisk //加上at risk的人数

 

label define drug_label 0 "安慰剂" 1 "试验药" //创建标签

label values drug drug_label //给变量的数值加上标签

sts graph, by(drug) ci atrisk xlabel(0(5)40) ylabel(0(0.2)1) legend(cols(2)) xtitle("x轴的标签") ytitle("y轴的标签") title("这里是标题") subtitle("这里是副标题") caption("这里是注释") note("这里可以写数据来源") //在图像上加上一些补充信息

 

sts test drug //默认使用Log-rank test检验两组是否有区别

############Do-file结束############

描述问题
选择一个标签 (请选择一个与您问题最相符的标签)
提交问题
我要提问
描述问题
选择一个标签 (请选择一个与您问题最相符的标签)
提交问题
描述问题
选择一个标签 (请选择一个与您问题最相符的标签)
    提交问题