李云霞，周杏杏

（浙江財(cái)經(jīng)大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，浙江 杭州 310018）

含有協(xié)變量的復(fù)發(fā)事件變點(diǎn)模型的參數(shù)估計(jì)

李云霞，周杏杏

（浙江財(cái)經(jīng)大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，浙江 杭州 310018）

針對(duì)復(fù)發(fā)事件數(shù)據(jù)協(xié)變量的重要作用，建立含有協(xié)變量的復(fù)發(fā)事件變點(diǎn)模型，考慮協(xié)變量作用于強(qiáng)度率函數(shù)的情形。對(duì)于此模型，使用最大似然方法得到變點(diǎn)及各參數(shù)估計(jì)，并得到了變點(diǎn)估計(jì)的相合性。最后對(duì)于同時(shí)存在待估參數(shù)和待估變點(diǎn)的似然函數(shù)，采用最速上升法進(jìn)行了數(shù)據(jù)模擬。

變點(diǎn)；復(fù)發(fā)事件；協(xié)變量；最大似然估計(jì)；最速上升法

一、引 言

復(fù)發(fā)事件數(shù)據(jù)是指對(duì)個(gè)體進(jìn)行觀察，某種感興趣事件重復(fù)發(fā)生的時(shí)間點(diǎn)、時(shí)間間隔和累積次數(shù)等所組成的數(shù)據(jù)，這類數(shù)據(jù)廣泛地應(yīng)用于生物統(tǒng)計(jì)、臨床醫(yī)學(xué)、工業(yè)可靠性控制和保險(xiǎn)精算等研究領(lǐng)域中，如觀測某些病人某種疾病的復(fù)發(fā)時(shí)間及復(fù)發(fā)次數(shù)，某些機(jī)器故障的多次發(fā)生時(shí)間等。因?yàn)槭录?fù)發(fā)的時(shí)間是有次序的并具有相依性，同時(shí)由于刪失時(shí)間的存在，以及刪失時(shí)間可能與事件發(fā)生的累積次數(shù)具有相依性，使得對(duì)復(fù)發(fā)事件數(shù)據(jù)的分析、建模及統(tǒng)計(jì)推斷變得十分困難。但由于復(fù)發(fā)事件數(shù)據(jù)具有廣泛的應(yīng)用性，近20年來，學(xué)術(shù)界對(duì)復(fù)發(fā)事件數(shù)據(jù)的研究引起了廣泛的重視，并得到了快速的發(fā)展，其研究結(jié)果不僅具有重要的理論意義，而且具有廣泛的應(yīng)用前景。越來越多的學(xué)者對(duì)復(fù)發(fā)事件數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析，如Wang，Pena等對(duì)復(fù)發(fā)事件數(shù)據(jù)進(jìn)行了深入的研究［1－3］。

另一方面，若生存數(shù)據(jù)的風(fēng)險(xiǎn)函數(shù)由于某種原因在某個(gè)未知時(shí)刻發(fā)生了變化，即生存數(shù)據(jù)的變點(diǎn)問題。Matthews和Farewell建立了最基礎(chǔ)的生存數(shù)據(jù)風(fēng)險(xiǎn)函數(shù)分段常數(shù)單變點(diǎn)模型［4］。這篇關(guān)于生存數(shù)據(jù)變點(diǎn)的文章，引發(fā)了人們對(duì)這類問題的研究興趣，發(fā)現(xiàn)了研究變點(diǎn)問題對(duì)于生存分析的重大意義。很多學(xué)者對(duì)生存數(shù)據(jù)中的變點(diǎn)進(jìn)行了估計(jì)，Chang，Chen和Hsiung提出了首先采用Nelson－Aalen非參數(shù)估計(jì)來估計(jì)變點(diǎn)，然后再使用極大似然函數(shù)來估計(jì)參數(shù)，并考慮了生存數(shù)據(jù)的刪失性，結(jié)合非參數(shù)和參數(shù)方法來提高變點(diǎn)和參數(shù)估計(jì)的精度［5］。Dupuy將常數(shù)變點(diǎn)模型進(jìn)行了推廣，加入了協(xié)變量對(duì)生存數(shù)據(jù)的影響，同時(shí)考慮了數(shù)據(jù)的刪失，提出了變點(diǎn)和參數(shù)的極大似然估計(jì)，并進(jìn)行了假設(shè)檢驗(yàn)［6－7］。Zhao，Wu和Zhou首先考慮了含有持久生存數(shù)據(jù)的變點(diǎn)風(fēng)險(xiǎn)函數(shù)模型，采用Kaplan－Meier估計(jì)，結(jié)合Chang研究中Nelson－Aalen非參數(shù)估計(jì)和極大似然參數(shù)估計(jì)來給出變點(diǎn)和參數(shù)估計(jì)，并證明了估計(jì)的相合性［8］。Li研究了含有協(xié)變量以及永久生存數(shù)據(jù)的變點(diǎn)模型，給出了變點(diǎn)及參數(shù)估計(jì)，并得到了估計(jì)的相合性［9］。以上這些都是關(guān)于生存數(shù)據(jù)中一次觀測數(shù)據(jù)，即非復(fù)發(fā)事件數(shù)據(jù)變點(diǎn)問題的研究。

