生存分析完整版本

第１章基本概念

第１節(jié)生存資料的特點

生存資料(SurvivalData)或失效時間資料(Failure-timeData)與多元線性回歸資料很相似，只不過因變量(或反應(yīng)變量)通常為觀測對象生存的時間，常用ｔ來表示。當(dāng)然，生存時間是廣義的，可以指在通常意義下生物體的生存時間、也可以指所關(guān)心的某現(xiàn)象（如疾病治愈后、合格品使用后）持續(xù)的時間。若生存時間是準(zhǔn)確觀測到的，則稱為完全數(shù)據(jù)。生存資料的一個明顯特點是：所收集的資料中常常包含不完全數(shù)據(jù)，也稱為截尾數(shù)據(jù)、刪失數(shù)據(jù)、終檢數(shù)據(jù)（CensoredData）。包括刪失數(shù)據(jù)的資料，稱為刪失資料。對于刪失數(shù)據(jù)，既不能簡單地棄之，踴能像對待完全數(shù)據(jù)那樣給予充分的信任，需要采取一些技術(shù)處理。專門處理這種資料的統(tǒng)計方法，稱為生存分析（SurvivalAnalysis）。

導(dǎo)致數(shù)據(jù)刪失有多種原因，最常見的有：失訪（病人因搬家、隨訪信件丟失、車禍等原因，導(dǎo)致醫(yī)生對他們的隨訪觀察中斷）和研究截止。由隨機因素引起的，稱為隨機刪失；若事先就定了截止日期，則稱為定時刪失（也稱Ⅰ型刪失）；若事先就定了觀察完多少例就截止研究，則稱為Ⅱ型刪失（也稱為定數(shù)刪失）。在表達刪失數(shù)據(jù)時，常在其右上角放一個“＋”號；而用ＳＡＳ軟件分析時，常在其前放一個“－”號或產(chǎn)生１個指示變量（如：Ｃ=０表示刪失數(shù)據(jù)、Ｃ=１表示完全數(shù)據(jù)，反過來也可以），便于計算時區(qū)別對待。為了使數(shù)據(jù)的表達與計算在形式上統(tǒng)一起來，本篇一律用負(fù)數(shù)表示刪失數(shù)據(jù)，因生存時間不可能為負(fù)值，故不會產(chǎn)生混淆。

第２節(jié)生存時間函數(shù)

描述生存時間規(guī)律的函數(shù)很多，統(tǒng)稱為生存時間函數(shù)。其中最主要的有生存函數(shù)、死亡概率函數(shù)、概率密度函數(shù)和危險率函數(shù)。

１．生存函數(shù)(SurvivalFunction)

生存函數(shù)也稱為生存概率或累積生存率，常用S(t)表示，它表示一個體生存時間長于ｔ的概率。在具體問題中，該函數(shù)在ｔ時刻的取值可用式(5.1.1)來估計∶

S(t)≈生存時間長于ｔ的病人數(shù)／病人總數(shù)(5.1.1)

２．死亡概率函數(shù)(FailureProbabilityFunction)

死亡概率函數(shù)簡稱為死亡概率，常用F(t)表示，它表示一個體從開始觀察起到時刻ｔ為止的死亡概率。它可以通過S(t)求得(詳后)。

３．概率密度函數(shù)(ProbabilityDensityFunction)

概率密度函數(shù)簡稱為密度函數(shù)，常用f(t)表示，它表示一個體死于(t，t+△t)小區(qū)間內(nèi)的概率的極限。在具體問題中，該函數(shù)在ｔ時刻的取值可用式(5.1.2)來估計∶

f(t)≈t時刻開始的區(qū)間內(nèi)死汀人數(shù)／(病人總數(shù)×區(qū)間寬度)(5.1.2)

４．危險率函數(shù)(HazardFunction)

危險率函數(shù)也稱為風(fēng)險函數(shù)、瞬時死亡率、年齡別死亡率、條件死亡率,常用h(t)表示,它表示已存活到t的一個體,死于(t,t+△t)小區(qū)間內(nèi)的概率的極限。在具體問題中，該函數(shù)在ｔ時刻的取值可用式(5.1.3)來估計∶

h(t)≈ｔ時刻開始的區(qū)間內(nèi)死汀人數(shù)／（生存到ｔ的病人數(shù)×區(qū)間寬度）(5.1.3)

５．上述幾個函數(shù)之間的相互關(guān)系(5.1.4)(5.1.5)(5.1.6)(5.1.7)

上述各函數(shù)中“＇”代表對ｔ求導(dǎo)數(shù)，“∫”代表積分。

第３節(jié)生存分析方法的分類

像普通統(tǒng)計分析一樣，生存分析也有一套完整的方法：統(tǒng)計描述（包括求生存時間的分位數(shù)、中數(shù)生存期、平均數(shù)、生存函數(shù)的估計、判斷生存時間的圖示法）；非參數(shù)檢驗（檢驗分組變量各水平所對應(yīng)的生存曲線是否一致，常用的方法有對數(shù)秩檢驗（Log-rankTest）、威爾科克森檢驗（WilcoxonTest）和似然比檢驗（LikelihoodratioTest））；ＣＯＸ模型(半?yún)?shù)模型)回歸分析（在特定的假設(shè)之下，建立生存時間隨多個危險因素變化的回歸方程）；參數(shù)模型回歸分析（已知生存時間服從特定的參數(shù)模型時，擬合相應(yīng)的參數(shù)模型，更準(zhǔn)確地刻劃變量之間的變化規(guī)律）。

