帶終止事件的復發(fā)事件數(shù)據(jù)加速均值模型

關(guān)楷諭，戴家佳

帶終止事件的復發(fā)事件數(shù)據(jù)加速均值模型

關(guān)楷諭，戴家佳*

(貴州大學理學院，貴州貴陽550025)

對帶終止事件的復發(fā)事件數(shù)據(jù)提出加速均值模型，運用廣義估計方程思想和逆概率加權(quán)方法給出參數(shù)估計，并證明所得估計的相合性和漸近正態(tài)性.

復發(fā)事件;終止事件;廣義估計方程;逆概率加權(quán);加速均值模型

生存分析廣泛應(yīng)用于醫(yī)療、保險、生物學等領(lǐng)域中，主要研究生存時間或比較不同個體之間的差異.如果對一個個體，某些感興趣的特定事件在一段時間內(nèi)重復發(fā)生，稱為復發(fā)事件，其相應(yīng)的數(shù)據(jù)稱為復發(fā)事件數(shù)據(jù).研究的個體中途退出研究過程或試驗結(jié)束稱為刪失事件，其相應(yīng)的時間稱作刪失時間.在研究復發(fā)事件數(shù)據(jù)時，假設(shè)當協(xié)變量給定時，刪失事件與復發(fā)事件獨立.文獻［1－2］考慮了基于歷史數(shù)據(jù)的條件強度模型.文獻［3－4］不考慮歷史數(shù)據(jù)，建立了比率和均值模型，相比于強度模型放寬了條件使之更加穩(wěn)健.文獻［5－6］提出了脆弱變量以描述復發(fā)事件間的關(guān)系.早期研究中，將死亡視為刪失事件.但死亡和疾病的復發(fā)可能有很強的相關(guān)性，比如腫瘤患者病情不斷復發(fā)，死亡概率也會變高，他們之間很有可能不獨立，于是提出終止事件概念，即個體死亡.對于帶終止事件的復發(fā)事件數(shù)據(jù)，一般采用邊際模型方法和隨機效應(yīng)模型方法.隨機效應(yīng)模型借助潛變量刻畫復發(fā)事件和終止事件的相關(guān)性，并假設(shè)在給定脆弱變量時復發(fā)事件與終止事件獨立，如文獻［7－11］;邊際模型方法則側(cè)重于復發(fā)事件與終止事件的邊際模型，不考慮之間的相關(guān)性，如文獻［12－16］，但還未有通過邊際模型方法對帶終止事件的加速均值模型的分析.加速均值模型是對基本均值函數(shù)中的時間做了尺度變換，具有結(jié)構(gòu)簡單、解釋性強的特點，同樣具有研究意義.

在本文中，將利用廣義方程估計思想和生存逆概率加權(quán)(IPSW)方法研究帶終止事件的復發(fā)事件數(shù)據(jù)加速均值模型.

1 模型與參數(shù)估計

設(shè)N*(t)為在［0，t］時間內(nèi)復發(fā)事件發(fā)生的次數(shù)，C表示刪失時間，D為終止事件發(fā)生的時間，Z表示協(xié)變量.終止事件的發(fā)生會導致復發(fā)事件過程停止，也就是說一般情況下，只能觀測到D和C的最小值，記X=D∧C，δ=I(D≤C)，其中a∧b=min (a，b)，I(·)是示性函數(shù).可觀測的復發(fā)事件次數(shù)為N(t)=N*(t∧X).假設(shè)在給定協(xié)變量條件下，C和(D，N*(·))是獨立的.可觀測的數(shù)據(jù)為{Ni(t)，Xi，δi，Zi，0≤t≤Xi，i=1，2，…，n}.

加速均值模型具有如下的形式其中，β0是未知回歸參數(shù)，μ0(t)是未知基本均值函數(shù).記珟Ni(t;β)=Ni(te－β'Zi)，Yi(t;β)=I(Ci≥te－β'Zi).當不考慮終止事件(D=∞)時，定義如下過程

如果模型(1)成立，則M0i(t;β0)是一個零均值過程.由估計方程思想［17］得

