治癒模式之文獻回顧

(1)

國立交通大學

統計學研究所

碩士論文

治癒模式之文獻回顧

Literature Review of the Cure Model

研究生：周欣茹

指導教授：王維菁教授

(2)

治癒模式之文獻回顧

Literature Review of the Cure Model

研究生：周欣茹 Student：Hsin-Ru Chou

指導教授：王維菁教授 Advisor：Dr. Weijing Wang

國立交通大學統計學研究所碩士論文 A Thesis Submitted to Statistics College of Science

National Chiao Tung University

in partial Fulfillment of the Requirements for the Degree of

Master in Statistics

June 2004

Hsinchu, Taiwan, Republic of China

(3)

摘要

傳統存活分析假設感興趣的事件一定會發生，然而當事件由“死亡＂推廣到如“發病＂、“復發＂… 的情形時，這個假設便不盡合理。文獻修正的方向是承認存在所謂的“免疫者＂永遠也不會發生此事件，這類方法統稱為“治癒模式＂。本論文回顧治癒模式的部份文獻，將其分成三個方向討論。最常見的治癒模式對免疫者無明確的定義；有一類模式以存活是否超過某固定時間做為免疫的指標；另一方向則以競爭風險發生順序為指標，後兩者對免疫均有明確的定義。所考慮的推論方法，則有迴歸分析、非迴歸分析以及無母數分析三種，論文的重點放在迴歸分析方法的回顧。

(4)

ABSTRACT

In classical survival analysis, it is implicitly assumed that the event of interest will occur eventually. However, this assumption may not implausible when there exist a proportion of subjects who will never experience the event despite of long-term follow-up.

In the thesis, we will review the literature on cure models. Three types of cure models are considered according to the definition of “immune” or “cure”. In the first type, there is no explicit definition for the immune. In other words, cured individuals are always mixed with susceptible but censored ones. The second class defines cure as being able to survive beyond a pre-specified time period. For the third type, whether a subject is immune is determined by the order of competing events. Inference methods for cure models include parametric, semi-parametric and nonparametric analysis. The focus here is on semi-parametric regression analysis.

(5)

致謝

這篇論文的完成，代表著我的學生生涯暫且告一個段落。將近二

十年的求學過程中，有快樂、有滿足，也有低潮。在這兩年裡，經歷

了許多，在此要感謝的人很多。

首先，要感謝我的指導老師王維菁老師，謝謝她這一年來在課業

上和生活上對我的指導與幫助，尤其在我經歷最低潮的時候，老師所

給予我的協助，使我能夠順利地度過，才有這篇論文的產生；感謝所

上每位教授在這兩年來的教導以及學長姐、郭姐和靜怡這段時間以來

對我們的照顧，還有班上同學們在日常生活上相互關心與支持；最後

要感謝我的家人，我們一起度過了最艱難的時期。

即將邁入人生中的另一階段，未來在新的環境底下，希望能將所

學有所發揮。

周欣茹僅誌中國民國九十三年仲夏于交通大學統計所

(6)

目錄

第一章序論 01

第二章第一類治癒模式-無法直接辨識免疫者 03

2.1 Logistic / Weibull 模式 04

2.2 Logistic / Cox 模式 05

1. Kuk and Chen 的邊際概似函數估計法 07

2. EM 演算法 07 2.3 無母數分析 11 2.3.1 充分追蹤時間的條件 12 2.3.2 Kaplan-Meier 估計量在尾端的性質 13 2.3.3 Kaplan-Meier 估計量尾端會降到 0 的機率 14 2.3.4 Kaplan-Meier 估計量的分配理論與收斂速度 15 第三章以固定時間為切點之治癒模式 17 3.1 Jung﹝1996﹞ 17 3.2 Subramanian (2001) 18 第四章考慮競爭風險之治癒模式 19 4.1 Larson-Dinse 的迴歸分析 20 4.2 Logistic/Weibull 模式 -- Taylor 24 4.3 母數分析─非迴歸模式 25 4.4 無母數分析 – Wang 27 第五章結論 30 附錄一、參考文獻 32 附錄二、治癒模式文獻一覽 35

(7)

第一章序論

