生存分析（survival analysis）即是將終點事件的出現與否和出現終點事件所經歷的時間結合起來的一種統計分析方法，生存分析通常研究的終點事件是死亡，生存分析由此得名。但生存分析可更廣泛地運用于惡性腫瘤、慢性疾病或其他情況的隨訪研究中的事件分析，比如疾病的發生、復發、轉移、傷口的愈合、某種癥狀的消失等。生存資料的分析主要特點就是考慮每個研究對象出現某一結局所經歷的時間。生存曲線即是以生存時間為橫軸，生存率為縱軸，將各個時間點對應的生存率連接在一起的曲線圖^[1-2]。

1   生存分析中幾個重要的基本概念

1.1   生存時間

生存時間（survival time）也是一個廣義概念，泛指所關心的某現象的持續時間，即隨訪觀察持續的時間，常用符號t表示。生存時間分為兩種類型：（1）完全數據（complete data）：指從觀察起點到發生“死亡”事件所經歷的時間。提供了觀察對象確切的生存時間。（2）截尾數據（censored data）：亦稱截尾值（censored value）或終檢值。指從觀察起點到發生非“死亡”事件所經歷的時間。

1.2   截尾

生存結局分為“死亡”與“截尾”兩類，“死亡”是感興趣的終點時間，其他終點事件或結局都歸為截尾。

1.3   死亡概率

死亡概率（probability of death）表示單位時間段開始存活的個體，在該段時間內死亡的可能性。符號q表示。“q=某年內死亡人數÷某年年初人口數”。

1.4   生存概率

生存概率（probability of survival）表示單位時間段開始存活的個體，到該段時間結束時仍存活的可能性。符號p表示。p=某年活滿一年人口數÷某年年初人口數。P=1-q。

1.5   生存率

生存率（survival rate,survival function）表示觀察對象經歷t_k個單位時間段后仍存活的可能性。若無截尾數據，則S（t_k）=P（T>t_k）=t_k。其中0<S（t）≤1。若有截尾數據，須分時段計算生存概率。假定觀察對象在各個時段的生存事件獨立，應用概率乘法定理：S（t_k）=P（T>t_k）=P₁·P₂…P_k^！，P_i為某時段的生存概率，故生存率又稱累積生存概率（cumulative probability of survival）。

1.6   生存曲線

生存曲線（survival curve）：生存時間為橫軸，將各時點所對應的生存率連接在一起的曲線圖，樣本量小時生存曲線呈階梯形，樣本量足夠大時，形成光滑的曲線。

1.7   中位生存時間

中位生存時間是指50%觀察對象能存活的時間。

2   生存分析的統計學方法

由于生存時間一般不呈正態分布，而且需要考慮截尾數據，生存分析有其獨特的統計學方法。常用的統計學方法有以下幾種。

2.1   描述性分析

根據樣本生存資料估計總體生存率及其他有關指標（如中位生存時間等）。常采用Kaplan-Meier法（乘積極限法）進行分析。對于頻數表資料則采用壽命表法進行分析。計算生存率需要考慮時間順序。

2.2   比較分析的方法

對不同組生存率進行比較分析，常采用非參數的log-rank檢驗，檢驗無效假設使兩組或多組總體生存時間分布相同。

2.3   影響因素分析

通過生存分析模型來探討影響生存時間的因素，常用的方法為COX比例風險模型。

3   基于Stata軟件的統計學實現生存分析

（筆者注：以下所舉實例數據全部來自于陳峰教授主編《現代醫學統計方法與Stata應用（第2版）》，相關Stata命令及結果解釋大部分來自于這本書，其中部分命令有少許改動。陳鋒教授主編的這本書通俗易懂，感興趣的讀者可以找來一讀。）

3.1   生存資料的定義

在對隨防資料進行生存分析之前，需先將該數據庫定義為生存資料數據庫，其命令是：