其中，τ是一個事先給定的常數(shù)，使得P(Ci≥τe－β0'Zi)＞0.Q(t;β)是一個給定的權(quán)函數(shù)，

當考慮終止事件時，部分Ci可能不能被觀測到，從而導致Yi(t;β)無法被觀測到.逆概率加權(quán)方法［18］是處理帶缺失的數(shù)據(jù)的常用方法之一:對完全情形下的估計方程的貢獻項進行加權(quán)，并當權(quán)取為選擇概率的逆時，定義的估計在通常情況下是相合的.對于模型(1)，可采用2種方法［14］，一是IPCW方法，對刪失數(shù)據(jù)進行建模.這種想法很直觀，建模完成后能較直接地給出缺失值的估計，但它的缺點是在多數(shù)情況下對發(fā)生刪失的情況并不關(guān)心，用過多的精力對刪失時間建模意義不大.所以考慮IPSW方法，即:通過對終止事件的生存函數(shù)進行建模，從而對Yi(t;β)給出合適的估計，同時生存函數(shù)也是感興趣的.

定義ωi(t;β)=I(Xi≥te－β'Zi)/S(te－β'Zi|Zi)，其中S(t|Zi)=P(D＞t|Zi).易知E(ωi(t;β))= E(Yi(t;β)).假設(shè)終止事件發(fā)生的時間Di滿足比例風險模型:λD(t|Zi)=(t)，其中γ0是未知回歸參數(shù)，(t)是未知基本風險函數(shù).記

其中

令

其中{G1，…，Gn}是相互獨立的服從標準正態(tài)分布的隨機變量.

固定觀測到的數(shù)據(jù)集{Ni(·)，Zi，Xi，δi，i=1，…，n}，重復產(chǎn)生{G1，…，Gn}，記是(3)的解.由文獻［26］可知具有相同的極限分布.為得到^β的方差估計，可用的經(jīng)驗協(xié)方差矩陣作為^β的協(xié)方差矩陣估計.也給出了(t)是μ0(t)相合估計和的漸近正態(tài)性的證明.

2 估計的漸近性質(zhì)

為了研究未知參數(shù)的漸近性質(zhì)，需要假設(shè)如下條件:

(C1){Ni(·)，Zi，Xi，δi}，i=1，…，n是獨立同分布的;

(C3)Ni(t)和Zi在［0，τ］有界，權(quán)函數(shù)Q(t; β0)有有界變差且在t∈［0，τ］一致依概率收斂到一個確定的函數(shù)q(t);

(C4)Cieβ0'Zi的密度函數(shù)有界，μ0(t)具有有界的二階導數(shù);

(C5)記u(β)表示n－1U(β)的極限，存在一個β0的緊鄰域N(β0)，滿足當β∈N(β0)，β≠β0時，有u(β)≠u(β0);

(C6)矩陣

定義如下過程:

所以可得Mi(t;β0)是一個零均值過程.

記

相關(guān)參數(shù)估計的漸近性質(zhì)由以下定理給出，證明參見附錄.

定理1在正則條件(C1)～(C3)下U(β0)是漸近正態(tài)的，其均值為0，協(xié)方差矩陣的相合估計

定理2在(C1)～(C6)的正則條件下是強相合的(β^－β0)是漸近正態(tài)的，其均值為0，協(xié)方差矩陣的相合估計是－1－1的正態(tài)分布，其中

定理3在(C1)～(C6)的正則條件下(t)在t∈［0，τ］上是μ0(t)的是強相合估計，［(t)－μ0(t)］是漸近正態(tài)的，其均值為0，在(t1，t2)∈(［0，τ］，［0，τ］)處的協(xié)方差函數(shù)的相合估計是

3 定理的證明

定理1的證明

由函數(shù)delta方法和鞅的中心極限定理可得:

其中

代入上式并交換積分次序可得

其中

其中

定理2的證明給定任意dn→0，對于‖ββ0‖≤dn有

由(5)式易得

再由文獻［27］的定理1的技術(shù)可知，(6)式的第1和第2部分都是o().第3部分:

由泰勒展開可得

代入上式得

綜上得到

由一致強大數(shù)定律可得n－1U(β)在N(β0)內(nèi)一致收斂到u(β)，易得u(β0)=0.結(jié)合(C5)可得^β是β的強相合估計.再由^β的定義與A可逆得

漸近收斂到均值為0，協(xié)方差陣為A－1ΣA－1的正態(tài)分布.

定理3的證明給定任意dn→0，對于‖ββ0‖≤dn有

由(5)式易得

結(jié)合文獻［27］的定理1的技術(shù)可知，(7)式右邊第一個部分是o().

其中

記PD(u，t)和ΓD(t)分別是D(u，t)和D(t)的極限.由^β的相合性，類似定理2證明可得^μ0(t)相合性和漸近正態(tài)性

其中

4 結(jié)論

在本文中對帶終止事件的復發(fā)事件數(shù)據(jù)提出了加速均值模型，通過生存逆概率加權(quán)(IPSW)方法對不可觀測值做出估計并帶入完全情形下的估計方程，得到未知參數(shù)的相合估計和其漸近性質(zhì).這個過程中包含對生存函數(shù)假設(shè)建模.首先這是感興趣的，另一方面，生存分析中對失效時間建模的理論較完善，可供參考的文獻也較多.在加速均值模型參數(shù)估計中，為簡化方差的估計，使用重抽樣方法.在估計方程中涉及權(quán)函數(shù)Q(t)，增加權(quán)函數(shù)的目的是使估計方程凸性增大利于求解及使估計方差減小但尋找困難，常用的權(quán)函數(shù)為log－rank權(quán)函數(shù)或Gehan型權(quán)函數(shù).

致謝貴州大學引進人才科研項目(2009－070)對本文給予了資助，謹致謝意.

［1］ANDERSEN P K，GILL R D.Cox’s regression model for counting processes:a large sample study［J］.Ann Statis，1982:1100－1120.

［2］ZENG D，LIN D Y.Efficient estimation of semiparametric transformation models for counting processes［J］.Biometrika，2006，93(3):627－640.

［3］PEPE M S，CAI J.Some graphical displays and marginal regression analyses for recurrent failure times and time dependent covariates［J］.J Am Statistical Association，1993，88(423):811－820.

［4］LIN D Y，WEI L J，YANG I，et al.Semiparametric regression for the mean and rate functions of recurrent events［J］.J Royal Statistical Society:Statistical Methodology，2000，B62(4):711－730.

［5］NIELSEN G G，GILL R D，ANDERSEN P K，et al.A counting process approach to maximum likelihood estimation in frailty models［J］.Scandinavian J Statistics，1992:25－43.

［6］VEKEMANS D，PROOST S，VANNESTE K，et al.Gamma paleohexaploidy in the stem－lineage of core eudicots:significance for mads－box gene and species diversification［J］.Molecular Biology and Evolution，2012，29(12):3793－3806.

［7］HUANG C Y，WANG M C.Joint modeling and estimation for recurrent event processes and failure time data［J］.J American Statistical Association，2004，99(468):1153－1165.

［8］LIU L，WOLFE R A，HUANG X.Shared frailty models for recurrent events and a terminal event［J］.Biometrics，2004，60(3): 747－756.

［9］YE Y，KALBFLEISCH J D，SCHAUBEL D E.Semiparametric analysis of correlated recurrent and terminal events［J］.Biometrics，2007，63(1):78－87.

［10］HUANG C Y，QIN J，WANG M C.Semiparametric analysis for recurrent event data with time－dependent covariates and informative censoring［J］.Biometrics，2010，66(1):39－49.

［11］戴家佳，關(guān)楷諭，吳歡.帶有終止事件的復發(fā)事件數(shù)據(jù)的加性加速比率模型［J］.應(yīng)用數(shù)學學報，2015，38(4):735－750.

［12］COOK R J，LAWLESS J F.Marginal analysis of recurrent events and a terminating event［J］.Statistics Medicine，1997，16(8): 911－924.

［13］GHOSH D，LIN D Y.Nonparametric analysis of recurrent events and death［J］.Biometrics，2000，56(2):554－562.

［14］GHOSH D，LIN D Y.Marginal regression models for recurrent and terminal events［J］.Statistica Sinica，2002，12(3):663－688.

［15］ZHAO H，ZHOU J，SUN L.A marginal additive rates model for recurrent event data with a terminal event［J］.Commun Statis: TM，2013，42(14):2567－2583.

［16］何穗，程希明，周潔.帶終止事件的多類型復發(fā)事件的一般加性乘積比例模型［J］.應(yīng)用數(shù)學學報，2012，35(5):804－816.

［17］LIANG K Y，ZEGER S L.Longitudinal data analysis using generalized linear models［J］.Biometrika，1986，73(1):13－22.

［18］王啟華，史寧中，耿直.現(xiàn)代統(tǒng)計研究基礎(chǔ)［M］.北京:科學出版社，2010.

［19］COX D R.Partial likelihood［J］.Biometrika，1975，62(2):269－276.

［20］BRESLOW N.Covariance analysis of censored survival data［J］.Biometrics，1974:89－99.

［21］GHOSH D.Accelerated rates regression models for recurrent failure time data［J］.Lifetime Data Analysis，2004，10(3): 247－261.

［22］LIN D Y，WEI L J，YING Z.Accelerated failure time models for counting processes［J］.Biometrika，1998，85(3):605－618.

［23］JIN Z，LIN D Y，WEI L J，et al.Rank－based inference for the accelerated failure time model［J］.Biometrika，2003，90(2): 341－353.

［24］JIN Z，LIN D Y，YING Z.On least－squares regression with censored data［J］.Biometrika，2006，93(1):147－161.

［25］LIN D Y，GEYER C J.Computational methods for semiparametric linear regression with censored data［J］.J Comput Graph Statis，1992，1(1):77－90.

［26］PARZEN M I，WEI L J，YING Z.A resampling method based on pivotal estimating functions［J］.Biometrika，1994，81(2): 341－350.

［27］YING Z.A large sample study of rank estimation for censored regression data［J］.Ann Statis，1993:76－99.

An Accelerated Mean Model for Recurrent Events Data with Informative Terminal Event

GUAN Kaiyu，DAI Jiajia
(College of Science，Guizhou University，Guiyang 550025，Guizhou)

In this paper，we propose an accelerated mean model for recurrent events data with informative terminal event.Based on generalized estimating equation and inverse probability weighting technique，the consistency and asymptotic normality properties of the proposed estimators are proved.

recurrent events;terminal event;generalized estimating equation;inverse probability weighting;accelerated mean model

O212.7

A

1001－8395(2016)03－0362－07

10.3969/j.issn.1001－8395.2016.03.011

(編輯周俊)

2015－11－20

國家自然科學基金(11361015)和貴州省科學技術(shù)基金(2009－2063)

*通信作者簡介:戴家佳(1976—)，女，教授，主要從事生存分析的研究，E－mail:jjdai@gzu.edu.cn

2010 MSC:62G05;62N01