第２章生存資料的非參數(shù)統(tǒng)計方法

第１節(jié)統(tǒng)計描述與非參數(shù)分析概述

１．統(tǒng)計描述

常用來反映一組生存時間平均水平的統(tǒng)計指標(biāo)有中位數(shù)、平均數(shù)２種，因生存資料多為正偏態(tài),故往往選用中位數(shù)更符合資料的特點。

對于壽命資料,首先需給出各時間點上生存函數(shù)的估計值，常用的方法有：乘積─極限法（Product-LimitMethod，簡稱PL法）和壽命表法（Life-TableMethod，簡稱LT法）。PL法是利用ti時刻之前各時間點上生存率的乘積來估計在時刻ti的生存函數(shù)S(ti)、而LT表法是通過計數(shù)落入時間區(qū)間[ti-1，ti]內(nèi)的失效和刪失的觀察例數(shù)來估計S(ti)。

若能知道壽命函數(shù)的具體，可有的放矢地去選用相應(yīng)的參數(shù)模型擬合資料，是非常有益的。實現(xiàn)這一目的途經(jīng)是圖解法，如:用(t,-logS(t))畫圖，若成一條直線，表明S(t)呈指數(shù)；又如：用(logt,log(-logS(t))畫圖，若成一條直線，表明S(t)呈圖爾。當(dāng)然，也有一些統(tǒng)計檢驗方法，如：判斷是否服從指數(shù)的Ｇ檢驗法、判斷是否服從圖爾的Mann-Scheuer-FertigTiku檢驗法和判斷是否服從對數(shù)正態(tài)的Ｗ檢驗法等，具體檢驗方法參見有關(guān)專著。

２．各層間生存曲線的齊性檢驗

設(shè)全部受試者接受了ｋ只同的處理，這ｋ種處理實際上就是一個名義分類變量或楓因素的ｋ個水平，于是，可按層估計生存函數(shù)。研究者常需比較ｋ條生存曲線之間是否有顯著差別，其方法有多種，ＳＡＳ中用了以下３種：對數(shù)秩檢驗（Log-rankTest）、威爾科克森檢驗（WilcoxonTest）和似然比檢驗（LikelihoodratioTest）。用它們來實現(xiàn)各層之間的齊性檢驗。

３．上述３種非參數(shù)檢驗的比較

當(dāng)生存時間的為指數(shù)、圖爾或?qū)儆诒壤ｋU模型時，Log-rank檢驗效率較高；當(dāng)生存時間的為對數(shù)正態(tài)等時，Wilcoxon檢驗效率較高；似然比檢驗是建立在指數(shù)模型上的，故當(dāng)資料偏離此模型時，其結(jié)果不如前２種檢驗方法穩(wěn)健。

４．協(xié)變量與生存時間聯(lián)系密切程度的檢驗

當(dāng)資料中還包含與生存時間有關(guān)的其他連續(xù)變量（即協(xié)變量）時,也可分析它們與生存時間聯(lián)系的密切程度。為實現(xiàn)此檢驗，LIFETEST過程中提供了２個分別建立在指數(shù)得分和威爾科克森得分基礎(chǔ)之上的刪失數(shù)據(jù)線性秩統(tǒng)計量─Log-rankTest和WilcoxonTest，這２種檢驗通過合并楓變量后進行計算，從而，校正了楓變量的影響。除了對重復(fù)（ties）生存時間的處理方法不同外，這里所說的２種檢驗與實現(xiàn)各層之間齊性檢驗中所提到的前２種檢驗是相同的。

為了不把讀者的注意力引向復(fù)雜的計算，特將上述各種方法的具體計算公逝在本章第３節(jié)中再介紹，以便必要時備查。

第２節(jié)用LIFETEST過程實現(xiàn)統(tǒng)計計算

[例5.2.1]某醫(yī)生收集到35例白血病患者治療后的生存時間ｔ(月),仔細觀察后發(fā)現(xiàn)這些病人中有一部分人出現(xiàn)了白細胞(WBC)倍增的現(xiàn)象?，F(xiàn)將他們按是否出現(xiàn)WBC倍增分成２組如下(注:負(fù)值代表刪失數(shù)據(jù))，試用生存分析方法分析患者有無WBC倍增，對其生存時間長短有無顯著影響。

A組(有WBC倍增)：2,-2.5,3.5,4,4,-5,6,-6,7,-7,8,-9,10.5,12.5,19;

B組(無WBC倍增)：2.5,5,7,-8.5,9,-10,11,-11,12,13,-14,15,-16,17,-18,19,-20,21,

24,32。

[SAS程序]──[D5P1.PRG]

DATAabc;PROCLIFETESTMETHOD=PL

INFILE'a:hlwbc.dat';PLOTS=(S,LS,LLS);