stset時間變量[,failure (截尾變量==#)]

其中，選擇項failure（截尾變量==#）規定截尾變量取值為“#”時研究對象出現預期結果，沒有該選擇項時，Stata以所有不等于0 的非缺失值為出現預期結果。對數據庫進行定義時必須注意變量順序，命令stset后的變量順序依次為時間變量、截尾變量。定義數據庫后，系統自動產生四個變量：

_st /* 數據庫中該條記錄是否被定義為生存資料

_d /* 該條記錄是否出現預期結果

_t /* 觀察對象被隨訪的時間

_t0 /* 觀察對象第一次被觀察到的時間（開始過程的時間為0）

例1 某醫院泌尿外科于1979～1982年間施行19例腎移植手術，擬了解腎移植后患者的生存時間（d）。規定隨訪開始時間為患者術后1 d，預期結果為該患者因與腎移植有關的各種原因的死亡。后改進手術方式，于1983～1986 年又施行14例，資料如下（有+的數據表示該患者截尾）。計算各組的生存率及可信區間（資料已存入文件例1.dta）。

一般手術: 3 9 15 20 20 26 30 41 46 64+64 135 223 365 450 596+680+900+900+

改進手術: 10 390+518+70+70+120 475+225 801+366 647+1001+1045+1045+

鍵入命令如下：

. stset time,failure（outcome）

得結果如下：

數據庫“例1”被定義為生存分析數據庫，變量“outcome”取值不等于0且不等于缺失值時，該記錄為完全數據，即出現預期結果。反之則為截尾值，表示未觀測到患者出現預期結果。完成上述定義后，即可用下面介紹的命令作進一步分析。

3.2   生存資料的描述

用于計算中位生存時間的命令是：

stsum [if表達式] [,by（分組變量）]

可用stci命令計算中位生存時間、平均生存時間、生存時間的百分位數，及其可信區間：

stci [if表達式],[by（分組變量）選擇項]

其中，選擇項有：

median /* 計算中位生存時間

emean /* 計算平均生存時間時，如果生存時間最長一例為截尾值，emean假設數據服從指數分布,并根據指數分布將該例后生存曲線部分延長至與橫軸相交，曲線下面積即為所求的平均生存時間。

rmean /* 計算平均生存時間時，如果生存時間最長一例為截尾值，rmean不對數據延長，曲線下面積即為所求的平均生存時間。此即為通常教科書上所教授的平均生存時間。

p（#）/* 生存時間的百分位數

level（#）/* 可信區間的可信度

也可用survsum命令計算中位生存時間的中位數。

繼續以例1數據為例，在命令窗口鍵入：

. stsum,by（treat）

數字化結果如下：

第2組（改進手術組）較早出現了截尾數據，故該組的中位生存時間無法進行估計，Stata用缺失值表示。

用stci命令可以計算平均生存時間及其可信區間：命令窗口鍵入命令如下：

. stci,rmean by（treat）

數字化結果如下：

第2組的平均生存時間明顯長于第1組。對于觀察隊列中最后1例為截尾值者，平均生存時間的估計值偏低(underestimated)。Stata在相應數值后加“*”表示。

3.3   生存率的估計

生存率的估計一般采用乘積極限（product-limit）法，又稱Kaplan-Meier法，其標準誤的計算用Greenwood近似法。根據生存率及其標準誤，可以繪制生存曲線，估計可信區間。

用于輸出生存率、生存率的標準誤等統計量的命令是：

sts list [if表達式] [,by（分組變量）strata（分層變量）adjustfor（校正變量）選擇項]

選擇項有：

failure /* 輸出死亡函數 (1-S(t))

na /* 輸出累積風險函數

level /* 規定所輸出可信區間的可信度

這里，by與strata選擇項的使用有所不同。使用by選擇項時，Stata對分組變量的不同水平分別計算生存函數和累積風險函數。而在使用strata選擇項時必須同時使用adjustfor選擇項，此時Stata將計算adjustfor選擇項中校正變量取值為0時的生存函數、累積風險函數，即計算基線生存函數、基線累積風險函數。

用于繪制Kaplan-Meier生存（死亡）曲線的命令是：

stsgraph [,by（分組變量）separate繪圖命令選擇項] /* 繪制Kaplan-Meier生存（死亡）曲線

stphplot [,by（分組變量）繪圖命令選擇項] /* 繪制log(-log（S（t）））與log(time）的線圖

stcoxkm [,by（分組變量）separate繪圖命令選擇項] /* 繪制Kaplan-Meier生存曲線與Cox預測曲線

在上述的“繪圖命令選擇項”中，可以選用xlab，ylab，xtick，ytick以及標題命令t1，t2，b1，b2，l1，l2，r1，r2等，但connect，symbol，pen等選擇項不能用，因y軸是從0開始，所以ylog等選擇項亦不能用。在stphplot命令中，可選用connect，symbol，pen等選擇項。

sts graph命令中的其他常用選擇項：

adjustfor（變量）/* 按指定變量進行調整

failure /* 指定繪制“死亡”曲線，缺失為繪制生存曲線

gwood /* 繪制生存（或死亡）曲線及其Greenwood可信區間

na /* 繪制Nelson-Aalen累積風險函數曲線

cna /* 繪制Nelson-Aalen累積風險函數曲線及其可信區間

lost /* 在曲線上標出該時間點截尾值例數

計算各組的生存率及標準誤，命令及結果如下：

. sts list if treat==1

繪制各組的生存曲線，命令及結果如下:

. sts graph,by（treat）lost

兩條曲線分別表示兩組的生存曲線，曲線上的數字表示在該時刻的截尾值例數。顯然，兩組的生存率不同。繪制各組的生存曲線及其可信區間，使用gwood選擇項。如對第1組，命令及結果如下：

. sts graph if treat==1,gwood

圖中，中間一條線是treat=1組的生存曲線，上、下兩條線分別表示生存率的可信區間的上下限。注意，率的可信區間是不對稱的。

3.4   生存率的比較

檢驗兩組或多組生存率是否相同一般采用Log-rank（Mantel-Cox）檢驗、Wilcoxon（Breslow）檢驗。命令如下：

sts test分組變量 [,選擇項]

選擇項有：

logrank/*進行Log-Rank檢驗

wilcoxon/*進行Wilcoxon（也叫Wilcoxon-Gehan檢驗或Breslow檢驗）

trend/*檢驗死亡）生存）率是否隨分組變量取值水平的增高而上升或下降

就例1資料，比較兩組患者的生存時間有無差別。鍵入命令如下：

. sts test treat,logrank

數字化結果如下：

這里的檢驗假設是第1處理組的生存率與第2組的相同。輸出結果中給出了兩組的實際數（events observed）及理論數（events expected）。本例中改進手術組的實際死亡數小于理論數，說明該組患者預后情況較好，經Log-rank檢驗，X²=6.71，自由度υ=1，P=0.0096，按α=0.05的檢驗水準認為兩組患者的生存時間有差別，以改進手術組為優。

3.5   Cox比例風險模型

惡性腫瘤患者生存時間的長短，不僅與治療有關，還受患者的年齡、性別、病情、心理、環境、社會等因素的影響，如果要確切地顯示治療措施的效果，所有的患者除了治療措施不同以外，其他影響因素必須相同（或相近），但這在實際上是不可能做到的。因此，我們最好能采用多因素分析方法，即分析包括治療措施在內的可能因素對生存時間長短的影響（大小和方向）。

但生存時間的分布往往不服從正態分布（大多為正偏態分布），有時不知道它的分布類型，又存在截尾數據（censored data），這樣，就不能用多元線性回歸方法來分析。而傳統的方法只能進行單因素分析，又不能利用截尾數據（censored data）。1972年，英國統計學家D. R. Cox提出了一種比例風險模型（Cox proportional hazard model），簡稱Cox模型。它可以分析多種因素對生存時間的影響，而且允許有“截尾”存在。是生存分析中最重要的模型之一。Cox模型主要用于腫瘤和其它慢性病的預后因素分析，也可以用于一般的臨床療效評價和隊列的病因探索。Cox比例風險模型的一般形式是：

$h\left({t,X} \right)={h_0}\left(t \right)\exp \left({{\beta _1}{X_1}+{\beta _2}{X_2}+ \cdots +{\beta _m}{X_m}} \right)$

h(t)表示研究對象在時點t暴露于協變量（x1,x2,…，xm）之下的風險函數，h0(t)表示所有協變量取值均為0時的基線風險函數。在Cox模型中h0(t)不能由樣本得出，因而不能估計生存率。但這并不妨礙對各協變量相對危險度的估計。

估計Cox比例風險模型的命令格式為：

stcox [協變量] [,nohr strata（分層變量）]

估計含有時依變量的Cox比例風險模型的命令格式為：

cox時間變量 [協變量] [,hr dead (截尾變量) strata（分層變量）選擇項]

進行逐步Cox回歸分析的命令為：

sw cox時間變量協變量[,cox命令選擇項逐步回歸選擇項]

這里選擇項有：

nohr /* 指定輸出回歸系數b而不是危險比exp（b）

tvc（時）依變量）/* 指定時依變量

texp（）/* 時依變量取值變化表達式

level（#）/* 可信區間的可信度

maximize-options /* 進行最大似然估計的控制選項

[注] 應用命令“cox”時無須事先應用stset對數據進行定義，且進行逐步回歸時只能使用cox命令。

用sw cox命令可以進行逐步Cox回歸分析。就例1資料進行Cox回歸分析。在應用stset對數據進行規定后，可直接用stcox命令進行Cox回歸分析。鍵入命令如下：

. stcox treat,nohr

或者也可以使用如下命令：

. cox time treat,dead（outcome)

結果如下：

或

風險函數一般用極大似然估計，用Newton-Rap-hson法迭代。結果中給出了每次迭代的似然函數之對數值（Log Likelihood），本例經4次迭代得極大似然估計變量treat的系數（coefficient）為-1.371 774，其標準誤（Std. Error）為.5708971，z值為-2.40，是對treat的系數是否為0的檢驗，P>|z|是相應的概率。本例的P=0.0016，故可認為treat的系數不為0。從而得Cox比例風險函數：

如果計算HR則可使用如下命令：

. stcox treat,hr

或者使用以下命令：

. cox time treat,dead(outcome) hr

結果如下：

或

以例2數據為例繼續演示Stata軟件實現Cox回歸。

例2 某臨床試驗比較A、B兩種治療方案對某病的治療效果，A組（group=0）12例，B組（group=1）13例。患者分組后檢驗其腎功能（kidney），功能正常者記0，不正常者記為1；治療后生存時間為stime（d）；數據已存入文件例2.dta。問不同治療方案及腎功能對患者的生存時間是否有影響？

這里，時間變量是stime，終檢變量是censor，治療方案（group）是研究因素，而腎功能（kidney）是混雜因素。例2數據如下圖所示：

鍵盤鍵入命令設置數據為生存數據，如下：

. stset stime,failure（censor）

結果如下：

繼續鍵入命令如下：

. stcox group kidney,nohr nolog

得結果如下：

計算HR，則輸入如下命令：

. stcox group kidney,hr

3.6   隨訪生存資料的壽命表法

當樣本含量較大或不能準確得知研究結果出現的時間時，可以將各研究對象的生存時間按年或月進行分組計算其生存率。Stata相應的命令是：

ltable時間變量 [dead(截尾變量)] [,by(分組變量) 選擇項繪圖選擇項]

ltable命令中大部分選擇項前面已經介紹過，未介紹過的有：

weight /* 指定權重變量

survival /* 指定輸出生存率

failure /* 指定輸出死亡率

hazard /* 指定輸出風險函數

interval /* 指定壽命表中生存時間組距

test /* 應用似然比檢驗、Log-rank檢驗對各總體生存率曲線是否相同進行檢驗

notab /* 不輸出壽命表

noconf /* 繪制生存率曲線時不繪制各時間點生存率的可信區間

graph /* 指定繪制生存率曲線

例3 隨訪某種惡性腫瘤患者生存情況如下圖所示，試作統計分析。這是一個分組資料，先將數據整理成下列形式，包括處理變量treat，生存年數year，是否截尾censor，以及頻數num。其中，生存年數輸入時“0～”輸為0.5，“1～”輸為1.5，其他依此類推。

計算壽命表并進行統計學檢驗，命令如下：

. ltable year censor [weight=num],test by(treat) interval(1)

得數字化結果如下：

Stata依次輸出各段生存時間起點及終點、期初人數、期內死亡人數、截尾例數、生存率及其標準誤和相應的95%可信區間。同時給出了兩組的齊性檢驗（Lawlsee，1982）及log-rank檢驗。

繪制第1組（group=1）患者的生存率曲線圖。命令如下：

. ltable year censor [weight=num],test by(treat) interval(1) graph

得數字化結果如下：

4   結語

生存分析應用廣泛，作為一個臨床醫生至少應該掌握使用一種統計學軟件實現生存分析，本文在參考了《現代醫學統計方法與Stata應用（第2版）》基礎上給大家演示了Stata軟件實現生存分析的過程，希望能對大家的科研工作有所幫助。

附作者簡介：周支瑞，復旦大學腫瘤醫院在讀博士。AME學術沙龍委員，愛好讀書，廣交友。業余時間擔任丁香園循證醫學討論版版主，自學臨床流行病學與循證醫學五年余。主要研究方向：惡性腫瘤的放射生物學研究、循證醫學與Meta分析方法學研究。目前以第一作者及共同作者發表SCI論文十余篇，參編學術著作兩部。

生存分析（survival analysis）即是將終點事件的出現與否和出現終點事件所經歷的時間結合起來的一種統計分析方法，生存分析通常研究的終點事件是死亡，生存分析由此得名。但生存分析可更廣泛地運用于惡性腫瘤、慢性疾病或其他情況的隨訪研究中的事件分析，比如疾病的發生、復發、轉移、傷口的愈合、某種癥狀的消失等。生存資料的分析主要特點就是考慮每個研究對象出現某一結局所經歷的時間。生存曲線即是以生存時間為橫軸，生存率為縱軸，將各個時間點對應的生存率連接在一起的曲線圖^[1-2]。

1   生存分析中幾個重要的基本概念

1.1   生存時間

生存時間（survival time）也是一個廣義概念，泛指所關心的某現象的持續時間，即隨訪觀察持續的時間，常用符號t表示。生存時間分為兩種類型：（1）完全數據（complete data）：指從觀察起點到發生“死亡”事件所經歷的時間。提供了觀察對象確切的生存時間。（2）截尾數據（censored data）：亦稱截尾值（censored value）或終檢值。指從觀察起點到發生非“死亡”事件所經歷的時間。

1.2   截尾

生存結局分為“死亡”與“截尾”兩類，“死亡”是感興趣的終點時間，其他終點事件或結局都歸為截尾。

1.3   死亡概率

死亡概率（probability of death）表示單位時間段開始存活的個體，在該段時間內死亡的可能性。符號q表示。“q=某年內死亡人數÷某年年初人口數”。

1.4   生存概率

生存概率（probability of survival）表示單位時間段開始存活的個體，到該段時間結束時仍存活的可能性。符號p表示。p=某年活滿一年人口數÷某年年初人口數。P=1-q。

1.5   生存率

生存率（survival rate,survival function）表示觀察對象經歷t_k個單位時間段后仍存活的可能性。若無截尾數據，則S（t_k）=P（T>t_k）=t_k。其中0<S（t）≤1。若有截尾數據，須分時段計算生存概率。假定觀察對象在各個時段的生存事件獨立，應用概率乘法定理：S（t_k）=P（T>t_k）=P₁·P₂…P_k^！，P_i為某時段的生存概率，故生存率又稱累積生存概率（cumulative probability of survival）。

1.6   生存曲線

生存曲線（survival curve）：生存時間為橫軸，將各時點所對應的生存率連接在一起的曲線圖，樣本量小時生存曲線呈階梯形，樣本量足夠大時，形成光滑的曲線。

1.7   中位生存時間

中位生存時間是指50%觀察對象能存活的時間。

2   生存分析的統計學方法

由于生存時間一般不呈正態分布，而且需要考慮截尾數據，生存分析有其獨特的統計學方法。常用的統計學方法有以下幾種。

2.1   描述性分析

根據樣本生存資料估計總體生存率及其他有關指標（如中位生存時間等）。常采用Kaplan-Meier法（乘積極限法）進行分析。對于頻數表資料則采用壽命表法進行分析。計算生存率需要考慮時間順序。

2.2   比較分析的方法

對不同組生存率進行比較分析，常采用非參數的log-rank檢驗，檢驗無效假設使兩組或多組總體生存時間分布相同。

2.3   影響因素分析

通過生存分析模型來探討影響生存時間的因素，常用的方法為COX比例風險模型。

3   基于Stata軟件的統計學實現生存分析

（筆者注：以下所舉實例數據全部來自于陳峰教授主編《現代醫學統計方法與Stata應用（第2版）》，相關Stata命令及結果解釋大部分來自于這本書，其中部分命令有少許改動。陳鋒教授主編的這本書通俗易懂，感興趣的讀者可以找來一讀。）

3.1   生存資料的定義

在對隨防資料進行生存分析之前，需先將該數據庫定義為生存資料數據庫，其命令是：

stset時間變量[,failure (截尾變量==#)]

其中，選擇項failure（截尾變量==#）規定截尾變量取值為“#”時研究對象出現預期結果，沒有該選擇項時，Stata以所有不等于0 的非缺失值為出現預期結果。對數據庫進行定義時必須注意變量順序，命令stset后的變量順序依次為時間變量、截尾變量。定義數據庫后，系統自動產生四個變量：

_st /* 數據庫中該條記錄是否被定義為生存資料

_d /* 該條記錄是否出現預期結果

_t /* 觀察對象被隨訪的時間

_t0 /* 觀察對象第一次被觀察到的時間（開始過程的時間為0）

例1 某醫院泌尿外科于1979～1982年間施行19例腎移植手術，擬了解腎移植后患者的生存時間（d）。規定隨訪開始時間為患者術后1 d，預期結果為該患者因與腎移植有關的各種原因的死亡。后改進手術方式，于1983～1986 年又施行14例，資料如下（有+的數據表示該患者截尾）。計算各組的生存率及可信區間（資料已存入文件例1.dta）。

一般手術: 3 9 15 20 20 26 30 41 46 64+64 135 223 365 450 596+680+900+900+

改進手術: 10 390+518+70+70+120 475+225 801+366 647+1001+1045+1045+

鍵入命令如下：

. stset time,failure（outcome）

得結果如下：

數據庫“例1”被定義為生存分析數據庫，變量“outcome”取值不等于0且不等于缺失值時，該記錄為完全數據，即出現預期結果。反之則為截尾值，表示未觀測到患者出現預期結果。完成上述定義后，即可用下面介紹的命令作進一步分析。

3.2   生存資料的描述

用于計算中位生存時間的命令是：

stsum [if表達式] [,by（分組變量）]

可用stci命令計算中位生存時間、平均生存時間、生存時間的百分位數，及其可信區間：

stci [if表達式],[by（分組變量）選擇項]

其中，選擇項有：

median /* 計算中位生存時間

emean /* 計算平均生存時間時，如果生存時間最長一例為截尾值，emean假設數據服從指數分布,并根據指數分布將該例后生存曲線部分延長至與橫軸相交，曲線下面積即為所求的平均生存時間。

rmean /* 計算平均生存時間時，如果生存時間最長一例為截尾值，rmean不對數據延長，曲線下面積即為所求的平均生存時間。此即為通常教科書上所教授的平均生存時間。

p（#）/* 生存時間的百分位數

level（#）/* 可信區間的可信度

也可用survsum命令計算中位生存時間的中位數。

繼續以例1數據為例，在命令窗口鍵入：

. stsum,by（treat）

數字化結果如下：

第2組（改進手術組）較早出現了截尾數據，故該組的中位生存時間無法進行估計，Stata用缺失值表示。

用stci命令可以計算平均生存時間及其可信區間：命令窗口鍵入命令如下：

. stci,rmean by（treat）

數字化結果如下：

第2組的平均生存時間明顯長于第1組。對于觀察隊列中最后1例為截尾值者，平均生存時間的估計值偏低(underestimated)。Stata在相應數值后加“*”表示。

3.3   生存率的估計

生存率的估計一般采用乘積極限（product-limit）法，又稱Kaplan-Meier法，其標準誤的計算用Greenwood近似法。根據生存率及其標準誤，可以繪制生存曲線，估計可信區間。

用于輸出生存率、生存率的標準誤等統計量的命令是：

sts list [if表達式] [,by（分組變量）strata（分層變量）adjustfor（校正變量）選擇項]

選擇項有：

failure /* 輸出死亡函數 (1-S(t))

na /* 輸出累積風險函數

level /* 規定所輸出可信區間的可信度

這里，by與strata選擇項的使用有所不同。使用by選擇項時，Stata對分組變量的不同水平分別計算生存函數和累積風險函數。而在使用strata選擇項時必須同時使用adjustfor選擇項，此時Stata將計算adjustfor選擇項中校正變量取值為0時的生存函數、累積風險函數，即計算基線生存函數、基線累積風險函數。

用于繪制Kaplan-Meier生存（死亡）曲線的命令是：

stsgraph [,by（分組變量）separate繪圖命令選擇項] /* 繪制Kaplan-Meier生存（死亡）曲線

stphplot [,by（分組變量）繪圖命令選擇項] /* 繪制log(-log（S（t）））與log(time）的線圖

stcoxkm [,by（分組變量）separate繪圖命令選擇項] /* 繪制Kaplan-Meier生存曲線與Cox預測曲線

在上述的“繪圖命令選擇項”中，可以選用xlab，ylab，xtick，ytick以及標題命令t1，t2，b1，b2，l1，l2，r1，r2等，但connect，symbol，pen等選擇項不能用，因y軸是從0開始，所以ylog等選擇項亦不能用。在stphplot命令中，可選用connect，symbol，pen等選擇項。

sts graph命令中的其他常用選擇項：

adjustfor（變量）/* 按指定變量進行調整

failure /* 指定繪制“死亡”曲線，缺失為繪制生存曲線

gwood /* 繪制生存（或死亡）曲線及其Greenwood可信區間

na /* 繪制Nelson-Aalen累積風險函數曲線

cna /* 繪制Nelson-Aalen累積風險函數曲線及其可信區間

lost /* 在曲線上標出該時間點截尾值例數

計算各組的生存率及標準誤，命令及結果如下：

. sts list if treat==1

繪制各組的生存曲線，命令及結果如下:

. sts graph,by（treat）lost

兩條曲線分別表示兩組的生存曲線，曲線上的數字表示在該時刻的截尾值例數。顯然，兩組的生存率不同。繪制各組的生存曲線及其可信區間，使用gwood選擇項。如對第1組，命令及結果如下：

. sts graph if treat==1,gwood

圖中，中間一條線是treat=1組的生存曲線，上、下兩條線分別表示生存率的可信區間的上下限。注意，率的可信區間是不對稱的。

3.4   生存率的比較

檢驗兩組或多組生存率是否相同一般采用Log-rank（Mantel-Cox）檢驗、Wilcoxon（Breslow）檢驗。命令如下：

sts test分組變量 [,選擇項]

選擇項有：

logrank/*進行Log-Rank檢驗

wilcoxon/*進行Wilcoxon（也叫Wilcoxon-Gehan檢驗或Breslow檢驗）

trend/*檢驗死亡）生存）率是否隨分組變量取值水平的增高而上升或下降

就例1資料，比較兩組患者的生存時間有無差別。鍵入命令如下：

. sts test treat,logrank

數字化結果如下：

這里的檢驗假設是第1處理組的生存率與第2組的相同。輸出結果中給出了兩組的實際數（events observed）及理論數（events expected）。本例中改進手術組的實際死亡數小于理論數，說明該組患者預后情況較好，經Log-rank檢驗，X²=6.71，自由度υ=1，P=0.0096，按α=0.05的檢驗水準認為兩組患者的生存時間有差別，以改進手術組為優。

3.5   Cox比例風險模型

惡性腫瘤患者生存時間的長短，不僅與治療有關，還受患者的年齡、性別、病情、心理、環境、社會等因素的影響，如果要確切地顯示治療措施的效果，所有的患者除了治療措施不同以外，其他影響因素必須相同（或相近），但這在實際上是不可能做到的。因此，我們最好能采用多因素分析方法，即分析包括治療措施在內的可能因素對生存時間長短的影響（大小和方向）。

但生存時間的分布往往不服從正態分布（大多為正偏態分布），有時不知道它的分布類型，又存在截尾數據（censored data），這樣，就不能用多元線性回歸方法來分析。而傳統的方法只能進行單因素分析，又不能利用截尾數據（censored data）。1972年，英國統計學家D. R. Cox提出了一種比例風險模型（Cox proportional hazard model），簡稱Cox模型。它可以分析多種因素對生存時間的影響，而且允許有“截尾”存在。是生存分析中最重要的模型之一。Cox模型主要用于腫瘤和其它慢性病的預后因素分析，也可以用于一般的臨床療效評價和隊列的病因探索。Cox比例風險模型的一般形式是：

$h\left({t,X} \right)={h_0}\left(t \right)\exp \left({{\beta _1}{X_1}+{\beta _2}{X_2}+ \cdots +{\beta _m}{X_m}} \right)$

h(t)表示研究對象在時點t暴露于協變量（x1,x2,…，xm）之下的風險函數，h0(t)表示所有協變量取值均為0時的基線風險函數。在Cox模型中h0(t)不能由樣本得出，因而不能估計生存率。但這并不妨礙對各協變量相對危險度的估計。

估計Cox比例風險模型的命令格式為：

stcox [協變量] [,nohr strata（分層變量）]

估計含有時依變量的Cox比例風險模型的命令格式為：

cox時間變量 [協變量] [,hr dead (截尾變量) strata（分層變量）選擇項]

進行逐步Cox回歸分析的命令為：

sw cox時間變量協變量[,cox命令選擇項逐步回歸選擇項]

這里選擇項有：

nohr /* 指定輸出回歸系數b而不是危險比exp（b）

tvc（時）依變量）/* 指定時依變量

texp（）/* 時依變量取值變化表達式

level（#）/* 可信區間的可信度

maximize-options /* 進行最大似然估計的控制選項

[注] 應用命令“cox”時無須事先應用stset對數據進行定義，且進行逐步回歸時只能使用cox命令。

用sw cox命令可以進行逐步Cox回歸分析。就例1資料進行Cox回歸分析。在應用stset對數據進行規定后，可直接用stcox命令進行Cox回歸分析。鍵入命令如下：

. stcox treat,nohr

或者也可以使用如下命令：

. cox time treat,dead（outcome)

結果如下：

或

風險函數一般用極大似然估計，用Newton-Rap-hson法迭代。結果中給出了每次迭代的似然函數之對數值（Log Likelihood），本例經4次迭代得極大似然估計變量treat的系數（coefficient）為-1.371 774，其標準誤（Std. Error）為.5708971，z值為-2.40，是對treat的系數是否為0的檢驗，P>|z|是相應的概率。本例的P=0.0016，故可認為treat的系數不為0。從而得Cox比例風險函數：

如果計算HR則可使用如下命令：

. stcox treat,hr

或者使用以下命令：

. cox time treat,dead(outcome) hr

結果如下：

或

以例2數據為例繼續演示Stata軟件實現Cox回歸。

例2 某臨床試驗比較A、B兩種治療方案對某病的治療效果，A組（group=0）12例，B組（group=1）13例。患者分組后檢驗其腎功能（kidney），功能正常者記0，不正常者記為1；治療后生存時間為stime（d）；數據已存入文件例2.dta。問不同治療方案及腎功能對患者的生存時間是否有影響？

這里，時間變量是stime，終檢變量是censor，治療方案（group）是研究因素，而腎功能（kidney）是混雜因素。例2數據如下圖所示：

鍵盤鍵入命令設置數據為生存數據，如下：

. stset stime,failure（censor）

結果如下：

繼續鍵入命令如下：

. stcox group kidney,nohr nolog

得結果如下：

計算HR，則輸入如下命令：

. stcox group kidney,hr

3.6   隨訪生存資料的壽命表法

當樣本含量較大或不能準確得知研究結果出現的時間時，可以將各研究對象的生存時間按年或月進行分組計算其生存率。Stata相應的命令是：

ltable時間變量 [dead(截尾變量)] [,by(分組變量) 選擇項繪圖選擇項]

ltable命令中大部分選擇項前面已經介紹過，未介紹過的有：

weight /* 指定權重變量

survival /* 指定輸出生存率

failure /* 指定輸出死亡率

hazard /* 指定輸出風險函數

interval /* 指定壽命表中生存時間組距

test /* 應用似然比檢驗、Log-rank檢驗對各總體生存率曲線是否相同進行檢驗

notab /* 不輸出壽命表

noconf /* 繪制生存率曲線時不繪制各時間點生存率的可信區間

graph /* 指定繪制生存率曲線

例3 隨訪某種惡性腫瘤患者生存情況如下圖所示，試作統計分析。這是一個分組資料，先將數據整理成下列形式，包括處理變量treat，生存年數year，是否截尾censor，以及頻數num。其中，生存年數輸入時“0～”輸為0.5，“1～”輸為1.5，其他依此類推。

計算壽命表并進行統計學檢驗，命令如下：

. ltable year censor [weight=num],test by(treat) interval(1)

得數字化結果如下：

Stata依次輸出各段生存時間起點及終點、期初人數、期內死亡人數、截尾例數、生存率及其標準誤和相應的95%可信區間。同時給出了兩組的齊性檢驗（Lawlsee，1982）及log-rank檢驗。

繪制第1組（group=1）患者的生存率曲線圖。命令如下：

. ltable year censor [weight=num],test by(treat) interval(1) graph

得數字化結果如下：

4   結語

生存分析應用廣泛，作為一個臨床醫生至少應該掌握使用一種統計學軟件實現生存分析，本文在參考了《現代醫學統計方法與Stata應用（第2版）》基礎上給大家演示了Stata軟件實現生存分析的過程，希望能對大家的科研工作有所幫助。

附作者簡介：周支瑞，復旦大學腫瘤醫院在讀博士。AME學術沙龍委員，愛好讀書，廣交友。業余時間擔任丁香園循證醫學討論版版主，自學臨床流行病學與循證醫學五年余。主要研究方向：惡性腫瘤的放射生物學研究、循證醫學與Meta分析方法學研究。目前以第一作者及共同作者發表SCI論文十余篇，參編學術著作兩部。