常用生存分析模型及其對時依性協(xié)變量效應的估計方法*

肖媛媛 許傳志 趙耐青

常用生存分析模型及其對時依性協(xié)變量效應的估計方法*

肖媛媛1許傳志1趙耐青2△

生存分析是利用統(tǒng)計學相關理論和方法探索研究因素與“事件時間”（time-to-event）關聯(lián)的問題。自20世紀50年代其研究方法初具規(guī)模以來，經過幾十年、尤其是近三十年來的蓬勃發(fā)展，生存分析已經成為現(xiàn)代統(tǒng)計學的一個重要分支，研究領域廣泛涉及醫(yī)學、生物學、工程學、經濟學及人口統(tǒng)計學等［1］。按照理論基礎的不同，和很多統(tǒng)計推斷方法一樣，自誕生之初起，生存分析也沿著頻率法（frequentistmethod）和貝葉斯法（Bayesian method）兩條道路分別演進，但發(fā)展速度和軌跡存在較大差別。

在貝葉斯生存分析方法中，由似然函數(shù)（likelihood function）和參數(shù)的先驗分布（prior distribution）所共同構建的后驗分布（posterior distribution）的表達式往往比較復雜，涉及高維重積分的運算，因此在計算機技術欠發(fā)達的年代，這一類生存分析方法的發(fā)展幾近停滯。20世紀80年代，有學者陸續(xù)提出了一些簡化后驗分布密度及各階矩計算的近似算法，如Lindley數(shù)值逼近法、Naylor-Smith逼近法、Tierney-Kadane逼近法等［2－4］，但這些算法的實現(xiàn)仍涉及十分復雜的近似求解技術。直到2000年以后，隨著計算機技術的發(fā)展和貝葉斯方法的改進，特別是馬爾科夫鏈蒙特卡洛法（Markov Chain Monte Carlo，MCMC）的成功運用后［5－6］，貝葉斯生存分析法才真正開始快速發(fā)展。盡管如此，與其他任何貝葉斯統(tǒng)計推斷方法一樣，如何準確定義先驗分布的問題仍然是限制其進一步發(fā)展的最大阻礙。由于貝葉斯生存分析法在實際運用中需要較為精深的數(shù)理統(tǒng)計理論知識作為支撐，其推廣運用必然受限?；谏鲜霰尘?，本文將不再討論這一類生存分析方法，而將重點放在常用的頻率生存分析法上。

與貝葉斯生存分析法不同，頻率生存分析法嚴格忠于數(shù)據本身，因此不會遭受“客觀性”的質疑。此外，其理論基礎沒有貝葉斯法抽象晦澀，計算方法也較貝葉斯法簡便得多。因此，不論是在當前的生存分析模型的方法學研究領域，還是在具體生存模型的運用領域，頻率生存分析法都占絕對主導地位。按照算法或模型的架構是否依賴特定參數(shù)分布，頻率生存分析法分為非參數(shù)法（non-parametric method）、半參數(shù)法（semi-parametric method）和參數(shù)法（parametric method）三類。

在研究某一因素與生存時間的關聯(lián)時，這一因素的取值或分類可以是恒定不變的，也可以是呈一定規(guī)律或無規(guī)律變化的。另外，該因素對結局變量的效應也可能隨時間發(fā)生變化，這一類因素可被統(tǒng)稱為時依性協(xié)變量（time-dependent covariates）。以醫(yī)學研究為例，性別、種族、成年人身高等因素可以認為在一定研究時期內保持不變，而各類血液化驗指標、體重、血壓等體格測量指標則很可能是不斷變化的。要評價這些變動的因素對生存時間的影響，必須將其變動的信息充分利用起來，才能得到更為準確的結論。絕大多數(shù)經典生存分析模型在提出之初都可以看作是“靜態(tài)生存分析模型”，因為為了簡化理論推導，一般都會給出協(xié)變量取值及對結局變量效應恒定的假設。在其后的運用中，再根據協(xié)變量的變動規(guī)律或近似變動規(guī)律對原始“靜態(tài)生存分析模型”進行拓展，形成各類“動態(tài)生存分析模型”。以Cox比例風險模型為例，在原始模型提出之后，Kalbfleisch和 Prentice［7］，以及 Cox本人和Oakes［8］就對該模型進行了改良，特別是討論了時依性協(xié)變量的分析方法，大大提升了原始模型的實際應用價值。時至今日，時依性協(xié)變量的分析方法在多類生存分析模型中均有涉及，本文將對常用生存分析方法及其對時依性協(xié)變量的效應估計做簡要梳理。