INPUTlt@@;TIMEt*censor(1);

IFlt<0THENcensor=1;STRATAgroup;

ELSEcensor=0;RUN;

IF_N_<16THENgroup='high-wbc';PROCLIFETESTMETHOD=LIFE

ELSEgroup='low-wbc';PLOTS=(S,H);

t=ABS(lt);TIMEt*censor(1);

STRATAgroup;

RUN;

(程序的第１部分)(程序的第２部分)

[程序修改指導(dǎo)]用全部35個數(shù)據(jù)建立的數(shù)據(jù)文件名為HLWBC.DAT，第１個IF語句產(chǎn)生１個指示變量CENSOR，其取值為１時為刪失數(shù)據(jù)、取值為０時為完全數(shù)據(jù)。第２個IF語句產(chǎn)生１個分組變量GROUP，前15個數(shù)據(jù)屬于有WBC倍增組、后20個數(shù)據(jù)屬于無WBC倍增組。對表示刪失和完全數(shù)據(jù)的變量lt取絕對值是為了保證參與計算的生存時間ｔ都是正值。

第１個過程步是選擇PL法計算（它也是隱含的方法）、第２個過程步是選擇LT法計算。PLOTS=要求繪圖，其中S表示生存函數(shù)、L表示取對數(shù)、H表示危險率函數(shù)，圖形的橫坐標(biāo)與縱坐標(biāo)分別為：

S─(t,S)、LS─(t,-log(S))、LLS─(log(t),log(-(log(S))))、H─(t,H)

生存時間ｔ與指示變量以乘法的形式寫在TIME語句中、分組變量寫在STRATA語句中。

當(dāng)用壽命表(LT)法分析資料時，程序會自動形成生存時間的區(qū)間，也可人為指定生存時間的分組區(qū)間。做法是：在PROC語句的分號之前加上INTERVALS=(atobbyc)，a、b、c分別為初值、終值、步長（必須是具體數(shù)值），步長的缺省值為１。

如果資料中還含有數(shù)值型的協(xié)變量，可將它們寫在TEST語句中，如：TESTx1x2x3;以便檢驗協(xié)變量與生存時間聯(lián)系的密切程度。當(dāng)然，若有PHREG和LIFEREG過程，用它們建立起因變量ｔ隨自變量（即危險因素）變化的回歸模型，可更好地揭示變量之間的內(nèi)在聯(lián)系。

[輸出結(jié)果及其解釋]

Product-LimitSurvivalEstimates

GROUP=high-wbc

①②③④⑤⑥

Survival

StandardNumberNumber

TSurvivalFailureErrorFailedLeft

0.00001.000000015

2.00000.93330.06670.0644114

2.5000*...113

3.50000.86150.13850.0911212

4.0000...311

4.00000.71790.28210.1198410

5.0000*...49

6.00000.63820.36180.130458

6.0000*...57

7.00000.54700.45300.140066

7.0000*...65

8.00000.43760.56240.148774

9.0000*...73

10.50000.29170.70830.155082

12.50000.14590.85410.129091

19.000001.00000100

*CensoredObservation

Quantiles75%12.5000Mean9.0775

50%8.0000StandardError1.6768

25%4.0000

這是用PL法對第１組生存資料進行統(tǒng)計描述的結(jié)果。標(biāo)號①～⑥分別是生存時間、生存概率、死亡概率、生存概率的標(biāo)準(zhǔn)誤差、已觀察到的不同失效時間的數(shù)目、尚未觀察到的不同失效或刪失時間的數(shù)目，打＊號的是刪失觀察值。接著，給出了生存時間的四分位數(shù)、均數(shù)及其標(biāo)準(zhǔn)誤差。結(jié)果顯示∶第１組患者中有25％的人(約４人)的生存時間短于４個月，即有75％的人的生存時間長于４個月；同理,可解釋P50=8(個月)、P75=12.5(個月)的含義。由此可知∶該組患者的中數(shù)生存期為８個月、平均生存期約為９個月。

Product-LimitSurvivalEstimates

GROUP=low-wbc

Survival

StandardNumberNumber

TSurvivalFailureErrorFailedLeft

0.00001.000000020

2.50000.95000.05000.0487119

5.00000.90000.10000.0671218

7.00000.85000.15000.0798317

8.5000*...316

9.00000.79690.20310.0908415

10.0000*...414

11.00000.74000.26000.1006513

11.0000*...512

12.00000.67830.32170.1095611

13.00000.61660.38340.1156710

14.0000*...79

15.00000.54810.45190.121488

16.0000*...87

17.00000.46980.53020.126896

18.0000*...95

19.00000.37590.62410.1317104

20.0000*...103

21.00000.25060.74940.1348112

24.00000.12530.87470.1113121

32.000001.00000130

*CensoredObservation

Quantiles75%24.0000Mean17.1618

50%17.0000StandardError2.2053

25%11.0000

這是用PL法對第２組生存資料的統(tǒng)計描述結(jié)果。各列的解釋同上，從略。第２組患者的

中數(shù)生存期為１７個月、平均生存期約為１７個月。

SummaryoftheNumberofCensoredandUncensoredValues