由于復(fù)發(fā)事件數(shù)據(jù)的存在更具廣泛性，學(xué)者們開始考慮對(duì)復(fù)發(fā)事件數(shù)據(jù)進(jìn)行變點(diǎn)建模，但相關(guān)的文獻(xiàn)較少。Frobish對(duì)含有協(xié)變量的復(fù)發(fā)事件變點(diǎn)模型進(jìn)行了研究，研究過程中認(rèn)為協(xié)變量作用于變點(diǎn)上，此時(shí)由于協(xié)變量的不同會(huì)導(dǎo)致變點(diǎn)不同［10］。因此，本文考慮了協(xié)變量作用于強(qiáng)度率上的情況，這與實(shí)際情況更相符，并通過變點(diǎn)估計(jì)值的標(biāo)準(zhǔn)差對(duì)協(xié)變量作用于強(qiáng)度率上的合理性進(jìn)行了說明。下面將對(duì)這種情況進(jìn)行建模。

下文的組織結(jié)構(gòu)如下：第二部分采用最大似然法估計(jì)變點(diǎn)及參數(shù)，并且對(duì)變點(diǎn)估計(jì)值的一致性進(jìn)行說明；第三部分采用最速上升法進(jìn)行模擬，對(duì)不同樣本量的估計(jì)效果進(jìn)行比較。

二、變點(diǎn)及相關(guān)參數(shù)的估計(jì)

三、模 擬

通過模擬來評(píng)價(jià)用最大似然法得到的變點(diǎn)估計(jì)值的效果。在模擬時(shí)，將聯(lián)合對(duì)數(shù)似然函數(shù) （8）式化簡為關(guān)于β和τ2的函數(shù)。記η＝ （τ2， β） ，有：

對(duì)于同時(shí)存在待估參數(shù)和待估變點(diǎn)的似然函數(shù)，需要用到多變量的極大化數(shù)值技術(shù)。常用的方法包括三種：一是最速上升法，這需要一階導(dǎo)數(shù)向量；二是推廣的Newton－Raphson方法，需要一階和二階導(dǎo)數(shù)；三是Marquardt′s方法，需要使用一個(gè)混合常數(shù)。由于強(qiáng)度率的二階導(dǎo)計(jì)算困難，因此采用第一種方法，即最速上升法。在模擬時(shí)，通過最速上升法得到未知變點(diǎn)和協(xié)變量系數(shù)的估計(jì)值，再由這兩個(gè)估計(jì)值得到強(qiáng)度率函數(shù)中參數(shù)的估計(jì)值。

在本次模擬中，假設(shè)研究開始時(shí)間為0，τ1＝0，τ2＝5，τu＝15，刪失時(shí)間是隨機(jī)的。為不失一般性，同時(shí)假設(shè)λ1＜λ2。當(dāng)λ1＝0．2時(shí)，λ2＝0．3，1．0；當(dāng)λ1＝0．6時(shí)，λ2＝1．0，由此可得三個(gè)不同模型。由于τ1＝0，所以令λ0＝0。同時(shí)令協(xié)變量X表示性別，取值為0或1，分別表示男性和女性，協(xié)變量系數(shù)取值為β＝5，假設(shè)變點(diǎn)前后取值不變。下面利用最速上升法來得到相應(yīng)待估參數(shù)和變點(diǎn)的估計(jì)值，并通過表格來對(duì)估計(jì)值進(jìn)行比較說明。

表1 模型（1）的參數(shù)估計(jì)值絕對(duì)偏差

表2 模型（2）的參數(shù)估計(jì)絕對(duì)偏差

表3 模型（3）的參數(shù)估計(jì)絕對(duì)偏差