傳統的存活分析假設感興趣的事件一定會發生，可以參考近年出版的教科書包含 Anderson et al.﹝1993﹞、Hougaard﹝2000﹞、Klein and Moeschberger ﹝1997﹞…等。當感興趣的事件由“死亡＂拓展到“發病＂、“復發＂ … 時，傳統存活分析的假設便不盡合理。因此文獻上出現所謂“治癒模式＂﹝cure model ﹞ 承認只有部分的人會發生感興趣的事件，稱之為 “ 可致病＂ ﹝susceptible﹞；但有部份的個體永遠不會發生此事件，稱之為“免疫＂ ﹝immune﹞。許多治癒模式的文獻採用混合模式﹝mixture model﹞做為分析的架構。混合模式的做法是將母體區分為兩類：免疫﹝immune 或 cure﹞和可致病 ﹝susceptible﹞。令T為到感興趣事件的發生時間，定義指標函數B，當B=1 代表可致病，此時T <∞；當B=0代表免疫，此時永遠觀察不到事件的發生，那麼可令T =∞。可以將T的分佈拆解為 ) 0 Pr( ) 1 Pr( ) 1 | Pr( ) Pr(T >t = T >t B= B= + B = 。上式的右邊將T 的分佈分成兩個部分：一部分是Pr(B=1)，稱之為致病模式

﹝cumulative incidence model﹞；另一部分為Pr(T > Bt| =1)，稱之為潛在發

病時間模式﹝latency model﹞。可發現當免疫者存在時，T的分佈並不符合一般隨機變數的機率性質，因為 0 1 ) 0 Pr( ) Pr( lim > = = = − > ∞ → T t B p t 。我們可令為 cure rate，稱之為治癒率或是免疫率，這個參數則是許多文獻感興趣的量。 p − 1 統計推論的目的是根據所觀測的資料，對感興趣的參數做推論，內容包含估 計和檢定…。若是在有限的時間點 t 已觀察到感興趣的事件發生，則可確知B=1 且。但是若是觀察不到事件發生，則有兩個可能情形 ─ 一是永遠免疫或是治癒，機率為；另外一個只是暫時尚未發生，機率為。文獻中最常見的治癒模式﹝稱之為第一類模式﹞並未直接定義如何區辨免疫的 t T < ) 0 Pr(B= Pr(T > Bt, =1)

(8)

觀測值和暫時設限的觀測值。第二類型的治癒模式則是將“治癒與否＂定義為是否存活超過一個固定的年限(假設M 年)，因此治癒率就是超過M 年的存活機率。第三類型的治癒模式把治癒與否定義為競爭風險的發生種類或是發生次序。我們將在第二章到第四章以這三個方向介紹不同類型的治癒模式，討論的重心在迴歸分析。混合模式所提供的架構可供研究者探討解釋變數對於“免疫與否＂和“發病時間＂是否存在不同的影響。在第五章我們提出對整篇論文的心得以及未來可能的研究方向。 ) Pr(T >M

(9)

第二章: 第一類治癒模式-無法直接辨識免疫者

第一類治癒模式只假設能觀察到部份有致病可能﹝susceptible﹞的個體，但無法由設限資料中分辨誰是免疫者。令發生感興趣事件發生的時間為T，但是只有部份的個體會發生此事件，以指標函數B=1代表這類人;免疫者以表示，其發生事件的時間可令為 0 = B ∞ = T 。在設限的情形下存在設限時間 C ，使得只 能觀察到X =T ∧C與指標函數δ = I(T ≤C)。假設樣本為隨機樣本，右設限的資料可表示為{(X_i,δ_i)(i=1,...,n)}，可知若δ_i =1則B_i =1;然而若δ_i =0則不知的值。然而若 i B i i C X = 的值很大，則當δ_i =0時可以猜測B_i =0，這個概念和之後會定義的“充分追蹤時間＂﹝sufficient follow-up﹞有關。根據前述混合模式: ) 0 Pr( ) 1 Pr( ) 1 | Pr( ) Pr(T >t = T >t B= B= + B= 。因為 lim_t_→_∞Pr(T >t|B=1)→0，所以lim_t_→_∞Pr(T >t)→Pr(B=0)。因為T的樣本符合傳統右設限資料，Maller and Zhou﹝1992﹞提議以無母數 Kaplan-Meier 估計量估計Pr(T >t)，所建議的無母數估計量可以表示為

∏

∑

≤ = = ≥ = = − = > t u n i i n i i i u X I u X I t T } ) ( ) 1 , ( 1 { ) r( Pˆ 1 1 δ 。以做為的估計量的優點早是存活分析領域被熟知的結果，然而以做為的估計量的合理性﹝ 為最大觀測到發生事件的值﹞，卻需要特殊條件才會成立。Maller 與 Zhou﹝1992﹞的文章強調此估計量只有在充分追蹤時間﹝sufficient follow-up﹞成立時才會合理，否則會有高估的情形。他們更在 1994 年的論文中提出檢驗充分追蹤時間是否成立的檢定方法，更在 1996 年的專書中將第一類治癒模式做有系統的討論。我們將在 2.3 節中摘錄 Maller & Zhou﹝1996﹞書中的部份結果。Li et al.﹝2001﹞提及在第一類模式下，以無母數方法處理免疫問題因為缺乏額外資訊一般有不可辨識 ) r( Pˆ T >t Pr(T >t) ) r( Pˆ T >t_max Pr(B=0) t_max ) 0 Pr(B=

(10)

的問題。Wang﹝2004﹞則具體指出當死亡這個無法避免的競爭風險存在時，無母數的估計方法有不可辨識性的問題。

當存在有解釋變數時，研究的重點在探討解釋變數對於“免疫與否＂和“發

病時間＂是否存在不同的影響。令Z為解釋變數，迴歸模式架構於 incidence 的

部份表示為： ;於 latency 部份可表示為: 。通常

以 binary regression 模式描述，如 logistic regression；針對的模式，有的作者給定母數分配的假設，有的採半母數模式，有的則以無母數方法分析之。例如 Farewell﹝1982﹞利用 logistic/Weibull 做為兩部分的模式假設。有數篇文章則假設 logistic/Cox 模式，如 Kuk & Chen ﹝1992﹞，Sy & Taylor﹝2000﹞，Peng & Dear﹝2000﹞。我們在以下數節中先介紹迴歸模式的推論方法，之後再討論無母數的推論問題。 ) | 1 Pr(B= Z Pr(T >t|B=1,Z) ) | 1 Pr(B= Z ) , 1 | Pr(T >t B= Z 2.1 Logistic/Weibull 模式

Farewell 在 1982 年的文章中針對致病模式的部分以 logistic regression 模式描述如下： ) exp( 1 ) exp( ) | 1 Pr( ) ( _' ' i i i i i Z Z Z B p β β β + = = = ，對於潛在發病模式的部分則是做 Weibull 的母數分配假設， ] ) ( exp[ ) , 1 | Pr(T >t_i B_i = Z_i = − λt_i γ 。值得注意的是 Farewell﹝1982﹞的論文假設Pr(T >t_i |B_i =1,Z_i)與解釋變數無關，但是到了 1986 的論文則將解釋變數的影響納入 Weibull 分配的參數之中。根據傳統治癒模式的資料型態{(δ_i,X_i,Z_i),i=1,...,n}，可以得到概似函數為

∏

= = − ₋ = n i I i i i i i F i x x p z x L 1 ) 1 ( 1 )] ) ( exp( ) ( ) ( [ ) , , , , (λ γ β β γλ λ γ λ γ δ } ))] ( 1 ( ) ) ( exp( ) ( [ − + − ( =0) × I i i i i x p p β λ γ β δ 。因為取了對數之後再對參數微分計算太過複雜，所以利用 EM 演算法來簡化計

(11)

算。首先建構“完整資料＂{(B_i,δ_i,X_i,Z_i)(i=1,...,n)}之下的概似函數：

∏

= = − ₋ = n i I i i i C i x x p L 1 ) 1 ( 1 )] ) ( exp( ) ( ) ( [ ) , , (λ γ β β γλ λ γ λ γ δ } ))] ( 1 [( )] ) ( exp( ) ( [ − ( =0, =1) − ( =0, =0) × i i I i Bi i B I i i x p p β λ γ δ β δ 。在 E-step 中對logLF(λ,γ,β)取條件期望值﹝在給定觀測值下﹞，其中牽涉到以下等式的計算： )] , , ( | ) 1 , 0 ( [I _i B_i _i x_i z_i E δ = = δ =Pr(B_i =1|δ_i =0,T_i >x_i,z_i) ) | , 0 Pr( ) | , 0 , 1 Pr( i i i i i i i i i z x T z x T B > = > = = = δ δ ) , 1 | Pr( ) | 1 Pr( ) | 0 Pr( ) , 1 | Pr( ) | 1 Pr( i i i i i i i i i i i i i i z B x T z B z B z B x T z B = > = + = = > = = ) ) ( exp( ) ( )] ( 1 [ ) ) ( exp( ) ( γ γ λ β β λ β i i i i i x p p x p − + − − = 。此外可得 )] , , ( | 0 , 0 ( [I _i B_i _i x_i z_i E δ = = δ =Pr(B_i =0|δ_i =0,T_i >x_i,Z_i) ) ) ( exp( ) ( )] ( 1 [ ) ( 1 γ λ β β β i i i i x p p p − + − − = 。

第二個步驟 M-step 為對 E[logL_C(λ,γ,β|data)] 取極大值，可將 )] | , , ( [logL data E _C λ γ β 對參數微分令其等於 0 求解以得估計值，這部分可以利用牛頓法以求得數值解。值得一提的是E[I(δ_i =0,B_i =1)]與E[I(δ_i =0,B_i =0)]會牽涉到未知參數，估計時為代入前階段之參數之估計值，在求極值的過程中視為定值。利用母數分配的假設給定致病模式以及潛在發病時間模式，在模式正確的情形下，所得之估計量具有效性，但伴隨的問題是模式假設可能不夠一般化，使得應用性可能不夠廣；此外若是假設錯誤，所得之估計量可能有很大的偏誤因此不夠穩健。所以檢驗模式適合度有其必要性。 2.2 logistic/Cox 模式 Farewell﹝1982﹞在致病時間模式上，並未將解釋變數的影響考慮其中，這

(12)

是一大缺失。為了彌補這個問題，後續文章除了於致病模式做同樣的 logistic 模式假設外， ) exp( 1 ) exp( ) | 1 Pr( ) ( _' ' i i i i i Z Z Z B p β β β + = = = ；

致病時間模式則是利用“等比例風險模式＂(Cox proportional hazards model) 做為假設： } ) ( exp{ ) , 1 | Pr( 0 ' 0

∫

− = = > i i t z i i i t B Z h u e du T α ，其中 h₀(t) 為 B=1 群體發病時間的基底風險函數 ﹝ baseline hazard

function﹞。值得一提的是在加入了免疫群體的 Cox model 下，解釋變數對合併的群體就不再是成比例(proportional)的影響。此外，將迴歸分析用在混合模式的好處是可將解釋變數在“發病與否＂和“發病時間＂的影響分開討論，例如有的解釋變數只對發病與否有影響，有的則可能改變發病時間。根據資料型態以及 logistic/Cox 模式的假設，概似函數可以寫成: ) 1 ( 1 ₀ ' 0 ' 0 0) [ ( ) ( ) exp( ( ) )] , , ( = =

∏

−

∫

= i i i i I n i x z z i i F h p h x e h u e du L δ α α β β α (2.1) } ))] ( 1 ( ) ) ( exp( ) ( [ 0 ' ( 0) = − ₊ ₋ − ×

∫

i i i I i x o z i h u e du p p β α β δ 。 (2.2) 一般 MLE 的求解方法是將概似函數取對數以後再予以微分求解。由上式可發現，概似函數L_F(α,β,h₀)取對數後，(2.2)式變成為 )] ) ( 1 ( ) ) ( exp( ) ( log[ ) 0 ( 1 0 ' 0

∑

∫

= − ₊ ₋ − × = n i i x z i i p h u e du p I i i β β δ α _。其中 log 裡為連加的情況，一旦再對參數微分後，整個計算將會變得複雜許多。問題所在來自當δ_i =0 時，概似函數必須傳遞兩種可能性 ﹝ 即和 ﹞，使得函數有連加的部份。可發現若是的值為已知，連加的部份會變為連乘，取對數後函數變為線性，求極值時計算亦得以簡化。以上的想法就是 0 = i B 0 = i B B_i

(13)

EM 演算法的概念：先假設完整資料以建構簡單的概似函數，對於缺失值的部份以條件期望值估計，再求極值。先前討論的 Farewell 就是以此想法提出估計方法，當模式拓展到半母數的 Cox 模式時，在 E-step 時會牽涉到未知的基底函數，推論的難度在於估計 ) ( 0 t h (α,β)時如何處理無窮維度的未知函數。對於 Cox 模式的推論﹝不包含免疫者﹞，一般以部份概似函數﹝partial likelihood function﹞估計，文獻發現雖然忽略h0(t)的影響但所得之α 估計量仍具有相當不錯的表現。然而當有免疫者存在時，若忽略基底風險函數的資訊，則會對 ) ( 0 t h α 的估計造成偏誤，代表部份概似函數估計法不可行。因此針對 logistic/Cox 模式有關致病時間模式部份的估計，便面臨到要估計基底風險函數的挑戰。文獻提出了數種方法針對基底風險函數的處理，我們在此簡述其概念。

2.2.1 Kuk and Chen 的邊際概似函數估計法

Kuk and Chen﹝1992﹞針對 logistic/Cox 模式估計的部分，考慮以邊際概似函數﹝marginal likelihood﹞的方法估計感興趣的迴歸參數(α,β)，這部份可以不需處理基底風險函數。概念簡述如下：根據所得的右設限資料型態 ) , , {(δ_i X_i Z_i ，，可將其分成兩部分，一部分為有設限的觀測值，另一部分為沒有設限的觀測值。將資料根據發生事件的時間做排序，再將排序資料做多重積分後建構邊際概似函數，這個做法可以巧妙的將基底函數去除，因此迴歸參數 n} 1,..., i= (C) ) (D ) , (α β 的估計量為對邊際概似函數求極值獲得，無須在此階段處理基底函數估計。然而這個方法的難度在於多重積分即使經過簡化，仍牽涉大數目的排列組合個數，使得概似函數變成是龐大組合個數的連加，因此作者建議以 Monte Carlo 模擬法求概似函數的近似值後再求解。後續的論文在估計(α,β) 時，選擇寧可直接面對基底函數存在的事實，而不願意為捨棄這個項目而引來更為複雜的計算問題。 2.2.2 EM 演算法

(14)

在 logistic/Cox 模式下的 EM 演算法的概念簡述如下：假設資料型態可表示為{(δ_i,X_i,z_i,B_i),i=1,...,n}，此時概似函數可表示為

∏

= = = ₋ = n i B I i B I i C i i _p p h L 1 ) 0 ( ) 1 ( 0) { ( ) (1 ( )) } , , (α β β β (2.3)

∏

∫

= = = = ₋ × n i B I x o z B I z i i i i i i i _h _u _e _du e x h 1 ) 1 ( ' 0 ) 1 , 1 ( ' 0( ) ] [exp( ( ) )] } {[ α δ α 。 (2.4) 則取對數以後可得，

∑

= − = + = = n i i i i i C h I B p I B p L 1 0) { ( 1)log ( ) ( 0)log(1 ( ))} , , ( log α β β β

∑

₌ = = + + = −

∫

+ n i x z i i i i i i i_du e u h B I z x h B I 1 0 ' 0 0( ) ' ] ( 1)[ ( ) ]} )[log 1 , 1 ( { δ α α 。可發現在logL_C(α,β,h₀)中的 log 函數裡不再有相加的問題，求極值的過程中對參數微分時也較好計算。和前述 Farewell﹝1982﹞的做法相同：E-step 須對設限的個體求其得病可能性的期望值來作為權數﹝weight﹞，推導如下： ] , , 0 | ) 1 ( [ ) 1 ( )] , , ( | ) 1 ( [I B_i _i x_i z_i I _i E I B_i _i T_i x_i z_i E = δ = δ = + = δ = > ，在 Cox 模式假設下可得 ] , , 0 | ) 1 ( [I B_i _i T_i x_i z_i E = δ = > =Pr(B_i =1|δ_i =0,T_i >x_i,z_i) ) | , 0 Pr( ) | , 0 , 1 Pr( i i i i i i i i i z x T z x T B > = > = = = δ δ ) , 1 | Pr( ) | 1 Pr( ) | 0 Pr( ) , 1 | Pr( ) | 1 Pr( i i i i i i i i i i i i i i z B x T z B z B z B x T z B = > = + = = > = = ) ' exp( 0 ) ' exp( 0 ) 1 | ( ) ( )] ( 1 [ ) 1 | ( ) ( i i z i i i i z i i i B x S p p B x S p α α β β β = + − = = ， E[I(B_i =0)|(δ_i,x_i,z_i)]=1−δ_i +E[I(B_i =0)|δ_i =0,T_i > x_i,z_i] ) ' exp( 0( | 1) ) ( )] ( 1 [ ) ( 1 1 i z i i i i i i B x S p p p α β β β δ = + − − + − = 。如前所述在 Cox 模式下這個權數包含未知的基底存活函數，需要額外處理。我們將在之後敘述兩篇文章的處理方法。EM 演算法的第二個部分為 M-step，將對數概似函數的期望值取極大值，其中權數部分經代入前階段的參數估計值可視為定

(15)

值，如此可使求極值的步驟不致太複雜。以下簡介兩篇文章以不同方式處理基底函數的估計問題。 (1) EM 演算法 -- Sy and Taylor Sy and Taylor﹝2000﹞的文章中提出兩種方法估計基底風險函數。 i. Breslow-type 估計量：令H₀(t)為基線累積風險函數，其定義為

∫

= t du u h t H 0 0 0( ) ( ) 。以 Breslow 的方法估計H0(t)如下:

∑ ∑

≤ ∈ ⎟⎟ ⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎜ ⎝ ⎛ = = t t i R l z l i i i l e w d B t H ) ( : ' 0( | 1) ~ α ，其中d_i為在t時間事件發生的個數，分母的部分則是在t時間所有處於風險之個體的加權期望機率，此加權的權數為風險比例eα' i w 為個體被觀察到得病的可能性，若為非設限資料， i z _， 1 = i w ；若為設限資料，w_i =π_i， ) 1 _i _i _i _i i T x z B I = = > 即 [ E ( |δ 0, , ] 。基線存活函數S₀(t|B=1) 的估計量可由 ) 1 | ( ~ 0 t B= H 求得如下: )) 1 | ( ~ exp( ) 1 | ( ~ 0 0 t B= = −H t B= S 。換句話說，Breslow-type 估計量是將未知的基底風險函數利用未知參數α 和來做表示，藉此降低未知參數的維度以簡化估計，代價是在 E-step 裡的權數是一個非常複雜的 i w i w α 和β的函數。在 M-step 中依然視為定值﹝代入前階段 i w α 、β 估計值﹞，利用 Cox 建議的部分概似函數將(2.4) 式簡化後，取極大值，經過疊代的過程，所收斂的解即為參數估計量。 ii.Product-limit 估計量：前述 Breslow 的方法將未知的表示成為 ) 1 | ( 0 t B= S α 與β 的“explicit＂函數，以此降低未知參數的維度，卻同時增加函

(16)

表示成風險機率參數的連乘，對這些額外的參數再以 NPMLE 的概念求解，其解可表示為α 、β 的函數。其概念簡述如下：根據 Kalbfleisch and Prentice(1980)的想法，將資料分為設限與非設限兩部分，(2.4) 式可以寫成 ) (C_i (D_i)

∏

= ∈ = ∈ − ₌ _× ₌ k i l C z B I i D l z i i i i l l i l _S _t _B B t S z t h 0 ) ' exp( ) 1 ( ) ( 0 ) ' exp( ) ( 0 ) ( ; ) ( | 1) ] ( | 1) } ( [ { α α 。 (2.5) 根據 PH 模式，基線存活函數可以 Product-limit 的形式表示為

∏

≤ = = t t j j j B t S ) ( : 0( | 1) γ ， (2.5)式即可利用 Product-limit 形式改寫為

∏

= ∈ ∈ − × − k i l R D z w i D l z i i i l l i l 0 ( ) ) ' exp( ) ' exp( } ] 1 [ { γ α γ α ， (2.6) 此處的定義與 Breslow-type 估計量中所提到的相同。換句話說，經過以上的整理，概似函數表示成為 l w α 、γ_i(i=1,...,k)與β的函數，其中γ_i為 “failure＂的個數。要“同時＂﹝simultaneously﹞針對這些函數求極值是似乎是難以達成的任務，因此 Sy and Taylor 建議以下的方式求解：固定α 、γ_i的估計函數可表示為 k i e w e i l i i l R l z l D l z i z ,..., 1 , 1 ' ) ' exp( ' = = ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ −

∑

∈ ∈ α α α γ 。 (2.7) 將以上具有 explicit-form 的 γ_i估計量後再代入概似函數(2.6)式中求 α。中間權數w_l部分，牽涉到α、β 與γi，經過層層疊代可以求得估計量。

(2) EM 演算法--Peng and Dear

Peng and Dear﹝2000﹞提出的做法與 Sy and Taylor﹝2000﹞的方法類似，兩者都是先根據完整資料建立對數概似函數，求其條件期望值後再求極值，兩篇論文的差異在 M-step 的處理三種參數α、β 和S₀(t|B=1)的方法不同，此外疊代步驟亦有差異。在 M-step 的部分，Peng and Dear 的做法和前述 Sy and Taylor

(17)

所提 Breslow-type 估計量類似，S₀(t|B=1)化簡為同樣的α 與β 的函數。但因為 (2.4)式較難處理，因此根據 Breslow 的做法，(2.4)式可寫為

∏ ∑

= ∈ k j R i d i i j j j z w s 1{ exp( ' )} ) ' exp( α α ， (2.8) j s 為所有非設限資料之解釋變數的總合， 為在時間事件發生的個數。M-step 的部分被簡化成對(2.8)式取極大值，(2.8)式中只包含了 j d t α 這個未知參數，並不含未知的基線風險函數，處理上來也較容易。

在 E-step 的部分，Peng and Dear 的做法則是和 Sy and Taylor 所提 Product-limit 估計量類似，根據 M-step 所估計出來的α、β以及 Kalbfleisch and Prentice(1980)的想法，(2.4)式可寫成(2.6)式的形式，可以得到γ _j ) ,..., 1 (j= k 的估計量

∑

∈ − ≈ j R i i i j j z w d ) ' exp( 1 ˆ α γ ，進而可以估計基底存活函數，將帶入 E-step 中的，如此疊代數次後即可求的估計量。 ) ( 0 t S S0(t) wi

2.3 無母數分析

此類的治癒模式的資料和右設限模式毫無區別，然而選擇以傳統存活分析的方法或是以治癒模式的角度卻需要事先做決定。統計分析所牽涉到的問題在於分析方法需要做額外的假設，先前提到的方法是依靠模式的假設，在此節中所提到的無母數分析則需要靠資料品質的假設，才能夠確保估計值的唯一性以避免所謂 “不可辨識＂﹝non-identifiability﹞的問題。Farewell﹝1986﹞建議研究者以應用領域的專業判斷免疫者存在與否做為抉擇的準則。無母數分析經常使用 Kaplan-Meier 估計量以估計Pr(T >t):

(18)

∏

∑

≤ = = ≥ = = − = > t u n i i n i i i u X I u X I t T } ) ( ) 1 , ( 1 { ) r( Pˆ 1 1 δ 。文獻中討論 Kaplan-Meier 估計量性質的文章相當多。我們在此節中整理 Maller and Zhou ﹝ 1996 ﹞ 書中的結果，討論在治癒模式的假設下使用

以估計 ) r( Pˆ T >tmax Pr(B=0)=1−p的可行性。經常使用存活分析方法的人都知道 K-M 估計量在估計時間尾端的分配時表現並不穩定。原因之一是當值很大時的個數變小，這部份位居分母，使得估計的變異性增大。另一個更為相關的問題在於的範圍只能估計限於能觀測到發生事件的最大時間點，如果觀察時間太短，會有一大部份本來是 susceptible 的人受到設限，此時若以估計會造成嚴重的高估。在下一小節中我們以較為嚴謹的數學表示法來說明所謂“充分追蹤時間＂﹝sufficient follow-up﹞的問題。 t

∑

= ≥ n i i t X I 1 ) ( ) r( Pˆ T >t ) r( Pˆ T >tmax 1−p 2.3.1 充份追蹤時間的條件

仿照 Maller and Zhou (1996)，定義以下幾個端點: } 1 ) Pr( : { inf ≤ = = _t t X t H τ ， } 1 ) Pr( : { inf ≤ = = _t t C t G τ ， } 1 ) Pr( : { inf ≤ = = _t t T t F τ 。假設觀察時間有限，τ_G <∞。在右設限下X =T∧C，可知τ_F ∧τ_G =τ_H。當不存在免疫者時，我們希望τ_G >τ_F，這樣τ_F =τ_H 才能有機會觀察到最大可能的事件發生時間。若τ_G <τ_F，代表觀察時間太短，會有 heavy censoring 出現，此時 K-M 曲線不會降到 0。然而當免疫者存在時，免疫者的發病時間定義為T =∞，此時，T 分配稱之為“improper＂，無論觀測時間如何拉長都會有 _G 1 ) Pr( ≤∞ < = T p 的 ∞ = ≤τF τ 。此時可以定義 } 1 ) 1 | Pr( : { inf 0 = t t T >t B= = τ ，

(19)

代表可觀察到 susceptible 者最大可能發生事件的時間。此時我們反而希望

0

τ

τ_G > ，代表觀察時間充裕到足以使所有 susceptible 的人都有足夠時間發生事件。

Maller and Zhou ﹝1996)﹞在書中的第二章先將如何檢驗“免疫者是否存在＂和“觀察時間是否充足＂量化成數學的符號。之後再討論以 K-M 估計量做為檢驗方法的可行性。第一步: 檢驗H₀₁:Pr(T ≤τ_G)=1。若是接受，代表觀察時間足以包含所有發生時間，可知 01 H ∞ < F τ 而且免疫者不存在，觀察時間亦充份。然而若是被拒絕，無法得知究竟是免疫存在使得 01 H ∞ = F τ 或者是觀察時間不夠充份，而使得τ_G的值太小。此時我們要進行第二步驟的討論。第二步: 在已知Pr(T ≤τ_G)<1的情形下，檢驗H₀₂:τ₀ ≤τ_G。若是成立，代表觀測時間夠長到足以使得所有 susceptible 的個體都發生事件。此時可以確認免疫的族群存在即。若被拒絕，代表研究時間不夠充份﹝follow-up

is not sufficient﹞，無法由資料判別免疫者是否存在。此時如 Farewell 所言，研究者依自己的專業判斷免疫者是否存在，分析時避免以無母數方法而是靠模式假設做分辨“真正免疫者＂和“暫時設限者＂的工作。 02 H 0 1− p> H₀₂ 2.3.2 Kaplan-Meier 估計量在尾端的性質

以下內容摘錄自 Maller and Zhou﹝1996﹞書中的第三章，我們討論如何檢驗前節所提出的兩個假說。討論估計量Pˆr(T >t_max)的性質可分為以下幾個方向:

a. Theorem 3.2: t_max →τ_H almost surely。

在 K-M 的估計式中t_max被用來估計τ ，很顯然 _H t_max ≤τ_H。這個定理說明當樣本數很大時，所觀察到的值會很接近真正被估計的端點。 b. Theorem 3.4: 當Pr(T ≤t)在τH 是連續的情形下， in probability。 ) Pr( ) r( Pˆ T ≤tmax → T ≤τH

(20)

這個結果保證了Pˆr(T ≤tmax)做為Pr(T ≤τH)的點估量具有一致性。值得注意的是最終目標是以估計的可行性，所以我們要討論與 在 t 發生於不同界限時的關係。 ) r( Pˆ T ≤tmax p p ) Pr(T ≤t c. Theorem 3.5: 當τ_G ≤τ_F ，則τ_G =τ_H，此時 in probability。 ) Pr( ) r( Pˆ T ≤t_max → T≤τ_G 當τ_G ≤τ_F 觀察時間未能包含最大可能觀測到的發病時間，值得一提的是這個結果未討論免疫者存在的情形。當免疫者不存在時，τ_G ≤τ_F 是不理想的狀態;然而當免疫者存在時，τG <τF 卻是必然，而我們希望以估計，所以有以下的結果。 ) r( Pˆ T ≤tmax 1 < p d. Theorem 3.6: 只有當τ₀ ≤τ_G，才有Pˆr(T ≤t_max)→Pˆr(T ≤τ₀)= p。這個結果說明了充分追蹤時間(τ₀ ≤τ_G)的重要性，此時以 Kaplan-Meier 尾端的估計量用來估計免疫的比例才會合理。值得注意的是若τ_G =τ_H <τ₀，則會高估，因為若將觀察時間延長會繼續記錄到事件發生，使得還有下降空間。 ) r( Pˆ T >t_max 1−p ) r( Pˆ T >t 2.3.3 Kaplan-Meier 估計量尾端會降到 0 的機率一旦，我們就會得到。當免疫者實際存在時﹝尤其當比例很小 ﹞ ，資料仍有可能得到的情形，雖然實際上。因此 Maller and Zhou 在 3.2 節的主題就以計算了在重複抽樣下

會得到的 frequency。熟悉存活分析的人都知道當最大的觀測值是觀察到的 failure 則 Kaplan-Meier 曲線會降到 0 即。令 1 ) r( Pˆ T ≤ tmax = Pˆr(B=0)=0 0 ) 0 r( Pˆ B= = 0 ) 0 Pr(B= > 1 ) r( Pˆ T ≤ tmax = 1 ) r( Pˆ T ≤ tmax = δmax代表t_max對應的指標函數，則

(21)

) 1 Pr( ) 0 ) 0 r( Pˆ Pr( B= = = δmax = 。

書中的 Theorem 3.9，推導了在p≤1的情形下Pr(δ_max =1)的理論值; Theorem 3.10，推導了在p≤1的情形下Pr(δ_max =0)的理論值。這些值都以積分或是連加的方式表示，與變數p值和T|B=1與 C 的母數分配有關。實際的計算在一般情形下需要用到數值方法。但是若分配是 exponential distribution，則可以得到 explicit 的結果。此外兩個定理的結果包含所有的樣本大小，所以是 finite-sample 的結果。 Theorem 3.11 和 3.12 則是針對τ₀,τ ,_F τ_G 大小的不同排列情形推導 ) 1 Pr(δ_max = 在大樣本的極限值。基本上我們要看當 p=1 時，是否 1 ) 1 Pr( lim_n_→_∞ δ_max = = ，因為當免疫者不存在且樣本很大時， K-M 曲線降到 0 的機率應該是 1 才合理。同理當 p<1代表免疫者存在，我們應該得到 1 ) 0 Pr( lim_n_→_∞ δ_max = = ，代表 K-M 曲線不會降到 0 的機率也應該是 1 才合理。我們重述其理論:

Theorem 3.11: 當p=1且觀測時間太短以致τ_G <τ_F，則lim_n_→_∞Pr(δ_max =0)=1。當Pr(B=1)<1，則亦會得到lim_n_→_∞Pr(δ_max =0)=1。

Theorem 3.12: 當p=1且觀察時間夠長到使得τ_G >τ_F，則lim_n_→_∞Pr(δ_max =1)=1。

書中約略提到τ_G =τ_F的情形，需要加入特殊的條件才會使得lim_n_→_∞Pr(δ_max =1)的極限存在。

2.3.4 Kaplan-Meier 估計量的分配理論與收斂速率

令。Maller and Zhou 的 Theorem 4.1 證明只有當觀察時間充份的情形下 ) r( Pˆ ˆ T t_max p= ≤ ) (τ₀ ≤τ_G ，。這個定理說明具有一致性的充分必要條件就是 p B p→p = = ) 1 Pr( ˆ pˆ G τ τ0 ≤ 。此時Pr(T ≤ Bt| =1)的無母數估計量為，在同樣的條件下具有 uniform consistency 的性質。 p t T )/ ˆ r( Pˆ ≤ p t T )/ ˆ r( Pˆ ≤

(22)

以上討論的是做為點估計量的性質，欲做後續的推論﹝如區間估計和假設

檢定﹞則需要有關於分配理論的推導結果。在 Maller and Zhou﹝1996﹞書中

第四章的 Theorem 4.1 證明了當 pˆ pˆ 1 0< p< 的情形下﹝即免疫者存在﹞， ) ˆ (p p n − 會收斂到常態分配，此外書中亦提供了Var{ n(pˆ−p)}的公式。然而當p=1，只知 n(pˆ−1)→0。換言之當免疫者不存在時﹝p=1﹞，會收斂到一個 non-degenerate distribution，只知的收斂速率比一般的 ) ˆ (p p n−q − pˆ n−1/2 速率來得快﹝即 ﹞，但是確切的值未知，而且的分配型態亦未知。文獻上稱時的性質討論為“boundary problem＂，因為是比例的上界。在一般統計分析中研究這個問題都是困難的， ﹝如泰勒展開式 …﹞在邊界無法使用。因為理論推導的問題牽涉到困難的數學分析能力，我們不做進一步討論。在此我們討論 2 / 1 > q q n−q(pˆ−p) 1 = p p=1 因為許多有用的數學工具 1 = p 的性質在實際問題的應用。點估計量pˆ的性質在所有p≤1 先前提到的第一步驟檢驗的情形都適用。但是做進一步的推論卻需要判斷是否免疫者存在。 H₀₁:Pr(T ≤τ_G)=1就包含檢定免疫者存在的情形。因為若是H₀₁被接受了，可以推得p=1; 若被拒絕則要繼續做

第二步驟的檢定。Maller and Zhou﹝1994﹞的論文提出討論檢定充份追蹤

時間的方法欲檢定，可以利用 02 H 。 H₀₁ pˆ−1的距離做為判別標準，然而的分配型態和值在卻是未知的。根據 Neyman-Pearson 法則，棄卻範圍由控制型一錯誤 ) 1 ˆ ( − − p n q q p=1 α 決定: Pr(pˆ−1>cn(α)|p=1)=α。換言之若是pˆ >1−cn(α)則拒絕，承認免疫不存在。然而 01 H c_n(α)的值卻因為pˆ的性質在p=1缺乏理論依據所以無法求得，使得以無母數的方法檢驗的目標遇到很大的挑戰。Maller and 的檢定。 01 H Zhou 以模擬分析討論H₀₁

(23)

第三章:以固定時間為切點之治癒模式

有一些文獻把存活超過某時間(M 年)定義為治癒，此時治癒比例為 ) Pr( ) 1 Pr(B= = T ≤M 。在解釋變數存在的情形下，模式建構於 incidence 部份: ) exp( 1 ) exp( ) ( ) | 0 Pr( _' ' ' Z Z Z Z B β β β π + = = = 。在沒有設限的情形下，資料可以表示為，可自動轉換成為，其中，則概似函數方法可表示為： ) ,..., (T₁ T_n (B₁,...,B_n) ) (T M I B_i = _i ≤ ) 0 ( ' ) 1 ( 1 ' ' ) exp( 1 1 ) exp( 1 ) exp( ) ( = = = ⎟⎟⎠ ⎞ ⎜⎜ ⎝ ⎛ + ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ + =

∏

i i I B i B I n i i i Z Z Z L β β β β ，對參數微分後可得 score equation： 0 ) ( ) ( ) ( } ) ( ) ( { _' _' ' 1 ' = − ≥

∑

= i i i i n i i i Z Z Z Z Z M T I β π β π β π β π φ ， (3.1) 其中 φ φ π φ π_φ ∂ ∂ = ( ) ) ( ，π = 1−π 。然而當設限存在時，我們只能觀察到{(X_i,.δ_i), ，其中 , } ,..., 1 n i= X_i =T_i ∧C_i δ_i =I(T_i ≤C_i)，代表設限變數，也因此的值有可能發生未知的情形。當 i C B_i 1 = i δ ，很顯然B_i =I(X_i ≤M)；當δ_i =0，此時的值是未知的。以此為方向的治癒模式並非我們主要的興趣，因此我們只選兩篇文章簡述其概念。 i B 3.1 Jung﹝1996﹞ Jung 發現E(I(X_i ≥t)|Z_i)=Pr(T_i ≥t|Z_i)G(t)，其中G(t)=Pr(C_i >t)。因此建議以I(X_i ≥M)/Gˆ(M)做為I(B_i =0)的“代理者＂(proxy)，修正(3.1)的估計函數可得： 0 ) ( ) ( ) ( } ) ( ) ( ˆ ) ( { _' _' ' 1 ' = − ≥

∑

= i i i i n i i i Z Z Z Z Z M G M X I β π β π β π β π φ 。 (3.2)

(24)

這個方法使用的原則是處理 missing data 常用的 “inverse probability weighting＂的方法。 3.2 Subramanian (2001) Subramanian 文章的概念是以E[I(T_i ≥M)|δ_i,X_i,Z_i]的無母數估計量(以表示之)取代(3.1)式中的。此文章所做的額外假設為解釋變數為離散變數，如此才可能針對個別可能的 i Eˆ ) (T M I _i ≥ Z_i Z值估計。所提出的估計函數可表示為 ) | Pr(X_i ≥M Z_i 0 ) ( ) ( ) ( } ) ( ˆ { _' _' ' 1 ' ₌ −

∑

= i i i i n i i i Z Z Z Z Z E β π β π β π β π φ ， (3.3) 其中 ] , , | ) ( [ ˆ ˆ i i i i i i E I T M X x Z E = ≥ δ = ) | r( Pˆ ) | r( Pˆ ) 0 , ( ) ( i i i i i i i Z x T Z M T M x X I M X I > > = < = + ≥ = δ 。 (3.4)

這個方法所應用的原則被稱為“imputation by conditional mean＂。這個方法的缺點在於Pr(T >t|Z)的估計與Z有關，所需要的資訊非的假設所能夠涵蓋。若是解釋變數為離散型，則可以根據 ) ( ) | 0 Pr(B= Z =π β'Z Z值將資料做切割，再以 Kaplan-Meier 的無母數估計量估計Pr(T >t|Z)。

(25)

第四章:考慮競爭風險之治癒模式

傳統的治癒模式並不直接定義哪個狀態為免疫，因此真正的免疫者是與非免疫的設限資料混合在一起。當模式假設不夠強時﹝指在無母數的情形下﹞，則強烈的需要依賴資料的“良好品質＂(指追蹤時間充份)，才能做正確的推論。第三章討論的治癒模式，分析的結果和人為決定的切點有關，所以應用有限。在本章中我們討論第三類型的治癒模式，免疫與否的定義取決於競爭風險發生的種類或是次序。我們先討論一個簡化的例子，是將 Betensky 與 Schoenfeld。2001﹞討論的新生兒因急性肺炎住院改編成因 SARS 住院的例子。令 1

T = time to hospital discharge (因 SARS 入院到活著出院的時間)

2

T = time to death (因 SARS 入院到在醫院死亡的時間)。

在此例中，“活著出院＂可視為痊癒。也就是說“免疫與否＂取決於競爭風險 (“死亡＂與“出院＂)發生的次序，因此若是在沒有外來設限的情形下，“痊癒＂是可被觀察到的事件。因為免疫與否有清楚的定義，這樣的架構不致發生前述不可辨識性的問題。因 SARS 入院死在醫院活著出院圖 4-1：SARS 入院的例子

(26)

以競爭風險的角度來看上述的問題，在這個簡化的例子中，只有兩種失敗型態(failure types)。令B~ =1，代表“活著出院＂的情形；B~ =2，代表“在醫院中死亡＂的情形。若令 “ 活著出院＂代表痊癒，則無法痊癒的比例為 ) 2 ~ Pr(B = 。然而對B~ =2的觀測值，的定義不明，有的做法是對死在醫院者時，令；或是只定義 1 T ∞ = 1 T T₁|B~=1。可得 ) 2 ~ Pr( ) 1 ~ Pr( ) 1 ~ | Pr( ) Pr(T₁ >t = T₁ >t B= B= + B= 。可以清楚看到，代表不會發生活著出院的比例。值得一提的是， ) 2 ~ Pr( ) Pr( lim ₁ > = = ∞ → T t B t ) Pr( ) Pr( 1− T1 >t = T1 ≤t 是競爭風險文獻裏經常討論的“累積發生函數＂( ，代表到時間所累積觀察到發生活著出院的比率。同理可得

cumulative incidence function)

t ) 1 ~ Pr( ) 2 ~ Pr( ) 2 ~ | Pr( ) Pr(T₂ >t = T₂ >t B= B= + B= ，其尾端機率Pr(B~=1)代表不會死於醫院的比例。

文獻中考慮競爭風險的治癒模式有 Greenhouse and Wolfe (1984)， Larsen and Dinse (1985)，Taylor (1995)，Ng and McLanchlan (1998)，Betensk and Schoenfeld (2001)，Maller and Zhou (2002)，Wang (2004)…等。由發表年代看來，這個研究方向似乎有變得熱門的趨勢。我們將統整這些文章的符號，以便做有系統的整理。這一系列文獻中，經常引用的文獻始自 Greenhouse and Wolfe (1984)，多數的文章討論不只兩種競爭風險。在此，我們由 Larsen & Dinse (1985)開始介紹，因為此文提供的思路為後續文章主要的脈絡。

4.1 Larson-Dinse 的迴歸分析

(27)

j p j B~= )= Pr( ， p₁+_L+p_J =1。再定義 T_j |B~= j為給定第個事件會發生的時間長度，則條件存活函數為 j ) ~ | Pr( ) (t T t B j Q_j = _j > = 。 Larson-Dinse﹝1985﹞以迴歸模式來描述與。令代表解釋變數，在失敗型態的部分用 logistic regression 做為模式的假設，則每一種失敗型態發生的機率可以表示成 j p Q_j(t) Z: p×1

∑

= = = _J j T j T j Z Z Z j B 1 ) exp( ) exp( ) | ~ Pr( β β ，其中β_j =(β_j₁,...,β_jp)；發生第個事件的時間長度，j T_j|Z,B~= j，則假設服從 Cox 模式，則條件存活函數可表示為

∫

− = = > = t T j j j t Z T t B j h u Z du Q 0 ) exp( ) ( 1 exp{ ) ~ | Pr( ) | ( γ ，其中風險函數的型態給定為 piece-wise exponential 模式，可將時間資料分成 ) (t h_j M 個區段(I₁,_L,I_M)，第個風險型態的風險函數為 j ) , , 1 , ( ) exp( ) (t t I m M h_j = α_jm ∈ _m = _L ，因此在每個區段的為常數。在風險函數給定的情況下，可用傳統的最大概似法以估計參數。然而若是基底風險函數不給定，估計將因未知參數的維度太高而變得複雜。然而根據第一類治癒模式所發展的方法概念，應可以類推到處理競爭風險的情形。 ) (t h_j h_j(t) 我們在此回顧基本的推論架構。在沒有設限的情形下，觀察到的資料型態為，，，則概似函數可以表示為： ) , , ~ (B_i x_i z_i B~ =_i 1,...,J i=1,...,n

∏∏

= = = = = = n i j B I i J j i i j j i z j B z x f J j h L 1 ) ~ ( 1 )} | ~ Pr( ) | ( ) ,..., 1 , , , (β γ ，

(28)

其中 t z t Q z t fj j ∂ ∂ − = ( | ) ) | ( 。而在有設限的情況下，有的觀測值無法得知真正的失敗型態，此時所觀測到的時間x_i =C_i，且。為了傳遞重要的概念，我們假設較簡單的情形，令 i Ji i T x T ,..., )> min( ₁ 2 = J ，所觀察到的資料型態為: ，，其中 ) , , (b_i x_i z_i 2 , 1 , 0 = i b (i=1,...,n) bi =1⇒B~ =i 1，xi = ，代表確知痊癒(or immune); T1i ，，代表確知未痊癒(or susceptible)；然而當時，的值未知，即不知道免疫與否，此時 2 = i b ⇒B~ =i 2 xi =T2i bi =0 i B~ xi =Ci <min(T1i,T2i)。設限下的概似函數可以表示為 ) 2 ( 2 2 ) 1 ( 1 1 1( | ) ( )] [ ( | ) ( )] [ = = =

∏

= i I bi i i i b I i i i n i F f x z p z f x z p z L 。 ) 0 ( 2 2 1 1( ) ( | ) ( ) ( | )} { + = × Ibi i i i i i i Q x z p z Q x z z p 和第二章的問題類似，最後一項在取對數後再對參數微分的計算會變得太過複雜，因此可用 EM 演算法來估計。假設Bi可觀察到完整的資料之概似函數為 ~ ) 2 ( 2 2 ) 1 ( 1 1 1( | ) ( )] [ ( | ) ( )] [ = = =

∏

= i I bi i i i b I i i i n i C f x z p z f x z p z L ) 2 , 0 ( 2 2 ) 1 , 0 ( 1 1( ) ( | )} { ( ) ( | )} { = = = = × i i Ibi Bi i i i B b I i i i Q x z p z Q x z z p 。取對數之後為

∑

= + = = n i i i i i C I b p z f x z L 1 1 1( ) log ( | )] )[log 1 ( { log )] | ( log ) ( )[log 2 (b_i p₂ z_i f₂ t_i z_i I = + + )] | ( log ) ( )[log 1 ~ , 0 (bi Bi p1 zi Q1 xi zi I = = + + 。 )]} | ( log ) ( )[log 2 ~ , 0 (bi Bi p2 zi Q2 xi zi I = = + + EM 演算法分為兩個部分，第一部分為先對取條件期望值，稱為 E-step，也就是求，其中牽涉到計算 C L log ] , , | [logL_C b_i x_i z_i E ) , | ~ Pr( ) 0 ( ) ( ] , , | ) ~ ( [I Bi j bi xi zi I bi j I bi Bi j X xi zi E = = = + = = =

(29)

) | ( ) 0 ( ) (b_i j I b_i w_j x_i z_i I = + = = ，在這裡w_j(x_i |z_i)可表示為 ) | ( ) ( ) | ( ) ( ) | ( ) ( ) | ( 2 2 1 1 i i i i i i i i j i j i i j z x Q z p z x Q z p z x Q z p z x w + = 。第二部分為對對數概似函數的期望值取極大值，稱為 M-step ，也就是求。經過疊代，就可以得到所求的估計值。值得一提的是牽涉到未知參數的部分為代入前階段之估計值，在求極大值的過程中視為定值。以上是以 2 種失敗型態來作為 Larsen-Dinse 在考慮競爭風險的治癒模式下的觀念闡述，至於種失敗型態也可以利用此觀念加以推廣。 ] [log maxE L_C w_j(x_i |z_i) J

文章分析了著名的心臟移植資料(Stanford Heart Transplant)，希望對於接受心臟移植後的病人(共 65 人)能探討解釋變數對於不同失敗型態的影響。共有兩個失敗型態，其中B~ =1代表因為發生排斥而死亡，B~ =2代表因為其他原因而死亡。在實際資料中， (29 人)， (12 人)，值得注意的是 (24 人)，代表缺失值頗嚴重。解釋變數的維度，除了截距項以外還包含“mismatch score＂ (越大代表捐贈者和接受者組織之相容度差)，“age＂代表接受移植時的年紀，“waiting time＂代表等到心臟移植的時間。其中“mismatch score＂和“age＂變數都經過標準化， time＂則是轉成是否超過 31 天的 0/1 變數。對於的模式採 piecewise-exponential 分配，試了三種配置法:當 45 . 0 ) 1 r( Pˆ b= = Pˆr(b=2)=0.18 37 . 0 ) 0 r( Pˆ b= = p=4 “waiting ) | (t Z Q_j 3 = M 時，將時間分成 3 個區段[0,45),[45,90),[90,∞);當M =2時，將時間分成 2 個區段 ;另一個選擇是 ) , 60 [ ), 60 , 0 [ ∞ 1 = M 。當發現“waiting time＂沒有太大影響時這個解釋變數便被捨棄，此時作者只考慮將“mismatch score＂和“age＂並以M =3的 piecewise-exponential

(30)

影響，倒是對於反映在上的排斥時間有顯著影響，年紀越輕且組織 match 較好的人到排斥發生的時間越長。換言之年紀大小和是否因排斥而死亡無關，但若是終究因排斥而死，年輕者之排斥現象會較年長者更晚發生﹝年輕者手術後可以撐比較久﹞。 ) ( 1 t Q 4.2 Logistic/Weibull 模式 -- Taylor

Larson and Dinse﹝1985﹞的分析，對於每一種失敗型態發生的機率( )

和給定第個事件會發生的條件存活函數( )皆做了母數模式的假設。而在 Taylor﹝1995﹞一文中，對於每一種失敗型態發生的機率( )仍做 logistic regression 的假設，以為例令 j p j Q_j(t) j p 2 = J ) exp( 1 ) exp( ) | 1 ~ Pr( _' ' 1 Z Z Z B p β β + = = = ，然而對於 Q_j(t|Z) (j=1,2) 這個部分，對於可能發病者 ( B~ =2 ) 假設並以無母數 product-limit 的方式分析之;但是對於免疫的人 ( ) ( ) | ( 2 2 t Z Q t Q = 1 ~ = B )視T =∞，所以避免了對做分配的假設。在前述假設下，可以得到概似函數 ) | (t Z Q1 ) 2 ( 2 2 ) 1 ( 1( )] [ ( | ) ( )] [ = = = i Ibi i i i b I i F p z f x z p z L ) 0 ( 2 2 1( ) ( | ) ( )] [ + = × Ibi i i i i Q x z p z z p ，

在這裡。和 Larson and Dinse (1985) 所提出的估計一樣，

最後一項在取對數後再對參數微分的計算會變得太過複雜，所以需要先做處理簡化，所得到的概似函數為 ) | ( ) ( ) | ( 2 2 t z h t Q t z f = F L ) 2 ( 2 2 ) 1 ( 1( )] [ ( ) ( | ) ( )] [ = = = i Ibi i i i i b I i C p z h x Q x z p z L ) 2 ~ , 0 ( 2 2 ) 1 ~ , 0 ( 1( )] [ ( | ) ( )] [ = = = = × i i I bi Bi i i i B b I i Q x z p z z p ，

(31)

取對數之後為

∑

= = = n i i i C I b p z L 1 1( )] )[log 1 ( { log )] ( log ) | ( log ) ( )[log 2 (b_i p₂ z_i Q₂ x_i z_i h x_i I = + + + )] ( )[log 1 ~ , 0 (bi Bi p1 zi I = = + 。 )]} | ( log ) ( )[log 2 ~ , 0 (bi Bi p2 xi Q2 xi zi I = = + + Taylor 也是建議以 EM 演算法來做估計，所需要的工作亦是針對簡化後概似函數中的和求條件期望值。對於的條件期望值 ) 1 ~ , 0 (b_i = B_i = I I(b_i =0,B~_i =2) I(b_i =0,B~_i =2) ) , | 2 ~ Pr( ) 0 ( ) 2 ( ] , , | ) 2 ~ ( [I Bi bi xi zi I bi I bi Bi Xi xi zi E = = = + = = = =I(b_i =2)+I(b_i =0)w₂(x_i |z_i)， Taylor 導出 ) | ( ) ( ) ( ) | ( ) ( ) , | 2 ~ Pr( ) | ( 2 2 1 2 2 2 i i i i i i i i i i i i i z x Q z p z p z x Q z p z x X B z x w + = = = = ，， ) | ( 1 ) | ( ₂ 1 xi zi w xi zi w = − 此處利用 Kaplan-Meier 所表示成的 product-limit 形式代入，透過 EM 演算法即得所求之估計值。 ) | ( 2 xi zi Q 利用以上 logistic / Kaplan-Meier 做為模式假設會發生的問題是過多的未知參數造成估計的繁雜，此外理論上lim_t_→_∞Q₁(t|Z)=0，但是在模擬分析上卻依然可能產生致病時間模式不會遞降到 0 的情形，尤其是當樣本數太小或是設限資料的比例太大的時候。針對這個問題 Taylor 亦提出修補的方法可強迫的無母數估計值遞降至 0。 ) | ( 1 t Z Q 4.3 母數分析 – 非迴歸模式:

Greenhouse and Wolfe(1984)提出以最大概似法做為推論方法，並針對常見模式做細節討論。Ng and McLanchlan (1998)亦是採取母數模式的架構，因為

(32)

提出的推論方法較為創新，因此回顧這篇文章。該文主要目的是估計T₁|B~=1 的分配與 p ，至於T₂ |B~=2的分佈較不重要。值得一提的是當 B~= j ，令。在 mixture 的表示法下： ∞ = k T (k ≠ j) ) , Pr( ) (t T1 t T2 t H = > > ) 2 ~ , Pr( ) 1 ~ , Pr( ₁ > = + ₂ > = = T t B T t B ) 1 )( ( ) ( ₂ 1 t p Q t p Q + − = 。論文中討論 likelihood-based 的推論方法。如果和的母數分配均為已知，則 log-likelihood 可以表示為 ) ( 1 t Q Q2(t) )} ; ( log{ ) 1 ( ) ( log ₁ 1 1 θ ψ n _j j j pf x b I L

∑

= = = ( 2)log{(1 ) ( ; ₂)} 1 2 j θ n j j p f x b I

∑

= − = + )} ; ( log{ ) 0 ( 1 ψ j n j j H x b I

∑

= = + ，其中 t t Q t f_j j ∂ ∂ − = ( ) ) ( ，ψ =(θ₁,θ₂,p)。假設我們只關心(θ₁,p)，視θ 為不感興₂ 趣的參數。以此為前提下作者認為 full-likelihood 需要對兩個分佈均做母數分配的假設並沒有必要，並進一步提出了 partial Maximum likelihood 的估計

方法。他們的想法是只對的函數型態做母數的假設，對於的資料給予密度函數，即 ) ( 1 t Q b_i =1 ) ; ( ₁ 1 tj θ pf ，但是對於bi =0或是bi =2的觀察值，都視為設限並假 設觀察到研究的結束(即 C )，此時在概似函數中給予 1_i ( −p)+ pQ₁(C_i;θ₁)。所以修正後的部份概似函數可以表示為 )} ; ( log{ ) 1 ( ) ( ~ log ₁ 1 1 θ ψ n _i i i pf x b I L

∑

= = = )} ; ( ) 1 log{( ) 2 ( * ₁ 1 1 i θ n i i p pQ x b I

∑

= + − = + )} ; ( ) 1 log{( ) 0 ( ₁ 1 1 i θ n i i p pQ x b I

∑

= + − = + ，其中最值得注意的是對b_i =2的觀測值，不能使用T₂ =x_i而是更大的。因為若是將以取代(即競爭風險視為外生的設限情形)，所得之估計量會有偏 * i x C= * i x x_i

(33)

誤。

Maller and Zhou (2002) 考慮更一般化的競爭風險架構，可以容許更多的失敗

型態，所以 Pr(~ ) 1。再定義以下存活函數 1 1 = = =

∑

= = J j j J j p j B ) ~ | Pr( ) (t T t B j Q_j = _j > = ，

∑

= = − − = = ≤ − = > > = J j J j j j j J t T t B j Q t p T t T t H 1 1 1 ) 1 {1 ( )} ~ , Pr( 1 ) ,..., Pr( ) ( 。令(T_ij,C_i) (i=1,...,n)為(T_j,C)的隨機樣本(j=1,...,J)，代表實際發生事件類型的指標函數定義為δ_ji =I(B~_i = j,min(T₁_i,...,T_Ji)<C_i)，

∑

，其中 = = J j ji i 1 δ δ 所以若是第i個人設限，δ_i =0，否則δ_i =1，此外令X_i =min(T₁_i,...,T_Ji,C_i)，資料可表示為{(δ₁_i,...,δ_Ji,X_i,δ_i),(i=1,...,n)}。概似函數就可表示為 i ji i J j j j n i J j i j jf x p Q x p Lψ δ −δ = = =

∑

∏ ∏

− − = 1 1 1 1 )]} ( 1 [ 1 }{ ] ) ( [ { ) ( ，其中−∂Q_j(t)/∂t = f_j(t)。以上討論的兩篇以母數分析為出發點的文章，如何在“最大概似法＂架構下，提出值得發表的創意。Ng and McLanchlan (1998) 提出較具穩健性的 partial ML 方法；Maller and Zhou (2002) 則是理論推導，做為檢定

的假說的基礎。在給定下，而且外生的設限存在，我們往往不知道 1 1 =

∑

= J j j p * J J = } 0 : {i δ_i = 的設限觀測值中究竟是否包含未出現過的失敗型態(指，此時 ) 。這個問題的方向在於希望得知參數的值是否發生在邊界 (boundary)，是非典型的推論問題，許多數學技巧都不能使用。 * J J > 1 1 <

∑

= J j j p 4.4 無母數分析 – Wang

Wang (2004) 只考慮 2 種競爭風險，並以 multi-state model 描述問題， 文中令 B~做為路徑的指標，並提出以無母數角度估計p,Q₁(t),Q₂(t)的方法。為求

(34)

符號的一致性，資料型態可以表示為{(bi,xi)(i=1,...,n)}。可知當b=0，B ~ 的值未知，在給定資料後B~ =1的條件機率為 ) , 0 | 1 ~ Pr( ) (c B b X c p = = = = ) Pr( ) Pr( ) Pr( ) 1 ~ , Pr( 2 1 2 1 c C c T T c C B c T T = > ∧ = = > ∧ =

∫

> = + ∈ = c v B dv v v T c H 1) ~ ), , [ Pr( ) ( 1 1 ，這裡H(t)=Pr(T₁ >t,T₂ >t)。在假設T₁、T₂和C為連續型隨機變數，則可得到 ) 1 ~ , ~ Pr( ) ) Pr( ) 1 ), , [ Pr( ) ( 1 ) ( + < 1 = > = + ∈ =

∫

> B T T v C b dv v v X c H c P c v ，其中T~=sup{t:Pr(T₁ >t,T₂ >t)Pr(C >t)>0} ) 1 ，代表可觀察時間的上界。值得一提的是Pr(T~<T₁,B~= 永遠無法被準確估計，除非是Pr(T~<T₁,B~=1)=0，這個條 Kaplan-Meier 估計量估計和如下: 件就是追蹤時間為充份。可利用 H(t) ) Pr( ) (t C t G = >

∏

∑

≤ = = ≥ ≠ = − = t u n i i n i i i u X I b u X I t H } ) ( ) 0 , ( 1 { ) ( ˆ 1 1 ，

∏

∑

≤ = = ≥ = = − = t u n i i n i i i u X I b u X I t G } ) ( ) 0 , ( 1 { ) ( ˆ 1 1 。對於Pr(X = bv, =1)，因為是可觀察變數之機率函數則可用 empirical function 來做估計。而Pr(T~<T₁,B~=1)可寫成p(T~)H(T~)，對於 p(T~)和H(T~)可以利用觀察到最大b≠0的值x_(n₎來做估計。因此可以得到p(c)的估計量為 ) ~ ( ˆ ) ~ ( ~ ) ( ˆ ) 1 , ( ) ( ˆ 1 ) ( ˆ ₍ ₎ ₍ ₎ 1 n n c v n i i i X H X p v G n b v X I c H c p =

∫ ∑

= = + > = ，其中~p(X~_(n₎)的值需要額外假設才能得知。

治癒模式之文獻回顧

國立交通大學

統計學研究所

碩 士 論 文