GROUPTotalFailedCensored?nsored

high-wbc1510533.3333

low-wbc2013735.0000

Total35231234.2857

這是２組患者的總?cè)藬?shù)、死亡數(shù)、刪失數(shù)和刪失百分比。

SurvivalFunctionEstimates

SSDF|

u1.0+*---HL

r|H*-HH-L---L

v||L---L---L

i|H---HL-L-L

v|H-HL---L

a0.5+H-HL---L---L

l|H----HL---L

|H---HL-----L

D|||

i|H------------HL---------------L

s0.0+HL

t---+----+----+----+----+----+----+----+----+----+----+----+----+----+---

r0.02.55.07.510.012.515.017.520.022.525.027.530.032.5

i

T

這是反映２組患者生存情況的生存曲線圖，H表示有WBC倍增組、L表示無WBC倍增組，從圖上可明顯看出：無WBC倍增患者比有WBC倍增患者的生存期長。

CensoredObservations

Strata

L+LLLLLLL

H+HHHHH

-------+------+------+------+------+------+------+------+-------

05101520253035

T

這幅圖反映了各組患者刪失時間的情況。

-Log(SurvivalFunction)Estimates

N-LOGSDF|

e2+H+L

g|+++

a|+++

t|++L

i|+H++

v1+++++L

e|+H+++L+

|+H++L+

L|H+++H++L+L+L+

o|+H++L+++L+++L+

g0+*+++H*++

-+----+----+----+----+----+----+----+----+----+----+----+----+----+

S0.02.55.07.510.012.515.017.520.022.525.027.530.032.5

D

T