非參數(shù)法

1.常用方法

非參數(shù)生存分析方法在理論架構上有共性，即:不使用原始數(shù)據中生存時間分布及具體時間長短等信息，而僅僅利用不同個體生存時間的秩次排序。最常見的非參數(shù)生存分析法有用以描述生存時間分布的乘積極限法（product limit method），也稱作 Kaplain-Meier法［9］，以及比較兩組或多組生存時間差異的一組非參數(shù)檢驗，如對數(shù)秩檢驗（log-rank test）、Wilcoxon檢驗及 Mantel-Haenszel檢驗［7，10－11］。

（1）Kaplain-Meier法

假設在一組觀察對象中，到觀察結束時，一共有m個個體在j個時間點發(fā)生了死亡（或其他任何所研究的結局事件，下文均以死亡為例），且死亡時間點的排序為:0≤t（1）≤t（2）…≤t（j）≤∞。假設恰好在某個特定死亡時間點t（j）之前，仍有 r（j）個觀察對象有死亡風險（意味著仍然存活且未發(fā)生截尾），而在t（j）時間發(fā)生了d（j）例死亡，則據此可估算t時刻的生存函數(shù)值為:

這一取值也被稱作K-M估計值（K-M estimator）。不難看出，K-M估計值從定義來看在時間維度上應該是一個離散函數(shù)。如假設所有死亡事件均恰好發(fā)生在死亡時點，而兩個離散的死亡時點間無死亡發(fā)生，則可根據K-M估計值將生存曲線繪制為連續(xù)的、逐步遞減的階梯函數(shù)（step function），只在發(fā)生死亡的時點變化數(shù)值。對于離散時點來說，可以證得K-M估計值實際上也是最大似然估計值（maximum likelihood estimate）［8］。

K-M估計值可被證明服從近似正態(tài)分布［12－13］，因此可以計算其可信區(qū)間。Greenwood提出了計算其可信區(qū)間的近似公式（Greenwood′s formula）［14］:

這一公式也可用來計算百分位數(shù)生存時間的可信區(qū)間。

（2）以秩次為基礎的非參數(shù)檢驗

對數(shù)秩檢驗、Wilcoxon檢驗及Mantel-Haenszel均是以秩次為基礎的非參數(shù)檢驗?？紤]其原理的大同小異，只簡略介紹最為常見的對數(shù)秩檢驗。

對數(shù)秩檢驗的基本思想為:假設有兩組觀察對象，將觀察時間內兩組發(fā)生的所有死亡個體提出，按照死亡時間從短到長混合排序0≤t（1）≤t（2）…≤t（j）≤∞。假設t（j）之前，兩組加在一起一共有 r（j）個對象有死亡風險，其中第一組有 r（1j）個，第二組有 r（2j）個。t（j）這一時刻一共發(fā)生了 d（j）例死亡，其中第一組 d（1j）例，第二組 d（2j）例。在無效假設成立的背景下，任意時刻發(fā)生死亡的風險在兩組間應相同（只存在隨機抽樣誤差），因此可以用在任意時刻的死亡風險 d（j）/r（j）來估算第一組的期望死亡數(shù)（e1j），并用每一時刻期望死亡數(shù)和觀測死亡數(shù)（d1j）的差值構建檢驗統(tǒng)計量UL:

