校正生存曲線在SPSS、Stata及R中的實現與方法比較*

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

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

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

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

資料與方法

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

表1 患者基本信息表

tcd：預處理是否清除T細胞

實例與程序

1.SPSS軟件中提供的方法

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

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

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

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

也可通過下列語句實現：

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).

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

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

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

將研究變量作為協變量建立Cox模型，需要研究變量也滿足比例風險假設，所以被認為是所有校正方法中最差的方法。可將研究變量作為分層變量建立Cox模型，此方法在SPSS中也很容易實現，只需要將研究變量填入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).

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

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

2.Stata中提供的方法

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

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

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

stset t_60,failure(s_60==1)

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

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

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

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

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

stset t_60,failure(s_60==1)

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

A根據研究變量分亞組；B研究變量作為分層變量

3.R軟件中提供的方法

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

(1)survival包中survfit命令

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

(2)survminer包中ggadjustedcurves命令

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

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

①條件方法

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)

由于前文中提到的功能包的開發者的錯誤，這里如果想要調用條件方法，在當前版本中，應該選擇method=”marginal”。用條件方法進行生存曲線校正的結果見圖3A。

②邊際方法

ggadjustedcurves中的邊際方法，應用逆概率加權方法平衡研究變量各組研究對象的分布。

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)

同樣，由于前文中提到的功能包的開發者的錯誤，這里如果想要調用邊際方法，在當前版本中，應該選擇method=″conditional″。用邊際方法進行生存曲線校正的結果見圖3B。

(3)RISCA包中ipw.survival命令

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

A 條件方法；B邊際方法

結 論

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