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

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

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

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

資料與方法

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

表1 患者基本信息表

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

實例與程序

1.SPSS軟件中提供的方法

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

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

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

Analyze→Survival→Cox Regression

Time框：t_60

Status框：s_60，Define Event選項Single value=1，Continue

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

Plots選項：勾選Survival，Separate Lines for框填入tcd，Continue

OK

也可通過下列語句實現(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)即可實現(xiàn)。

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

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

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àn)的[7]，具體方法如下：

stset t_60,failure(s_60==1)

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

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

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

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

將Stata中根據(jù)研究變量分五組進(jìn)行校正生存曲線的程序中的by()命令改為strata()命令，就可以將研究變量作為分層變量，擬合1個Cox模型，協(xié)變量在每個分層中的系數(shù)相同，但是每個分層有不同的基礎(chǔ)風(fēng)險函數(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)計研究者的推崇。R軟件的統(tǒng)計分析主要是基于功能包(package)來實現(xiàn)的，所以R軟件進(jìn)行生存曲線校正的方法，將按照不同的功能包來介紹。

(1)survival包中survfit命令

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

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

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

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

①條件方法

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

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ā)者的錯誤，這里如果想要調(diào)用條件方法，在當(dāng)前版本中，應(yīng)該選擇method=”marginal”。用條件方法進(jìn)行生存曲線校正的結(jié)果見圖3A。

②邊際方法

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

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ā)者的錯誤，這里如果想要調(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檢驗的命令ipw.log.rank，可以根據(jù)校正后的生存曲線用log rank檢驗比較生存曲線的差異。

A 條件方法；B邊際方法

結(jié) 論

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