在每一死亡時點上，每一組死亡數(shù)均服從參數(shù)為r（j）和 d（j）的超幾何分布（hypergeometric distribution），據此可計算得到UL的方差為:

在無效假設成立的條件下，UL服從正態(tài)近似分布N（0，VL），即可參照正態(tài)分布做出統(tǒng)計推斷。

Mantel-Haenszel檢驗與對數(shù)秩檢驗的主要差別只體現(xiàn)在，其檢驗統(tǒng)計量值和方差均為對數(shù)秩檢驗統(tǒng)計量和方差的平方，因此依據卡方分布做出統(tǒng)計推斷。而Wilcoxon檢驗與對數(shù)秩檢驗的區(qū)別是，其計算檢驗統(tǒng)計量及其方差時乘以每一時刻存活病人總數(shù)（r（j））作為權重。近年來，國內有學者討論了無刪失存在時Wilcoxon檢驗與對數(shù)秩檢驗的I型錯誤率和檢驗效能，模擬結果表明兩法存在一定差別［15］。

2.對時依性協(xié)變量效應的估計

K-M法一般用來繪制生存時間分布較多，而基于秩次的幾種非參數(shù)檢驗常用來初步比較單個分類變量或少數(shù)幾個分類變量的聯(lián)合分布對生存時間的影響。一旦分類變量數(shù)量較多且每個變量分類較多，由于分層后數(shù)據過于稀疏，其檢驗效能往往會大打折扣。此外，更是無法分析連續(xù)性變量對生存時間的影響。在如此多的固有缺陷下，想要處理時依性協(xié)變量就顯得更不實際了。在未限定發(fā)表時間，以三個關鍵詞“nonparametric”、“survival analysis”、“time dependent”組合檢索Web of Science后，經過標題和摘要篩選，我們只尋找到了一篇切題文獻，但經過全文閱讀后，發(fā)現(xiàn)作者采用的方法并非嚴格意義的非參數(shù)法，而是半參數(shù)法［16］。

半參數(shù)法

1.常用方法

（1）半參數(shù)比例風險模型（Cox proportional hazards model）

比例風險模型的構建滿足如下假設:在基礎風險和其他協(xié)變量固定不變的前提下，某一協(xié)變量取值每增加一個單位，得到的風險函數(shù)的取值等于原來的風險函數(shù)取值乘以一個固定系數(shù)。其表達式十分簡潔:

上式中，h0（t）是基礎風險函數(shù)。在半參數(shù)比例風險模型中，對基礎風險函數(shù)（baseline hazard function）不做任何限定，只是對各研究因素（x1～xp）對基礎風險的作用做出參數(shù)限定（βT）。這也是“半參數(shù)模型”這一名稱的由來，這一模型最初是由英國著名統(tǒng)計學家Cox提出的，因此也稱作Cox比例風險模型［17］。

由于未限定基礎風險函數(shù)，也就是未限定生存時間分布，不可能得到概率密度值（考慮生存函數(shù)、風險函數(shù)、概率密度函數(shù)三者間的轉化關系，后文將會提及），因此在估計參數(shù)時，傳統(tǒng)的極大似然法無法使用。Cox采用了一個非常巧妙的處理來化解這個問題，假設每個死亡時點只有一個死亡對象，則其概率密度估算值可用風險函數(shù)表示為:

P（i對象在 t（j）時刻死亡｜在 t（j）時刻存在風險的R（j）個對象至少有一個死亡）

由于上式分子分母均同樣含有t（j）時刻的基礎風險，因此可以約除。假設一共有n個死亡事件，則似然函數(shù)可表示為:

上式的實質是通過風險函數(shù)所構建的子集似然函數(shù)（profile likelihood function），并不是真正意義上的似然函數(shù)。因為不包含所有的未知參數(shù)，故也被稱作偏似然函數(shù)（partial likelihood function）［18］。在偏似然函數(shù)中，未作限定的基礎風險值已經不存在，不需要對其進行估計。與普通似然函數(shù)求極值相似，仍然可通過迭代法（如New ton-Raphson法）尋找其最大值，從而得到β′的極大偏似然估計。