分別對(duì)樣本容量為100，200，500的樣本進(jìn)行變點(diǎn)和參數(shù)的估計(jì)，分別給出了參數(shù)λ1，λ2，系數(shù)β和變點(diǎn)τ2的估計(jì)值的絕對(duì)偏差。m如前所述，表示樣本量，bias表示估計(jì)值與真實(shí)值之間的絕對(duì)偏差，此值越小說明估計(jì)值越接近真實(shí)值。通過表1至表3可以發(fā)現(xiàn)，不論參數(shù)取值多少，樣本容量越大，估計(jì)效果越好。當(dāng)樣本容量為500時(shí)，各估計(jì)值的偏差要明顯小于樣本容量為100時(shí)的偏差，但是與樣本容量為200時(shí)的偏差不大。這說明，最大似然的估計(jì)方法得到的估計(jì)值不需要很大的樣本量就能得到比較好的估計(jì)效果。當(dāng)λ1＝0．2，λ2＝1．0，樣本容量為500時(shí)變點(diǎn)的絕對(duì)偏差只有0．01，估計(jì)效果明顯優(yōu)于λ1＝0．2，λ2＝0．3及λ1＝0．6，λ2＝1．0這兩種情況，且此時(shí)λ1和λ2的估計(jì)值的絕對(duì)偏差可以忽略不計(jì)，因此λ1＝0．2，λ2＝1．0時(shí)的估計(jì)效果最好。當(dāng)λ1＝0．6，λ2＝1．0時(shí)各估計(jì)值的偏差要小于λ1＝0．2，λ2＝0．3時(shí)的偏差，這說明變點(diǎn)前后強(qiáng)度率函數(shù)相差越大，估計(jì)效果越好，此時(shí)變點(diǎn)也越明顯。同時(shí)可以發(fā)現(xiàn)，協(xié)變量系數(shù)β的估計(jì)效果沒有參數(shù)λ1，λ2好，在λ1＝0．2，λ2＝1．0時(shí)的＾β的偏差比較小，說明用極大似然法得到的估計(jì)值是合理的。綜上所述，含有協(xié)變量的復(fù)發(fā)事件變點(diǎn)模型用最大似然法進(jìn)行變點(diǎn)估計(jì)是可行的。

四、結(jié) 論

復(fù)發(fā)時(shí)間間隔服從的分布最常見的是指數(shù)分布，相對(duì)應(yīng)的條件強(qiáng)度函數(shù)為常數(shù)，因此強(qiáng)度函數(shù)的分段常數(shù)變點(diǎn)模型是最基礎(chǔ)和最具應(yīng)用背景的變點(diǎn)模型。而復(fù)發(fā)事件的發(fā)生往往與眾多影響因素相關(guān)聯(lián)，例如年齡、性別、身高、體重等因素，本文在假設(shè)強(qiáng)度率函數(shù)為分段常數(shù)形式的基礎(chǔ)上，考慮了含有協(xié)變量的復(fù)發(fā)事件變點(diǎn)模型，對(duì)其進(jìn)行建模和參數(shù)估計(jì)。應(yīng)用最大似然法進(jìn)行變點(diǎn)及參數(shù)估計(jì)，同時(shí)得到了變點(diǎn)估計(jì)的一致性。對(duì)于多變量的似然函數(shù)，避免較大的估計(jì)誤差，采用最速上升法來進(jìn)行模擬，發(fā)現(xiàn)最大似然法可以很好地對(duì)變點(diǎn)進(jìn)行估計(jì)。對(duì)于協(xié)變量作用于強(qiáng)度率函數(shù)的機(jī)理上的原因如下。協(xié)變量可以是定性的，如性別、種族等，不隨時(shí)間變化；也可以是定量的，如年齡，身高等，隨時(shí)間變化。在生存分析中，更多的是協(xié)變量隨時(shí)間而變化，這類協(xié)變量稱為時(shí)間相依協(xié)變量。對(duì)于定性的協(xié)變量，由于其不隨時(shí)間變化，不論作用于變點(diǎn)上還是復(fù)發(fā)事件的強(qiáng)度率函數(shù)上，都表示了協(xié)變量對(duì)強(qiáng)度率函數(shù)造成影響。對(duì)于時(shí)間相依的協(xié)變量，如果只作用與變點(diǎn)上，就會(huì)忽略變點(diǎn)前后協(xié)變量對(duì)強(qiáng)度率函數(shù)的影響 ，只有作用與強(qiáng)度率函數(shù)上，才能充分利用協(xié)變量所包含的信息。并且在生存分析中，對(duì)于協(xié)變量一般采用Cox模型來處理，這類模型將協(xié)變量作用與風(fēng)險(xiǎn)函數(shù)上。由此可知，協(xié)變量作用與復(fù)發(fā)事件的強(qiáng)度率函數(shù)上更能充分利用已知信息，避免估計(jì)偏差過大。