這是按(t,-log(S(t))繪出的圖，２條線都不呈直線趨勢，說明生存時間不呈指數(shù)。

Log(-Log(SurvivalFunction))Estimates

LL(-L(S))|

o2+

g|

|++H+L

N0++H++++H++L+L+L+

e|++++H++H++LL++L+

g|+H+++++L+++L

a-2++++++H++++++L++++++L++

t|H+++*++++++

i|

v-4+

e|

---+----+----+----+----+----+----+----+----+----+----+----+----+---

L0.500.751.001.251.501.752.002.252.502.753.003.253.50

o

LogT

LegendforStrataSymbols

H:GROUP=high-wbcL:GROUP=low-wbc

這是按(logt,log(-log(S(t)))繪出的圖，２條現(xiàn)別近似呈直線趨勢，說明生存時間近似呈圖爾。

TestingHomogeneityofSurvivalCurvesoverStrata

RankStatistics

①GROUPLog-RankWilcoxon

high-wbc5.073946117

low-wbc-5.07395-117

②CovarianceMatrixfortheLog-RankStatistics

GROUPhigh-wbclow-wbc

high-wbc3.36249-3.36249

low-wbc-3.362493.36249

③CovarianceMatrixfortheWilcoxonStatistics

GROUPhigh-wbclow-wbc

high-wbc2161.30-2161.30

low-wbc-2161.302161.30

④TestofEqualityoverStrata

Pr>

TestChi-SquareDFChi-Square

Log-Rank7.656510.0057

Wilcoxon6.333710.0118

-2Log(LR)2.834710.0922

這是關(guān)于各層生存曲線之間齊性檢驗的結(jié)果。①用２種檢驗方法算得公式(5.2.20)中的向量v，即v=(5.073946,-5.07395)'(Log-Rank法)、v=(117,-117)'(Wilcoxon法)。②、③分別用２種檢驗法算得此式中的協(xié)方差矩陣Ｖ，它們都是計算④中卡方值的中間結(jié)果，讀者最需要的是第④部分。這里給出了３種檢驗法的檢驗結(jié)果∶P值依次為0.0057(Log-Rank法)、0.0118(Wilcoxon法)、0.0922(似然比檢驗法，-2Log(LR))。

LifeTableSurvivalEstimates

GROUP=high-wbc

⑴⑵⑶⑷⑸⑹

Conditional

EffectiveConditionalProbability

IntervalNumberNumberSampleProbabilityStandard

[Lower,Upper)FailedCensoredSizeofFailureError

054114.50.27590.1174

510348.00.37500.1712

1015203.00.66670.2722

1520101.01.00000

⑺⑻⑼⑽⑾

SurvivalMedianMedian

IntervalStandardResidualStandard

[Lower,Upper)SurvivalFailureErrorLifetimeError

051.0000009.12702.4177

5100.72410.27590.11746.50002.1213

10150.45260.54740.14403.75002.1651

15200.15090.84910.1322..

EvaluatedattheMidpointoftheInterval

⑿⒀⒁⒂

PDFHazard

IntervalStandardStandard

[Lower,Upper)PDFErrorHazardError

050.05520.02350.0640.031588

5100.05430.02630.0923080.051855

10150.06030.03120.20.122474

15200.03020.02640.40

這是用LT法對第１組資料進行統(tǒng)計描述的結(jié)果。編號⑴～⒂所代表的含義分別為：⑴按區(qū)間寬度=5將生存時間自動劃分成若干區(qū)間、⑵死亡數(shù)、⑶刪失數(shù)、⑷有效樣本大小、⑸死亡的條件概率、⑹第⑸列數(shù)據(jù)的標(biāo)準(zhǔn)誤差、⑺區(qū)間左端點處生存概率、⑻區(qū)間左端點處死亡概率、⑼第⑺列數(shù)據(jù)的標(biāo)準(zhǔn)誤差、⑽中數(shù)剩余生存壽命(即在時刻ti活著的人有一半可望生存到的時間)、⑾第⑽列數(shù)據(jù)的標(biāo)準(zhǔn)誤差、⑿區(qū)間中點概率密度函數(shù)的估計值、⒀第⑿列數(shù)據(jù)的標(biāo)準(zhǔn)誤差、⒁區(qū)間左端點處危險概率的估計值、⒂第⒁列數(shù)據(jù)的標(biāo)準(zhǔn)誤差。

此處是用LT法對第２組資料進行統(tǒng)計描述的結(jié)果（從略），解釋方法同上。用３種檢驗方法給出的結(jié)果與用PL方法算得的結(jié)果相同，從略。

PL法與LT法的區(qū)別：在用計算機處理數(shù)據(jù)時，計算麻煩的苦惱已不存在，故PL法可適用于各種情況；當(dāng)用手工計算且樣本含量較大時，用LT法更方便一些。哩的計算結(jié)果基本上是一致的。PL法可看成是LT法的特殊情況，每個生存時間的區(qū)間寬度都為１。

[專業(yè)結(jié)論]因前面的圖示結(jié)果已表明此資料不服從指數(shù)，近似服從圖爾，故宜選用Log-Rank法或Wilcoxon法檢驗的結(jié)果。均拒絕H0、接受H1，即２組生存曲線之間差別顯著，無WBC倍增患者的生存期顯著地長于有WBC倍增的患者。

【例5.2.2】1965年某市腫瘤醫(yī)院總結(jié)隨訪了15年曾在該醫(yī)院住院手術(shù)的乳腺癌患者607例，整理后的資料如下，試分析該醫(yī)院乳腺癌患者手術(shù)后的生存率。

術(shù)后年數(shù):0～1～2～3～4～5～6～7～8～9～10～

期內(nèi)死亡人數(shù):59694330137144300

期內(nèi)失訪人數(shù):6371553831262111151222

【分析與解答】此生存資料以分組的形式給出，SAS程序可按如下方式編寫。

【SAS程序】━━【D5P2.PRG】

DATAabc;6359

KEEPftc;7169

RETAINt-0.5;5543

INPUTwithdrawfail;3830

t=t+1;3113

c=0;267

f=fail;2114

OUTPUT;114

c=1;153

f=withdraw;120

OUTPUT;220

CARDS;;

(程序的第１部分)(程序的第２部分)

PROCLIFETESTPLOTS=(S,LS,LLS,H,P)

INTERVALS=(0TO10)METHOD=LT;

TIMEt*c(1);

FREQf;

RUN;

(程序的第３部分)

因篇幅所限，程序修改指導(dǎo)、輸出結(jié)果及其解釋等項內(nèi)容從略。

第３節(jié)生存資料非參數(shù)統(tǒng)計方法中的有關(guān)計算公式

１．乘積─極限法(Product-Limit，PL法或稱為Kaplan-Meier，KM法)

讓t1<t2<…<tk代表離散的失效(死亡或復(fù)發(fā)等)時間，設(shè)ni為第ｉ個時刻開始之前生存的個體數(shù)目,即危險集的大小(i=1,2,…,k),再設(shè)di是在時刻ti失效的個體數(shù)目、si=ni-di。則在時刻ti的生存函數(shù)的PL估計值是ti時刻之前各時間點上生存率的乘積,即(5.2.1)

式中的估計量屬于右連續(xù)的，即在ti時刻發(fā)生的失效事件已包括在S^(ti)的估計中。與它對

應(yīng)的標(biāo)準(zhǔn)誤的估計值可用Greenwood的公式來計算：(5.2.2)

生存時間的第１樣本四分位數(shù)定義如下：(5.2.3)

第２、３樣本分位數(shù)可用類似的方式計算。q.50(即第２樣本分位數(shù))就是中位數(shù)，也就是中數(shù)生存期。

平均生存時間的估計值為∶(5.2.4)

式中t0=0。若最后一個是刪失數(shù)據(jù)，此式就低估了平均數(shù)。μ＾的標(biāo)準(zhǔn)誤差被定義為：(5.2.5)

式中。

２．壽命表法（Life-Table,LT法）

壽命表估計量通過計數(shù)落入時間區(qū)間[ti-1，ti]內(nèi)的失效和刪失的觀察例數(shù)來計算，這

里i=1,2,…,k+1,t0=0、tk+1=∞。令ni為進入?yún)^(qū)間[ti-1，ti]內(nèi)的個體數(shù)目、di和wi分別為發(fā)生在此區(qū)間內(nèi)的事件(指死亡或失效,下同)數(shù)目和刪失數(shù)目、bi=ti-ti-1、n'i=ni-wi/2；

n'i被稱為此區(qū)間內(nèi)有效的樣本大小。再令tmi為此區(qū)間的中點、p^i=1-q^i。事件在此區(qū)間內(nèi)發(fā)生的條件概率及其標(biāo)準(zhǔn)誤差的估計值分別由式(5.2.6)和(5.2.7)定義：(5.2.6)(5.2.7)

在時刻ti生存函數(shù)的估計值及其標(biāo)準(zhǔn)誤差的估計值分別由式(5.2.8)和(5.2.9)定義：(5.2.8)(5.2.9)

在tmi處的密度函數(shù)及其標(biāo)準(zhǔn)誤差的估計值分別由式(5.2.10)和(5.2.11)定義：(5.2.10)

(5.2.11)

在tmi處的危險率函數(shù)及其標(biāo)準(zhǔn)誤差的估計值分別由式(5.2.12)和(5.2.13)定義：(5.2.12)

(5.2.13)

設(shè)在區(qū)間[tj-1，tj]內(nèi)滿足關(guān)系式S^(tj-1)≥S^(ti)/2>S^(tj)，則在ti處的中數(shù)剩余壽命(即在ti活著的人有一半可望生存到的時間)及其標(biāo)準(zhǔn)誤差的估計值由式(5.2.14)和(5.2.15)

定義：(5.2.14)(5.2.15)

３．各層間生存曲線的齊性檢驗

(１)對數(shù)秩檢驗和Wilcoxon檢驗

檢驗各層之間齊性所用的秩統(tǒng)計量為：

v'V-v～χ2，df=R(V)(5.2.16)

式中R(V)（即協(xié)方差矩陣V的秩），從而可獲得近似的概率水平。這里，v是一個c×1的向量(v1,v2,…,vc)'，其具體表達式為：(5.2.17)]

這里c是分層變量的層數(shù)；協(xié)方差矩陣V=（Vjl）由下式來估計：(5.2.18)

這里下標(biāo)i表示離散的失效時間，當(dāng)j=l時，δjl=1;其他情況下，δ=0。nij和dij分別是第i個失效時間第j層中危險集的大小與事件數(shù)目，、、。Vj可以解釋為：在生存曲線相同的假設(shè)之下，觀察的與期望的失效數(shù)目之差的加權(quán)和。

當(dāng)權(quán)wi=1時為對數(shù)秩檢驗；當(dāng)權(quán)wi=ni時為Wilcoxon檢驗。V-代表V的廣義逆矩陣。由此可知：數(shù)值較大的失效時間在對數(shù)秩檢驗統(tǒng)計量中所起的作用大；而數(shù)值較小的失效時間在威爾科克森檢驗統(tǒng)計量中所起的作用大（因通常壽命短的頻數(shù)較大）。

(２)似然比檢驗

當(dāng)資料各層之間服從指數(shù)的假設(shè)成立時，檢驗各層之間齊性（即檢驗指數(shù)的尺度參數(shù)相等）的似然比檢驗統(tǒng)計量由式（5.2.19）定義：(5.2.19)

這里Nj是第j層內(nèi)事件的總例數(shù)，，是第ｊ層中用于檢驗的總時間、mj是第j層中觀察的總例數(shù)、。把Z視為服從自由度為c-1的卡方分布，從而求得近似的概率水平。

４．協(xié)變量與生存時間聯(lián)系密切程度的檢驗

用于檢驗協(xié)變量與生存時間聯(lián)系密切程度的秩檢驗是用于齊性檢驗的秩檢驗的更一般推廣，這種秩檢驗統(tǒng)計量具有如下的形式：(5.2.20)

式中v由式(5.2.21)定義、V-為V的廣義逆矩陣、V分別由式(5.2.22)與(5.2.24)定義：(5.2.21)

當(dāng)此式中時，為對數(shù)秩得分檢驗；

當(dāng)此式中時，為Wilcoxon得分檢驗。

式(5.2.21)中下標(biāo)及符號的含義如下：

α是觀測對象的編號，α=1,2,…,n；n是觀察的總數(shù)；i，j=1,2,…,k(k為不同的時間點數(shù))；t(j)代表有序的事件時間；Z(j)代表相應(yīng)的協(xié)變量向量；tα代表有序時間(含刪失與事件時間)。δα=1(如果觀察到事件發(fā)生)、δα=0（如果觀察到刪失發(fā)生）；得分Cα,δα取決于刪失的類型、且對全部觀察值求和。

用于對數(shù)秩統(tǒng)計量的協(xié)方差矩陣的估計量是：(5.2.22)

式中Vi是與t(j)時危險集所對應(yīng)的校正的平和與交叉乘積和構(gòu)成的矩陣，即(5.2.23)

式中

用于Wilcoxon統(tǒng)計量的協(xié)方差矩陣的估計量為下式(其編號為∶(5.2.24))∶此式中ai、a*i、Si、xi、si分別為：、、

、

第３章ＣＯＸ模型回歸分析

第１節(jié)ＣＯＸ回歸模型(半?yún)?shù)回歸模型)