（2）半參數(shù)加速失效模型（sem i-parametric accelerated failure timemodel）

加速失效模型是另一類常見的生存分析模型，但其在實際研究中的應用卻遠沒有比例風險模型廣泛。加速失效模型的目標函數(shù)直接選擇為生存時間T，因此較之比例風險模型，其協(xié)變量系數(shù)的解釋更為明確。加速失效模型構建的假設為:在基線生存函數(shù)及其他協(xié)變量恒定的情況下，某一協(xié)變量每增加一個單位，生存時間等于原來的生存時間乘以一個固定系數(shù)。

模型的表達式一樣不失簡潔:

上式中，ε是隨機誤差，對應基線生存函數(shù)的對數(shù)。xl～xp是研究對象的協(xié)變量取值，a1～ap是協(xié)變量的系數(shù)。

在半參數(shù)加速失效模型中，對基線生存函數(shù)的分布未做任何限定，只限定協(xié)變量的系數(shù)。在滿足加速失效模型假設的前提下，風險函數(shù)的表達式為:

前面提到，在Cox比例風險模型中，雖然一樣未限定基線風險函數(shù)，但比例風險的假設可以在用風險函數(shù)的比值構建偏似然函數(shù)時將基線風險順利約除。而在上式中，誤差項exp（ε）的風險函數(shù)值λ（·）完全未知。因此無法在半參數(shù)加速失效模型中用風險函數(shù)如法炮制偏似然函數(shù)。正是由于這一限制，參數(shù)加速失效模型反而比半參數(shù)模型應用更為廣泛。不過在過去三十年里，許多統(tǒng)計學家提出了一系列以秩次和最小二乘為基礎的半參數(shù)加速失效模型參數(shù)估計和推斷方法，其中比較有代表性的有 Prentice，Buckley和James，Ritov，Tsiatis，Lai和 Ying等［19－24］。除去計算極其復雜之外，這些方法存在一個共同問題，那就是均為離散函數(shù)，可能存在多重根，難以對函數(shù)值及其方差做出估計。

（3）其他半參數(shù)生存分析模型

其他較為常見的半參數(shù)生存分析模型還包括比例優(yōu)勢模型（proportional odds model）和相加風險模型（additive hazardsmodel）。其模型表達都比較簡潔，但參數(shù)估計時所涉及的運算卻極具挑戰(zhàn)性，在此不做詳述。如比例優(yōu)勢模型參數(shù)估計的常用方法，是在給出一系列分布限制條件的基礎上（從這個意義上說，半參的假定已經被部分違背）構建似然函數(shù)，并且該似然函數(shù)也只是在進一步的限制條件下才存在最大值［25］。再如相加風險模型，雖然可以順利構建似然函數(shù)，但似然函數(shù)的最大值卻不存在，而隨后提出的一些近似估計方法的表現(xiàn)也差強人意［26－27］。

2.對時依性協(xié)變量效應的估計

（1）Cox比例風險模型

在Cox比例風險模型中，目前有兩大類處理時依性協(xié)變量的方法。第一類不妨稱其為“函數(shù)法”。其思路很直接:如果一個協(xié)變量在生存時間的維度隨時間的變化而變化，則可以在原始Cox比例風險模型中加入該變量和時間的交互作用項來描述其變化對基線風險的影響。調整后的風險函數(shù)表達式為:

上式中，xj為時依性協(xié)變量，其在t時刻的取值為xj（t＝0）×g（t）。從該表達式不難看出，在控制了 xj的變化對基礎風險的影響之后，其他協(xié)變量仍滿足比例風險假定。

