肖媛媛 許傳志 趙耐青

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

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

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

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

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

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

非參數(shù)法

1.常用方法

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

（1）Kaplain-Meier法

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

這一取值也被稱作K-M估計值（K-M estimator）。不難看出，K-M估計值從定義來看在時間維度上應(yīng)該是一個離散函數(shù)。如假設(shè)所有死亡事件均恰好發(fā)生在死亡時點，而兩個離散的死亡時點間無死亡發(fā)生，則可根據(jù)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）以秩次為基礎(chǔ)的非參數(shù)檢驗

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

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

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

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

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

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

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

半?yún)?shù)法

1.常用方法

（1）半?yún)?shù)比例風(fēng)險模型（Cox proportional hazards model）

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

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

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

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

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

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

（2）半?yún)?shù)加速失效模型（sem i-parametric accelerated failure timemodel）

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

模型的表達(dá)式一樣不失簡潔:

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

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

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

（3）其他半?yún)?shù)生存分析模型

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

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

（1）Cox比例風(fēng)險模型

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

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

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

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

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

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

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

（2）半?yún)?shù)加速失效模型

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

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

參數(shù)法

1.常用方法

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

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

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

參數(shù)比例風(fēng)險模型和半?yún)?shù)的Cox比例風(fēng)險模型在表達(dá)式上基本相同:

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

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

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

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

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

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

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

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

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

總 結(jié)

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

［1］林靜.基于MCMC的貝葉斯生存分析理論及其在可靠性評估中的應(yīng)用.南京理工大學(xué)，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ù)據(jù)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模型在生存分析中的應(yīng)用.中國衛(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］高峻，董偉，高爾升，等.多結(jié)局生存分析模型與Cox模型的隨機(jī)模擬比較.中國衛(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.

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

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

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

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

（責(zé)任編輯:郭海強(qiáng)）