像通常的回歸分析一樣，人們也希望能建立起生存時間（因變量或反應(yīng)變量）隨危險因素（自變量或協(xié)變量）變化的回歸方程，以便對危險因素的作用大小有一個全面的了解和掌握、并根據(jù)危險因素的不同取值對生存概率(或危險率)進行預(yù)測。由于生存時間的準(zhǔn)確很難獲得，前述目的很難直接實現(xiàn)。1972年ＣＯＸ提出了比例危險模型（ProportionalHazardModel），簡稱為COX模型。由于此模型在表達形式上與參數(shù)模型相似，但在對模型中各參數(shù)進行估計時踴依賴于特定的假設(shè)，所以又有半?yún)?shù)模型之稱。此模型的實用面很寬，在生存分析中占有特殊的地位。其模型的具體形式為：

hi(t)=h0(t)exp(β1xi1+β2xi2+…+βmxim)(5.3.1)

式中hi(t)為第i名受試者生存到ti時刻的危險率函數(shù)，h0(t)是當(dāng)所有危險因素(即xij=0)不存在時的基礎(chǔ)危險率函數(shù)，X=(xi1,xi2,…,xim)'是可能與生存時間有關(guān)的ｍ個危險因素所構(gòu)成的向量。將式(5.3.1)變形如下：