這種方法的運用需要滿足兩個前提:一是可以準確找到協(xié)變量xj隨時間變化的函數(shù)關系式g（t）。在實際研究中，一般都是通過獲取的數(shù)據來尋找變量間的函數(shù)關系。當兩者關系為一次或者二次函數(shù)時，找準關系或許不難。而一旦多次函數(shù)出現(xiàn)，在數(shù)據本身就夾雜了抽樣誤差和各類偏倚的情況下，準確定位其函數(shù)關系顯然無法實現(xiàn)。二是協(xié)變量xj每改變一個單位，對基礎風險的影響保持恒定不變。這同樣是一個極強的假設，以醫(yī)學研究為例，往往某一指標發(fā)生變化后，其作用機理也會發(fā)生改變，由此導致效應強度發(fā)生變化。

另一種處理方法是借鑒counting process法的思路調整生存數(shù)據后，再運用Cox模型進行分析，并在模型參數(shù)估計的時候運用多結局生存分析的思路校正觀測間的組內相關性。生存分析中有著廣泛應用的counting process法是由 Fleming和 Harrington［28］于 20世紀90年代初在Aalen［29］的研究基礎上發(fā)展成熟的。這一方法給生存分析理論的拓展帶來了極大推進，如基于counting process法計算得到的鞅殘差（martingale residual）估計值可以作為模型擬合的評判標準，再如運用鞅論的相關公理，可證得協(xié)變量系數(shù)的極大似然估計值服從近似正態(tài)分布，且其協(xié)方差矩陣可以從觀測到的信息矩陣（information matrix）中估計獲得［30］。這里提到的方法，僅僅是借鑒了counting process的思路來重構生存數(shù)據，并不實際運用原方法，為了加以區(qū)別，不妨稱其為“多結局counting process法”。

多結局counting process法重構生存數(shù)據的基本思想是:假設變動的因素為x，研究對象A結局事件發(fā)生時間為t（d），在結局事件發(fā)生前，一共隨訪了3次（t（1）＜t（2）＜t（3）＜t（d）），每次測得 x的取值分別為 x（1），x（2）和 x（3）。假設 x（1）＝x（2）≠x（3）。由此可知研究因素x取值在t（3）發(fā)生了變化，則以t（3）為分界點，將原來數(shù)據集中研究對象A所對應的一條記錄（t（d），1）（數(shù)據結構為（T，C），其中 T代表隨訪時間，C是指示變量，1代表死亡，0代表刪失），裂解為新生存數(shù)據集中的兩條記錄（t（1），t（3），0）和（t（3），t（d），1）（數(shù)據結構為（T（1），T（2），C），其中 T（1）是隨訪開始時點，T（2）為隨訪結束時點，C是指示變量，1代表死亡，0代表刪失）。如此一來，每個研究對象從原始生存數(shù)據集中的一條記錄變?yōu)樾律鏀?shù)據集中的一組記錄，可以視為同一研究對象的多結局事件。

在重構的新生存數(shù)據集中，每一條記錄所對應的x為恒定取值，因此可根據不同假設采用Cox比例風險模型進行參數(shù)估計。常見的假設有兩種，一是假設各隨訪開始時點（T（1））的基線風險均相同，最后得到研究因素的唯一效應估計值；二是假設各隨訪開始時點的基線風險不相同，最后得到的是不同時刻研究因素效應的一組估計值。在進行參數(shù)可信區(qū)間估計時，考慮到新生存數(shù)據集中來自同一個觀察對象的多條記錄存在內部相關，采用“夾心方差估計”（sandw ich variance estimator）調整方差［31－34］。

雖然多結局counting process法在應用時也有較強的假設作為前提，最為突出的一點就是假設每個研究對象的各個“結局”間互相獨立。但較之前述及的“函數(shù)法”，顯然靈活和強大得多。在完成數(shù)據結構調整后，分析方法與傳統(tǒng)Cox比例風險模型相同，僅僅需要額外調整組內相關性。因此采用常用統(tǒng)計分析軟件（如SAS或R）的生存分析模塊即可輕松實現(xiàn)，具體操作方法和程序編碼目前也有相關文獻可供參考［35－36］。