此外，本文考慮了復(fù)發(fā)事件變點(diǎn)模型在有協(xié)變量情形中的應(yīng)用，在以后的工作中可以進(jìn)一步研究變點(diǎn)前后協(xié)變量系數(shù)的變化問題。Li等對(duì)風(fēng)險(xiǎn)函數(shù)的變點(diǎn)模型做了變點(diǎn)前后協(xié)變量系數(shù)發(fā)生變化的假設(shè)，并進(jìn)行了相關(guān)的統(tǒng)計(jì)推斷，在以后的工作中，可以將其應(yīng)用在復(fù)發(fā)事件中進(jìn)行相應(yīng)的研究。

［1］ Wang M C，Qin J，Chiang C T．Analyzing Recurrent Event Data with Informative Censoring［J］．Journal of the American Statistical Association，2001，96（12）．

［2］ Pena E，Stocker R．A General Class of Parametric Models for Recurrent Event Data［J］．Technometrics，Technometrics，2007，49（2）．

［3］ Pena E，Slate E，Gonzalez J．Semiparametric Inference for a General Class of Models for Recurrent Events［J］．Journal of Statistical Planning and Inference，2007，137（6）．

［4］ Matthews D E，F(xiàn)arewell V T．On Testing for Constant Hazard Against a Change－point Alternative［J］．Biometrics：1982，38（2）．

［5］ Chang I S，Chen C H，Hsiung C A．Estimation in Change－point Hazard Rate Models with Random Censorship［C］∥Carlstein E，Muller H G，Siegmund D．Change－point Problems．Hayward：Inst．Math．Statist．，1994．

［6］ Dupuy J F．Estimation in a Change－point Hazard Regression Model［J］．Statistics ＆probability letters，2006，76（2）．

［7］ Dupuy J F．Detecting Change in a Hazard Regression Model with Right－censoring［J］．Journal of Statistical Planning and Inference，2009，139（5）．

［8］ Zhao X，Wu X，Zhou X．A Change－point Model for Survival Data with Long－term Survivors［J］．Statistica Sinica，2009，19（1）．

［9］ Li Y，Qian L，Zhang W．Estimation in a Change－point Hazard Regression Model with Long－term Survivors［J］．Statistics ＆Probability Letters，2013，83（7）．

［10］Frobish D J．Estimation of Change－points in Recurrent Events Models［D］．Doctoral Dissertation：Northern Illinois University，2006．

［11］Andersen P，Gill R，Keiding N．Statistical Models Based on Counting Provesses［M］．New York ：Springer－Verlag，1993．

［12］Cook R J，Lawless J F．The Statistical Analysis of Recurrent Events［M］．New York：Springer，2007．

［13］Frobish D，Ebrahimi N．Parametric Estimation of Change－points for Actual Event Data in Recurrent Events Models［J］．Computational Statistics ＆ Data Analysis，2009，53（3）．

［14］Van Der Vaart A．Asymptotic Statistics［M］．Cambridge：Cambridge University Press，1998．

Parametric Estimation of Change－points with Covariates in Recurrent Events Model

LI Yun-xia，ZHOU Xing-xing
（School of Mathematics and Statistics，Zhejiang University of Finance and Economics，Hangzhou 310018，China）

Recurrent event data is widely applied to biostatistics，clinical medicine，industry reliability and insurance actuarial．The model of change－points in recurrent events has certain significance in the abstract and foreground in the application．We consider the change－points with covariates in recurrent events and adopt MLE to estimate the change－points and parameters．Moreover we illustrate the consistence of the estimators．

change-points；recurrent events；covariate；MLE；steepest ascent

O211

A

1007－3116（2014）07－0011－05

2014－02－18；修得日期：2014－05－31

國家自然科學(xué)基金項(xiàng)目《時(shí)間序列的極限理論及在變點(diǎn)問題上應(yīng)用研究》（10901136）；全國統(tǒng)計(jì)科學(xué)研究計(jì)劃項(xiàng)目《變點(diǎn)問題在生存分析中的應(yīng)用》（2012LY161）；浙江省自然科學(xué)基金項(xiàng)目《生存分析中若干變點(diǎn)模型的研究及其應(yīng)用》（LY14A010022）；浙江省社會(huì)科學(xué)界聯(lián)合會(huì)研究課題《生存數(shù)據(jù)的變點(diǎn)模型研究及其在可靠性理論和生物醫(yī)學(xué)中的應(yīng)用》（2013Z56）；浙江財(cái)經(jīng)大學(xué)校級(jí)研究生科研項(xiàng)目《含有長期生存者的復(fù)發(fā)事件變點(diǎn)問題》（2013YJS075）

李云霞，女，浙江杭州人，統(tǒng)計(jì)學(xué)博士，教授，研究方向：概率極限理論，時(shí)間序列分析，計(jì)量經(jīng)濟(jì)，生存分析；

周杏杏，女，山東濟(jì)南人，碩士生，研究方向：變點(diǎn)問題，生存分析。

（責(zé)任編輯：張治國）