ln[hi(t)/h0(t)]=β1xi1+β2xi2+…+βmxim(5.3.2)

此式表明：各危險因素與回歸系數(shù)的線性組合就是第ｉ名受試者的相對危險率函數(shù)的自然對數(shù)值。再設(shè)有i、j２個受試者,其危險因素向量分別為X1與X2，由式(5.3.1)不難得出他們的相對危險率的自然對數(shù)為：

ln[hi(t)/hj(t)]=β1(xi1-xj1)+β2(xi2-xj2)+…+βm(xim-xjm)(5.3.3)

即利用“具有某預(yù)后因素向量的受試者的死亡風(fēng)險與不具有該預(yù)后因素向量的受試者的死亡風(fēng)險在所有時間上都保持一個恒定比例”的假設(shè)，巧妙地獲得了各時間點上２個受試者相對危險率函數(shù)的估計值。

然而，當(dāng)資料不滿足上述假設(shè)時，即有些危險因素作用的強度是隨時間而變化的，２個受試者的危險率函數(shù)之比(相對危險)隨時間而改變，就應(yīng)改用時變協(xié)變量模型，也稱為非比例危險模型(Nonproportionalhazardmodel)。當(dāng)只有一個危險因素時，其模型的具體形式為：

hi(t)=h0(t)exp[βxi+γ(xiti)](5.3.4)

式中ti為第i個受試者的生存時間。

上述各式中的回歸系數(shù)需用最大似然法進行估計，一旦有了危險率函數(shù)的估計值，再利用生存時間函數(shù)之間的相互關(guān)系，可獲得其他生存時間函數(shù)的估計值。

第２節(jié)ＣＯＸ模型回歸分析應(yīng)用舉例

[例5.3.1]某醫(yī)院腫瘤科提供的一份關(guān)于肺癌病人的失效時間資料，因變量（或反應(yīng)變量）為病人治療后的生存時間ｔ(天)，當(dāng)ｔ為刪失數(shù)據(jù)時，用前面加一個負(fù)號來表示；考察的協(xié)變量（即危險或預(yù)后因素）如下。

①癌細胞的類型(Cell$)，它有４個水平，即腺癌細胞(adeno)、鱗癌細胞(squamous)、小細胞肺癌(small)和大細胞肺癌(large)；

②治療類型(THERAPY$)，它有２個水平，即標(biāo)準(zhǔn)的方法(standard)和試驗的方法(test)；

③療前處理(PRIOR$)，它有２個水平，即采取了療前處理(yes)和未采取療前處理(no)；

④病人的年齡(age)(歲)；

⑤從診斷到治療的等待時間(diagtime)；

ⅰ人的行動狀態(tài)用Karnofsky率來度量,其取值用KPS表示，10≤KPS≤30表明病人完全靠醫(yī)院護理、40≤KPS≤60表明病人的行動部分地受到限制、70≤KPS≤90表明病人的行動可以自理。

前３個變量被當(dāng)作分類變量，后３個變量被當(dāng)作連續(xù)性變量。資料的形式為：

各組病人的治療方法、癌細胞類型,同一組中的樣本含量

生存時間、KPS值、diagtime值、年齡、與療前處理對應(yīng)的指示變量PR值

（注：PR=0等價于令PRIOR='YES'，即表示采取了療前處理；PR=10等價于令PRIOR='NO'，即表示未采取療前處理）。全部數(shù)據(jù)如下（文件名：lung.dat）：

STANDARDSQUAMOUS15

7260769041170564102286033801266096310

118701165101020549082401069101108029680

3145018430-10070670042604810840586310

144304630-2580952101170114810

STANDARDSMALL30

306036103846094204402350548046310

13604560-123403550-9760567015360146310

593026501178034601630453101515012690

226046805680124310214025510182015420

139802640203056503175365052702550

2876025661018304600516016701228028530

27608620547016707507720635011480

3924046801040236710

STANDARDADENO9

82019611092701060035406620117802380

1328055001250463101628056403303430

95804340

STANDARDLARGE15

177501666101628056202165015520553702470

2786012630124012681026080545020080124110

156702660-1829026201439086001058011660