（2）半參數(shù)加速失效模型

前面已經提到，在估計恒定不變的協(xié)變量的效應時，以秩次和最小二乘為基礎的參數(shù)估計方法所涉及的計算已經極其復雜且表現(xiàn)不佳，在此基礎上，不可能再將時依性協(xié)變量整合其中。直到2007年，Zeng和Lin提出了一個在半參數(shù)加速失效模型中估計時依性協(xié)變量效應的可行方法［37］。他們首先將時依性變量整合進半參數(shù)加速失效模型中，得到:

上式中，λ（·）是誤差項exp（ε）的未知基線風險函數(shù)，同樣，偏似然函數(shù)仍無法構建。他們提出了采用核平滑（kernel smoothing）來構建子集似然函數(shù)的平滑近似函數(shù)，再以此平滑近似函數(shù)為基礎進行參數(shù)估計。隨后的數(shù)據模擬發(fā)現(xiàn)，據此方法得到的近似估計值較為準確和穩(wěn)定。雖然這一方法的提出為半參數(shù)加速失效模型中時依性協(xié)變量的效應估計提供了可能，但如此復雜的運算顯然無法推廣應用。

參數(shù)法

1.常用方法

參數(shù)生存分析模型和半參數(shù)生存分析模型的最大區(qū)別在于，限定生存函數(shù)（survival function，S（t））滿足一定概率分布，由此，風險函數(shù)（hazard function，h（t））及結局事件概率密度函數(shù)（probability density function，f（t））也相應確定。因為三個函數(shù)滿足以下轉換關系:

目前常用的參數(shù)生存分析模型分為兩大類:比例風險模型和加速失效模型。

（1）參數(shù)比例風險模型（parametric proportional hazards model）

參數(shù)比例風險模型和半參數(shù)的Cox比例風險模型在表達式上基本相同:

唯一的區(qū)別為:參數(shù)比例風險模型假設基礎風險h0（t）滿足一定的概率分布。可以看出，該模型仍然假設協(xié)變量對基礎風險的改變等于固定的乘積比例。滿足“比例風險”這一理論假設的生存分布有指數(shù)分布（exponential distribution）、韋伯分布（Weibull distribution）和 Gompertz分布（Gompertz distribution），因此參數(shù)比例風險模型也分為這三種。

一旦生存分布選定，即可根據生存函數(shù)和概率密度函數(shù)之間的轉換關系，構建似然函數(shù)，得到模型參數(shù)的極大似然估計。

（2）參數(shù)加速失效模型（parametric accelerated failure time model）

參數(shù)加速失效模型與半參數(shù)加速失效模型的本質區(qū)別也在于限定了基線生存函數(shù)的分布，即函數(shù)表達式中εi項的分布。其表達式為:

滿足加速失效模型假設的生存函數(shù)分布類型較滿足比例風險假設的分布類型多，有指數(shù)分布、韋伯分布、對數(shù)-logistic分布（log-logistic distribution）、對數(shù)-正態(tài)分布（log-normal distribution）和伽馬分布（Gamma distribution）。

同樣，生存時間分布一旦限定，即可構建似然函數(shù)，從而得到模型參數(shù)的極大似然估計。

2.對時依性協(xié)變量效應的估計

在估計時依性協(xié)變量x的效應時，兩類參數(shù)生存分析模型往往在假設x的單位增量對應變量（風險或生存時間）的效應保持不變、且x遵循一定的隨時間變化規(guī)律（f（t））的情況下，在原生存分析模型中加入x與f（t）的交互作用項。再以調整后的函數(shù)為基礎構建似然函數(shù)，采用極大似然法得到協(xié)變量效應估計值。前面述及的“多結局counting process法”當然也可以順利應用于參數(shù)生存分析模型，但由于Cox模型在實際運用中的巨大優(yōu)越性，counting process法在參數(shù)比例風險模型中的探討較少。一些學者討論了其在參數(shù)加速失效模型中的運用［38－39］。

