校正生存曲線在SPSS、Stata及R中的實(shí)現(xiàn)與方法比較*

劉 爽 徐澤鋒 董文革 崔 壯 王俊峰

【提 要】 目的 比較校正生存曲線圖在SPSS、Stata和R軟件中實(shí)現(xiàn)方法的異同，指導(dǎo)臨床醫(yī)師實(shí)現(xiàn)校正生存曲線圖的繪制與解釋。方法 利用歐洲血液和骨髓移植學(xué)會(huì)注冊(cè)數(shù)據(jù)，展示SPSS、Stata和R統(tǒng)計(jì)軟件實(shí)現(xiàn)校正生存曲線的操作過程和結(jié)果，詳細(xì)解析不同統(tǒng)計(jì)軟件采用的校正方法。結(jié)果 SPSS、Stata所提供的校正生存曲線繪制方法默認(rèn)為均值法，R軟件提供的逆概率加權(quán)法對(duì)生存曲線有更好的校正作用。結(jié)論 不同統(tǒng)計(jì)軟件在實(shí)現(xiàn)生存曲線校正時(shí)采用了不同的方法，研究者需要知悉這些方法差異，建議選用R軟件中提供的逆概率加權(quán)法，可以得到更為合理的校正后生存曲線。

繪制生存曲線是生存分析中最基礎(chǔ)的方法，通過將研究對(duì)象分組繪制生存曲線，可以比較生存概率在組間的差異。但是在非隨機(jī)對(duì)照試驗(yàn)設(shè)計(jì)的情況下，生存概率可能受到在各組間分布不均衡的其他協(xié)變量的影響。隨著分子診斷和靶向治療的發(fā)展，以血液腫瘤為典型的惡性腫瘤不同個(gè)體在病理及臨床方面呈高度異質(zhì)性，患者對(duì)治療的反應(yīng)及預(yù)后也相差很大。原始的、未經(jīng)校正的生存曲線圖與Cox比例風(fēng)險(xiǎn)回歸模型的相對(duì)風(fēng)險(xiǎn)結(jié)論有較大差異。近年來，校正后的生存曲線因其可以直觀展示比較結(jié)果和各組生存概率，越來越受到研究者們的重視。

很多學(xué)者提出了校正生存曲線的方法，或?qū)ΜF(xiàn)有的校正生存曲線方法進(jìn)行了歸納和比較[1-4]。如谷鴻秋等介紹了6種校正生存曲線的方法，并比較了這些方法的優(yōu)勢(shì)和局限[4]。然而，大部分臨床研究人員尚未掌握校正生存曲線的繪制，或在實(shí)施統(tǒng)計(jì)方法進(jìn)行數(shù)據(jù)分析時(shí)，僅依賴于統(tǒng)計(jì)軟件提供的功能，在研究中采用哪種方法進(jìn)行生存曲線的校正，很大程度上取決于研究者常用的統(tǒng)計(jì)軟件中提供哪種方法。但是，不同的統(tǒng)計(jì)軟件，在進(jìn)行校正生存曲線的時(shí)候，所用的方法并不相同。本文的目的就是理清常用軟件(SPSS,Stata,R)中提供的生存曲線校正方法，詳細(xì)描述各軟件操作過程，并展示不同調(diào)整方法得到結(jié)果的差異，供臨床研究者掌握校正生存曲線的繪制與解釋。

資料與方法

本研究演示數(shù)據(jù)來自歐洲骨髓和血液移植學(xué)會(huì)(EBMT)開源數(shù)據(jù)，由R軟件中mstate包中加載獲取，數(shù)據(jù)具體信息可見文獻(xiàn)[5-6]，本研究使用的是數(shù)據(jù)名為ebmt2的數(shù)據(jù)。該數(shù)據(jù)中共有8966例接受骨髓移植的患者，時(shí)間起點(diǎn)為患者接受移植的時(shí)間，事件定義為患者死亡，最長(zhǎng)觀察時(shí)間設(shè)置為60個(gè)月。我們選取預(yù)處理是否清除T細(xì)胞(tcd)作為本次研究關(guān)注的分組變量，選取患者年齡、移植年份、供體受體性別匹配作為需要校正的協(xié)變量。在排除分組變量缺失的患者后，共有6110例患者納入最終分析，患者基本信息見表1。該數(shù)據(jù)僅供演示目的，不代表任何真實(shí)情況，也不應(yīng)由此數(shù)據(jù)得出任何臨床結(jié)論。

表1 患者基本信息表

tcd：預(yù)處理是否清除T細(xì)胞

實(shí)例與程序

1.SPSS軟件中提供的方法

SPSS作為簡(jiǎn)單易用的統(tǒng)計(jì)分析軟件，常為臨床研究者的首選。SPSS中提供的校正生存曲線方法是協(xié)變量均值替代法，也可以輸入具體數(shù)值替代均值。這種方法首先建立一個(gè)Cox模型，然后將需要調(diào)整的協(xié)變量的均值帶入Cox模型中，得到不同研究組的生存曲線。根據(jù)在Cox模型中如何處理研究變量(作為協(xié)變量或者作為分層變量)，又可細(xì)分為兩種方法。