103805380250708531010060133710

TESTSQUAMOUS20

99990125410112806600-87803480-2315085210

2425017009917075010111703620120216510

58760358038990262033306640252036630

3577013580467902640201802852101507350

30701163044601370102839025101550134010

TESTSMALL18

25302690-103702236102120471013302620

8760260024036441020309541072011660

2460849099703720880268099854620

617027102570270095701610805017710

513087591029408670

TESTADENO18

24402600184056910-8399357031803390

5160562090602250105260343073603700

850566036708610481048107404580

140703630186903600848046210195010420

4540369080404630

TESTLARGE12

5260445016470156810193043910536012660

15305630436011491034080106410133751650

1116056402317018671037880465049303370

[SAS程序]──[D5P3.PRG]

DATAvalung;PROCPHREG;

RETAINtherapycell;MODELt*censor(1)=kpsagediagtime;

LENGTHprior$3;RUN;

INFILE'a:lung.dat';PROCPHREG;

INPUTtherapy$cell$n;MODELt*censor(1)=kpsagediagtime;

DOi=1TOn;STRATAcell;

INPUTtkpsdiagtimeagepr@@;RUN;

censor=(t<0);PROCPHREG;

t=ABS(t);MODELt*censor(1)=kpsagediagtime;

IFpr=10THENprior='yes';STRATAcelltherapyprior;

ELSEprior='no';RUN;

OUTPUT;

END;

（程序的第１部分）（程序的第２部分）

[程序修改指導(dǎo)]調(diào)用PHREG過程時，MODEL語句等號右邊必須是連續(xù)性變量。這里所寫的３個過程步的區(qū)別在于STRATA語句中所包含的楓變量的個數(shù)：第１個過程步不含STRATA語句，就是把所有資料看成來自１個總體；第２個過程步要求按CELL的４個水沏分析資料；第３個過程步要求按CELL、THERAPY、PRIOR的16種水平組合楓分析資料。顯然，在分層的條件下，COX所作的比例危險假設(shè)容易得到滿足，但各層的樣本含量不應(yīng)太小。

另外，若連續(xù)變量很多時，可在MODEL語句最后增加選擇項/SELECTION=方法名，進行變量篩選。方法名有如下幾種：

BACKWARD或B(后退法)、FORWARD或F(前進法)、STEPWISE或S(逐步回歸法)、SCORE(最優(yōu)回歸子集法)。

[輸出結(jié)果及其解釋]ThePHREGProcedure

TestingGlobalNullHypothesis:BETA=0

WithoutWith

CriterionCovariatesCovariatesModelChi-Square

-2LOGL397.545363.88133.665with3DF(p=0.0001)

Score..34.158with3DF(p=0.0001)

Wald..31.510with3DF(p=0.0001)

此結(jié)果表明：含３個自變量的COX回歸模型是顯著的（P<0.0001）。

AnalysisofMaximumLikelihoodEstimates

ParameterStandardWaldPr>Risk

VariableDFEstimateErrorChi-SquareChi-SquareRatio

KPS1-0.0361270.0065530.446860.00010.965

AGE1-0.0215710.011573.476020.06230.979

DIAGTIME10.0064260.012450.266240.60591.006

這是用最大似然法對模型中各參數(shù)估計并檢驗的結(jié)果，顯然，只有變量KPS是顯著的（P<0.0001）。這里給出的只是最后１個過程步輸出的結(jié)果，因前２個過程步輸出的結(jié)果與此結(jié)果相似，從略。詳細的討論請看下面的例子。

[例5.3.2]若用LIFETEST過程預(yù)處理此資料會發(fā)現(xiàn):不同細胞類型的生存曲線之間差別非常顯著，生存曲線從左到右依次為：ADENO、SMALL、LARGE、SQUAMOUS。在[例5.3.1]中，請設(shè)法將變量CELL引入COX模型中來，重建模型。

[分析與解答]因分類變量無法直接放入回歸方程，這對模型的擬合是不利的。補救的辦法是：引入啞變量，使分類變量轉(zhuǎn)變成數(shù)值袖量后再用PHREG過程（請看本例）;若已知生存時間近似服從某特定的參數(shù)模型時，可直接用LIFEREG過程擬合參數(shù)模型,因參數(shù)模型中可包含２類變量（請看[例5.4.1]）。于是,可對變量CELL作變換，使它變成３個啞變量。對數(shù)據(jù)文件作如下修改（改程序較困難）:把ADENO改成100；把SMALL改成010；把LARGE改成001；把SQUAMOUS改成000。如第１行

STANDARDSQUAMOUS15應(yīng)改為：STANDARD00015，其他如法炮制。

程序[COXLUNG.PRG]的數(shù)據(jù)步中第１個INPUT語句和過程步需作如下修改

INPUTtherapy$adenosmalllargen;

程序[COXLUNG.PRG]中的過程步需作如下修改(因其他連續(xù)變量的作用不顯著，故未將它們寫入下面的模型語句之中)：

PROCPHREG;

MODELt*censor(1)=kpsadenosmalllarge;

STRATAtherapyprior;

RUN;

[輸出結(jié)果及其解釋]