可能是由于假設生存分布已知，掃除了生存分析模型參數(shù)估計時的理論障礙，目前參數(shù)生存模型中針對時依性協(xié)變量效應估計的研究并不多。近年來有學者開始討論在弱假設前提下時依性協(xié)變量的效應估計，如Sparling和Hamid等提出了將樣條函數(shù)（Spline function）整合入參數(shù)生存分析模型，并結合不同數(shù)據刪失類型，估計非線性變化的時依性協(xié)變量的效應［40－41］。

總 結

雖然在個別生存分析模型，如半參數(shù)加速失效模型中，缺乏可行的時依性協(xié)變量效應的估計方法。目前對于大多數(shù)常見的半參數(shù)和參數(shù)生存分析模型來說，借助常用統(tǒng)計分析軟件中已有的生存分析模塊，可順利實現(xiàn)對時依性協(xié)變量的效應估計。現(xiàn)有研究數(shù)量的不足及理論瓶頸決定了時依性協(xié)變量效應估計新方法的探索，以及現(xiàn)有方法的簡化和拓展將是未來很長一段時期內生存分析方法學研究的重點和難點之一。此外，一些現(xiàn)有復雜算法在常用統(tǒng)計分析軟件中的實現(xiàn)也將是統(tǒng)計開發(fā)人員需直面的技術挑戰(zhàn)。

［1］林靜.基于MCMC的貝葉斯生存分析理論及其在可靠性評估中的應用.南京理工大學，2008.

［2］Lindley DV.“Approximate Bayesian methods”in Bayesian statistics.Valencia，Spain:Valencia Press，1980.

［3］Naylor JC，Sm ith AFM.Applications of a method for the efficient computation of posterior distributions.Applied Statistics，1982，31（2）:214-225.

［4］Tierney L，Kadane JB.Accurate approximations for posterior moments and marginals.Journal of American Statistical Association，1986，81（1）:82-86.

［5］Congdon P.Bayesian statistical modeling.England:John Wiley and Sons，2001.

［6］Congdon P.Applied Bayesian modeling.England:John Wiley and Sons，2003.

［7］Kalbfleisch JD，Prentice RL.The statistical analysis of failure time data.New York:John W iley＆Sons，1980.

［8］Cox DR，Oakes D.Analysis of survival data.London:Chapman and Hall，1984.

［9］Kaplan E，Meier P.Nonparametric estimation from incomplete observations.Journal of the American Statistical Association，1958，53（282）:457-481.

［10］Mantel N.Evaluation of survival data and two new rank order statistics arising in its consideration.Cancer Chemotherapy Report，1966，50（3）:163-170.

［11］Gehan EA.A generalized Wilcoxon test for comparing arbitrarily singly censored samples.Biometrika，1965，52:203-223.

［12］Andersen PK，Borgan O，Gaill RD，et al.Statistical methods based on counting processes.New York:Springer，1993.

［13］Flem ing TR，Harrington DP.Counting processes and survival analysis.New York:John W iley＆Sons，1991.

［14］Greenwood M.A report on the natural duration of cancer.Reports on Public Health and Medical Subjects 33.London:HM Stationary Office，1926:1-26.

［15］陳靖，何春拉，潘建紅，等.無刪失生存數(shù)據Wilcoxon秩和檢驗與logrank檢驗的比較.中國衛(wèi)生統(tǒng)計，2012，29（5）:654-660.

［16］Zucker DM，Karr AF.Nonparametric survival analysis with time-dependent covariate effects:a penalized partial likelihood approach.The Annals of Statistics，1990，18（1）:329-353.

［17］Cox DR.Regression models and life-tables.Journal of the Royal Statistical Society.Series B，1972，34（2）:187-220.

［18］Cox DR.Partial likelihood.Biometrika，1975，62（2）:269-276.