(1)研究變量作為協(xié)變量

此方法將研究變量與其他協(xié)變量一起，作為協(xié)變量建立Cox模型。在SPSS(版本25)中的實(shí)現(xiàn)的方法如下：

Analyze→Survival→Cox Regression

Time框：t_60

Status框：s_60，Define Event選項(xiàng)Single value=1，Continue

Covariates框：dissub、match、tcd、year、age

Plots選項(xiàng)：勾選Survival，Separate Lines for框填入tcd，Continue

OK

也可通過下列語句實(shí)現(xiàn)：

COXREG t_60

/STATUS=s_60(1)

/PATTERN BY tcd

/CONTRAST (age)=Indicator

/CONTRAST (year)=Indicator

/CONTRAST (tcd)=Indicator

/CONTRAST (match)=Indicator

/METHOD=ENTER year match age tcd

/PLOT SURVIVAL

/CRITERIA=PIN(.05) POUT(.10) ITERATE(20).

根據(jù)此方法生成的校正后生存曲線見圖1B。

如果不想用均值，而是給定的值，只需要在PATTERN命令中加入variable(value)即可，比如固定年齡組為1(≤20歲)，在/PATTERN與BY tcd之間加入age(1)即可實(shí)現(xiàn)。

(2)研究變量作為分層變量

將研究變量作為協(xié)變量建立Cox模型，需要研究變量也滿足比例風(fēng)險(xiǎn)假設(shè)，所以被認(rèn)為是所有校正方法中最差的方法。可將研究變量作為分層變量建立Cox模型，此方法在SPSS中也很容易實(shí)現(xiàn)，只需要將研究變量填入Cox模型中的Strata項(xiàng)即可。或運(yùn)行下列語句：

COXREG t_60

/STATUS=s_60(1)

/STRATA=tcd

/CONTRAST (age)=Indicator

/CONTRAST (year)=Indicator

/CONTRAST (match)=Indicator

/METHOD=ENTER year match age

/PLOT SURVIVAL

/CRITERIA=PIN(.05) POUT(.10) ITERATE(20).

根據(jù)此方法生成的校正后生存曲線見圖1C。

A原始生存曲線；B研究變量作為協(xié)變量；C研究變量作為分層變量

2.Stata中提供的方法

與SPSS類似，Stata(版本15)中校正生存曲線方法也是協(xié)變量均值替代法，但是在如何處理研究變量上略有不同，也可細(xì)分為兩種方法。

(1)根據(jù)研究變量分亞組

Stata中校正生存曲線是通過sts graph命令中的adjustfor選項(xiàng)來實(shí)現(xiàn)的[7]，具體方法如下：

stset t_60,failure(s_60==1)

sts graph,by(tcd) adjustfor(year match age)

這種方法類似于SPSS中將研究變量作為協(xié)變量的方法，但是不是擬合一個(gè)Cox模型，而是根據(jù)研究變量的分類，將數(shù)據(jù)分為N個(gè)亞組(N=研究變量分類數(shù))，在每個(gè)亞組中分別擬合一個(gè)Cox模型[7]。在每一個(gè)亞組模型中，協(xié)變量的系數(shù)可以不同。根據(jù)此方法生成的校正后生存曲線見圖2A。

需要注意的是，Stata中協(xié)變量并不是取均值，而是取值為0。對(duì)于連續(xù)型協(xié)變量，如果想要根據(jù)具體參考值進(jìn)行校正，需要提前生成一個(gè)新的協(xié)變量generate var_new=var_old - reference_value，這里reference_value為參考值。

(2)研究變量作為分層變量

將Stata中根據(jù)研究變量分五組進(jìn)行校正生存曲線的程序中的by()命令改為strata()命令，就可以將研究變量作為分層變量，擬合1個(gè)Cox模型，協(xié)變量在每個(gè)分層中的系數(shù)相同，但是每個(gè)分層有不同的基礎(chǔ)風(fēng)險(xiǎn)函數(shù)。根據(jù)此方法生成的校正后生存曲線見圖2B。

stset t_60,failure(s_60==1)

sts graph,strata (tcd) adjustfor(year match age)

A根據(jù)研究變量分亞組；B研究變量作為分層變量

3.R軟件中提供的方法

R軟件由于其強(qiáng)大的功能和擴(kuò)展能力，備受統(tǒng)計(jì)研究者的推崇。R軟件的統(tǒng)計(jì)分析主要是基于功能包(package)來實(shí)現(xiàn)的，所以R軟件進(jìn)行生存曲線校正的方法，將按照不同的功能包來介紹。

(1)survival包中survfit命令

與前面所介紹的其他軟件相同，在R軟件最常用的survival包的survfit命令，也提供了協(xié)變量均值替代法進(jìn)行生存曲線校正。在利用coxph命令建立Cox模型之后，只需在survfit命令中newdata選項(xiàng)給出新的輸入數(shù)據(jù)，就可以得到協(xié)變量均值替代法校正后的生存曲線。

