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

肖媛媛 許傳志 趙耐青

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

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

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

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

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

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

非參數(shù)法

1.常用方法

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

（1）Kaplain-Meier法

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

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

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

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

（2）以秩次為基礎(chǔ)的非參數(shù)檢驗(yàn)

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

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

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

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

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

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

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

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

1.常用方法

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

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

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

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

P（i對(duì)象在 t（j）時(shí)刻死亡｜在 t（j）時(shí)刻存在風(fēng)險(xiǎn)的R（j）個(gè)對(duì)象至少有一個(gè)死亡）

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

參數(shù)法

1.常用方法

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

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

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

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

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

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

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

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

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

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

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

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

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

總 結(jié)

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

［1］林靜.基于MCMC的貝葉斯生存分析理論及其在可靠性評(píng)估中的應(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］陳靖，何春拉，潘建紅，等.無(wú)刪失生存數(shù)據(jù)Wilcoxon秩和檢驗(yàn)與logrank檢驗(yàn)的比較.中國(guó)衛(wèi)生統(tǒng)計(jì)，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)用.中國(guó)衛(wèi)生統(tǒng)計(jì)，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ī)模擬比較.中國(guó)衛(wèi)生統(tǒng)計(jì)，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.

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

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

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

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

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