［19］徐英，駱福添.Buckley-James模型在生存分析中的應用.中國衛(wèi)生統(tǒng)計，2007，24（1）:69-70.

［20］Prentice RL.Linear rank tests with right-censored data.Biometrika，1978，65（1），167-179.

［21］Buckley J，James I.Linear regression with censored data.Biometrika，1979，66（3），429-436.

［22］Ritov Y.Estimation in a linear regression model with censored data.The Annals of Statistics，1990，18（1），303-328.

［23］Tsiatis AA.Estimating regression parameters using linear rank tests for censored data.The Annals of Statistics，1990，18（1），354-372.

［24］Lai TL，Ying Z.Rank regression methods for left-truncated and right censored data.The Annals of Statistics，1991，19（2），531-556.

［25］Murphy SA，Rossini AJ，Van der Vaart AW.Maximum likelihood estimation in the proportional oddsmodel.Journal of the American Statistical Association，1997，92（439）:968-976.

［26］Lin DY，Ying Z.Semiparametric analysis of the additive riskmodel.Biometrika，1994，81（1）:61-71.

［27］Zeng D，Cai J.Asymptotic results for maximum likelihood estimates in joint analysis of repeated measurements and survival time.The Annals of Statistics，2005，33（5）:2132-2163.

［28］Flem ing TR，Harrington DP.Counting process and survival analysis.New York:John W iley＆Sons，1991.

［29］Aalen O.Nonparametric inference for a family of counting processes.The Annals of Statistics，1978，6（4）:701-726.

［30］Hosmer DW，Lemeshow S，May S.Applied survival analysis:regression modeling of time-to-event data，Second Edition.New York:John Wiley＆Sons，2008.

［31］高峻，董偉，高爾升，等.多結局生存分析模型與Cox模型的隨機模擬比較.中國衛(wèi)生統(tǒng)計，2007，24（3）:248-254.

［32］Kelly PJ，Lim LL.Survival analysis for recurrentevent data:an application to childhood infectious disease.Statistics in Medicine，2000，19（1）:13-33.

［33］Eric WL，Wei LJ，Amato DA，et al.Cox-type regression analyses for large numbers of small groups of correlated failure time observations.Survival analysis:state of the art.Netherlands:Springer Netherlands，1992:237-247.

［34］Finkelstein DM，Schoenfeld DA，Stamenovic E.Analysis of multiple failure time data from an AIDS clinical trial.Statistics in Medicine，1997，16（8）:951-961.

［35］Thomas L，Reyes EM.Tutorial:survival estimation for Cox regression models with time-varying coefficients using SAS and R.Journal of Statistical Software，2014，61:1-23.

［36］Powell TM，Bagnell ME.Your“survival”guide to using time-dependent covariates.SASGlobal Forum，2012:168.

［37］Zeng D，Lin DY.Efficient estimation for the accelerated failure time model.Journal of the American Statistical Association，2007，102（480）:1387-1396.

［38］Lin DY，Wei LJ.Accelerated failure time models for counting processes.Biometrika，1998，85（3）:605-618.

［39］Huang Y，Peng L.Accelerated recurrence time models.Scandinavian Journal of Statistics，2009，36（4）:636-648.

［40］Sparling YH，Younes N，Lachin JM.Parametric survival models for interval-censored data with time-dependent covariates.Biostatistics，2006，7（4）:599-614.

［41］Ham id HA.Flexible parametric survival models with time-dependent covariates for right censored data.University of Southampton，2012.

國家自然科學基金（81460519，H2611）；云南省自然科學基金（2013FZ064）

1.昆明醫(yī)科大學公共衛(wèi)生學院流行病與衛(wèi)生統(tǒng)計學系（650500）

2.復旦大學公共衛(wèi)生學院生物統(tǒng)計學教研室

△通信作者:趙耐青，E-mail:nqzhao1954@163.com

（責任編輯:郭海強）