(2)survminer包中g(shù)gadjustedcurves命令

在R軟件的survminer包中，提供了更多對(duì)生存曲線進(jìn)行校正的方法[8]。在早期的survminer包中，這個(gè)命令叫做ggcoxadjustedcurves，從0.4.1版本以后，舊的命令被ggadjustedcurves命令取代，并且加入了更多新的功能選項(xiàng)。survminer包中g(shù)gadjustedcurves命令主要是根據(jù)Therneau,Crowson和 Atkinson文章中討論的生存曲線校正方法[9]，在R軟件中進(jìn)行實(shí)施。在當(dāng)前版本中(0.4.3)，共提供了四種校正方法：?jiǎn)我磺€方法(single)，平均方法(average)，條件方法(conditional)和邊際方法(marginal)。需要特別注意的是，由于survminer包的開發(fā)者的失誤，在ggadjustedcurves命令中條件方法和邊際方法的調(diào)用被互換了，與Therneau原文中的定義剛好相反。本文作者已經(jīng)與survminer開發(fā)者取得聯(lián)系，并得到確認(rèn)，這一現(xiàn)象確實(shí)是開發(fā)過程中發(fā)生的錯(cuò)誤，將在之后的版本中進(jìn)行更正。所以當(dāng)使用ggadjustedcurves命令時(shí)，請(qǐng)注意版本號(hào)和解釋文檔，以免誤用。本文為保持與Therneau文章的一致性，將以該文中的方法名稱為準(zhǔn)。

單一曲線方法僅繪制一條生存曲線，無法比較研究變量的各個(gè)分組，平均方法與前文中的均值替代法相同，所以這里只著重介紹條件方法和邊際方法。

①條件方法

ggadjustedcurves中所謂的條件方法，來源于Therneau的文章，代表先計(jì)算每個(gè)研究對(duì)象的生存曲線，然后計(jì)算每組平均生存概率，并不是谷鴻秋文章提到的條件概率校正法。

mod1<-coxph(Surv(t_60,s_60>0) ～tcd+year+match+age,data=data)

ggadjustedcurves(mod1,data=data,method=″conditional″,variable=″tcd″,ylab=′Survival probability′,xlab=′Months′,size=1.5)

由于前文中提到的功能包的開發(fā)者的錯(cuò)誤，這里如果想要調(diào)用條件方法，在當(dāng)前版本中，應(yīng)該選擇method=”marginal”。用條件方法進(jìn)行生存曲線校正的結(jié)果見圖3A。

②邊際方法

ggadjustedcurves中的邊際方法，應(yīng)用逆概率加權(quán)方法平衡研究變量各組研究對(duì)象的分布。

mod1<-coxph(Surv(t_60,s_60>0)～tcd+year +match+age ,data=data)

ggadjustedcurves(mod1,data=data,method=″marginal″,variable=″tcd″,ylab=′Survival probability′,xlab=′Months′,size=1.5)

同樣，由于前文中提到的功能包的開發(fā)者的錯(cuò)誤，這里如果想要調(diào)用邊際方法，在當(dāng)前版本中，應(yīng)該選擇method=″conditional″。用邊際方法進(jìn)行生存曲線校正的結(jié)果見圖3B。

(3)RISCA包中ipw.survival命令

R軟件中的RISCA包(0.8.2)也提供了逆概率加權(quán)法進(jìn)行生存曲線校正的功能，可以用ipw.survival命令繪制生存曲線[10]，由于與4.2.2中的邊際方法相同，這里就不做演示了。值得一提的是，在RISCA包中，除了校正生存曲線，開發(fā)者還提供了相應(yīng)的校正log rank檢驗(yàn)的命令ipw.log.rank，可以根據(jù)校正后的生存曲線用log rank檢驗(yàn)比較生存曲線的差異。

A 條件方法；B邊際方法

結(jié) 論

從以上的展示和比較中，我們可以看到，不同的統(tǒng)計(jì)軟件，在進(jìn)行生存曲線校正時(shí)使用的方法并不相同。如果不加說明地使用校正生存曲線匯報(bào)研究結(jié)果，可能會(huì)產(chǎn)生誤解。研究者在論文中使用校正生存曲線時(shí)，應(yīng)報(bào)告使用的軟件，以及具體的方法，以便于審稿人、編輯和讀者更好地理解。如使用SPSS或Stata進(jìn)行均值法校正生存曲線，應(yīng)盡量避免使用均值，而是使用有實(shí)際意義的特定值。然而，均值法校正生存曲線由于方法本身的缺陷，并沒有起到很好的校正作用[4]，得到的校正后曲線的意義有限。所以在有條件的情況下，筆者建議選用 R軟件中提供的逆概率加權(quán)法，可以得到更為合理的校正后生存曲線。本研究為繪制校正的生存曲線圖提供了很好參考，可以解決我國(guó)相當(dāng)一部分科研人員的實(shí)際問題。