張思琪,韋程東,何國源

(南寧師范大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院,廣西 南寧 530100)

某一事件隨著時間的推移而重復(fù)發(fā)生的過程稱為復(fù)發(fā)事件過程,這種過程提供的數(shù)據(jù)稱為復(fù)發(fā)事件數(shù)據(jù).復(fù)發(fā)事件數(shù)據(jù)廣泛出現(xiàn)在醫(yī)藥公共衛(wèi)生、工商業(yè)、社會科學(xué)和保險等領(lǐng)域,如呼吸病患者的哮喘發(fā)作、新型冠狀病毒肺炎的重復(fù)感染、抑郁癥患者的再次發(fā)作、肝癌肝移植術(shù)后復(fù)發(fā)、同一片區(qū)域多次發(fā)生地震、重要路段交通事故重復(fù)出現(xiàn)等[1].

對于復(fù)發(fā)事件數(shù)據(jù),學(xué)者們根據(jù)協(xié)變量對復(fù)發(fā)事件過程的影響建立不同的模型,主要有兩類.一類是強度模型,針對危險率函數(shù)和強度函數(shù),將死亡事件作為復(fù)發(fā)事件的刪失進行處理,但并未明顯區(qū)分復(fù)發(fā)事件和終止事件.另一類是基于邊際均值或邊際比率模型,不再將終止事件看作獨立刪失,將終止事件與復(fù)發(fā)事件聯(lián)系起來.例如Lin等[1]在不考慮歷史條件的情況下研究了計數(shù)過程的比例均值和比率模型,該模型與強度模型相比更為穩(wěn)健.對于帶終止事件的復(fù)發(fā)事件數(shù)據(jù),一般采用邊際模型方法和隨機效應(yīng)模型方法來處理.

邊際模型方法主要研究復(fù)發(fā)事件和終止事件的邊際模型,不考慮兩者之間的相關(guān)性,如Ghosh[3]通過一個任意計數(shù)過程的模型,提出協(xié)變量改變基線速率函數(shù)的時間尺度的加速速率模型.隨機效應(yīng)模型借助潛變量刻畫復(fù)發(fā)事件和終止事件的相關(guān)性,并假設(shè)在給定脆弱變量時復(fù)發(fā)事件與終止事件獨立,如Wang等[4]通過潛在變量的非平穩(wěn)泊松過程對復(fù)發(fā)事件的發(fā)生進行模擬,將乘法強度模型擴展到半?yún)?shù)回歸模型,提出了估算乘法強度下累積速率函數(shù)和回歸參數(shù)的估算程序.然而在某些情況下,復(fù)發(fā)事件的終止時間可能與復(fù)發(fā)事件的過程相關(guān),不符合獨立刪失的假設(shè),這時就需要考慮復(fù)發(fā)事件過程和失效時間聯(lián)合建模.如Huang等[5]使用共同的潛在變量對復(fù)發(fā)事件過程的強度和失效時間危險之間的關(guān)聯(lián)進行建模,提出“borrow-strength”估計方法,從復(fù)發(fā)事件數(shù)據(jù)中估計潛在變量的值.戴家佳等[6]對帶有終止事件的復(fù)發(fā)事件數(shù)據(jù),提出了加性加速比率回歸模型,運用廣義估計方程的思想給出了參數(shù)的估計.Huang等[7]提出一種既不依賴于脆弱變量的分布,也不依賴于失效時間的一種半?yún)?shù)估計方法,允許協(xié)變量將刪失時間與復(fù)發(fā)事件聯(lián)系在一起,這種方法雖然很靈活,但它是基于復(fù)發(fā)事件過程的泊松分布的,在某些方面無法應(yīng)用.因此Kalbfleisch等[8]在完全沒有泊松過程假設(shè)的情況下進行參數(shù)估計,所得結(jié)果具有極大的穩(wěn)健性,其中的共享脆弱性的估計僅與復(fù)發(fā)事件過程和終止事件過程之間的關(guān)聯(lián)有關(guān).Xu等[9]用共享脆弱性變量對復(fù)發(fā)事件和失效時間建立聯(lián)合尺度變換模型,不需對復(fù)發(fā)事件過程進行泊松假設(shè),回歸參數(shù)具有邊際解釋.Xu等[10]提出了一種既不需要脆弱性分布的參數(shù)假設(shè),也不需要重復(fù)事件過程的泊松假設(shè)估計模型參數(shù).在處理復(fù)發(fā)事件數(shù)據(jù)時,復(fù)發(fā)事件的均值函數(shù)比風(fēng)險函數(shù)具有更強的解釋能力[11].據(jù)文獻所知,目前還未有學(xué)者研究協(xié)變量對于復(fù)發(fā)事件過程的改變時間指數(shù)的尺度效應(yīng).本文所提出的模型與文獻[5]中考慮的聯(lián)合模型相比,并未對給定Z的N(t)條件分布施加泊松型假設(shè).本文的復(fù)發(fā)事件模型與文獻[10]中提出的模型相比,共同點是均包含兩種類型的協(xié)變效應(yīng),不同之處是本文中模型的協(xié)變效應(yīng)是改變時間指數(shù)尺度效應(yīng)和改變風(fēng)險的乘法變化效應(yīng).

本文提出一種帶有指數(shù)變換效應(yīng)的半?yún)?shù)估計,放寬了應(yīng)用條件,無需對復(fù)發(fā)事件過程進行泊松假設(shè),不需要了解未觀察到的脆弱變量的信息,在適當(dāng)?shù)恼齽t條件下,建立了估計量的相合性和漸近正態(tài)性.估計程序使用MATLAB,數(shù)值模擬結(jié)果表明本文提出的估計量幾乎是無偏的,真實值與估計值曲線擬合度高.因此本文提出的半?yún)?shù)模型是有效的.

1 模型建立

設(shè)[0,τ]是觀察復(fù)發(fā)事件發(fā)生的時間段,N(t)為時間區(qū)間[0,t]內(nèi)復(fù)發(fā)事件發(fā)生的次數(shù).將刪失時間Y定義為Y=C∧D,∧表示最小值,C為非信息性刪失時間,D為信息性刪失時間.在研究復(fù)發(fā)事件數(shù)據(jù)時,可以假設(shè)在給定協(xié)變量的基礎(chǔ)上,刪失時間Y與個體所發(fā)生的事件次數(shù)是獨立的.

設(shè)Z是取非負(fù)值的潛在脆弱變量,用來表示復(fù)發(fā)事件過程N(·)和刪失事件(C,D)之間的關(guān)聯(lián).設(shè)X是改變風(fēng)險的乘法效應(yīng)的p維協(xié)變量,W是改變時間指數(shù)尺度效應(yīng)的q維協(xié)變量,注意這兩者可能相同.設(shè)N(·)、C、D是相互獨立的.

以(Z,X,W)為條件,由λ(t)dt=E[dN(t)|Z,X,W]定義的計數(shù)過程N(·)的速率函數(shù)采用以下形式:

(1)

在模型(1)的條件下,定義N(·)的累積均值函數(shù)為Λ(t)=E[N(t)|Z,X,W],即

Λ(t)=ZΛ0(tWTβ+1)exp(XTα).

(2)

模型(2)顯示了W和X的不同效應(yīng).

2 參數(shù)估計

2.1 參數(shù)β的估計

假設(shè)對i=1,…,n,觀察數(shù)據(jù)Xi、Wi、Ni(t),Yi(0≤t≤Yi)是獨立同分布的.設(shè)mi是個體i在時間Yi之前的事件數(shù),如果mi>0,則Ni(t)的跳躍時間為觀測事件時間tij,j=1,…,mi.對于p維向量b,考慮變換

t*(b)=tWTb+1,Y*(b)=YWTb+1.

變換后的復(fù)發(fā)事件的計數(shù)過程為

在模型(1)的基礎(chǔ)上有

特別地,當(dāng)b=β時有

(3)

首先考慮的Λ0(t)一致估計量.對于任意固定的b,定義兩個經(jīng)驗過程Qn(t;b)和Rn(t;b):

Qn(t;b)和Rn(t;b)的極限函數(shù)分別記為Q(t;b)和R(t;b).在模型(1)和(2)的條件下有

由Λ0(τ)=1得

(4)

由此得到Λ0(t)的一致估計量

(5)

前面的討論假設(shè)β是已知的.現(xiàn)在構(gòu)造β的估計方程.定義[0,τ]上的隨機過程

根據(jù)模型(3),對于t∈[0,τ]有

(6)

由(6)的第一式得到H的一致估計量

(7)

2.2 參數(shù)α的估計

在參數(shù)β的估計的基礎(chǔ)上進一步對參數(shù)α進行估計.由(1)和(2)可得

當(dāng)β和Λ0為已知條件時,有估計方程

(8)

3 估計量的漸近性質(zhì)

為了研究所提出的估計量的大樣本性質(zhì), 本文加了如下的正則性條件:

條件2協(xié)變量W和X是一致有界的,并且E(Z2)<∞.

條件3在(X,W,Z)的條件下,Y條件概率密度函數(shù)是連續(xù)且一致有界的.

條件4速率函數(shù)λ0(t)在區(qū)間[0,τ]上恒取非負(fù)值,且有有界的二階導(dǎo)函數(shù).

條件5定理中定義的J和J2是非奇異矩陣.

在上述條件下有以下漸近結(jié)果.

于是有

其中

由泰勒展開公式可得

(9)

我們有以下事實(證明參考[9,Theorem3.1]):

又由于

故由泰勒展開公式得

因此對于‖b-β‖≤dn→0,t∈[0,τ],有如下表示:

進一步有

因為

故由Glivenko定理([12,Theorem 2.4.3])得

其中

又因為

故

故有

證明定義

因此,下列結(jié)果對于考慮的結(jié)果是一致的:

其中

進一步推導(dǎo)Un(γ;β)的漸近正態(tài)性.Un(γ;β)可寫為

這里V(x,y,w,m)是(X,Y,W,m)的聯(lián)合分布函數(shù).因此

給出了n1/2Un(γ;β)的漸近正態(tài)性,其漸近均值為0,協(xié)方差矩陣為E[ci(γ;β)ci(γ;β)T].

4 數(shù)值模擬

本節(jié)通過數(shù)值模擬研究前面提出的估計的樣本性質(zhì),其中樣本偏差(Bias)、樣本均方誤差(MSE)、平均標(biāo)準(zhǔn)誤差(SE)和95%的經(jīng)驗覆蓋概率(CP)用來評估估計結(jié)果.

4.1 數(shù)值模擬參數(shù)設(shè)置

在本次數(shù)值模擬中總共生成500個數(shù)據(jù)集, 每個數(shù)據(jù)集的樣本量n分別設(shè)置為n=100和n=400.對于研究的個體i,設(shè)置協(xié)變量Xi=(Xi1,Xi2),其中Xi1為成功概率為0.5的二項分布,Xi2由均勻分布U=(0,1)獨立生成.脆弱變量Zi設(shè)置為由均值為0.4,方差為0.08的伽馬分布生成Z～Gamma(2,5).假設(shè)協(xié)變量Wi=Xi,復(fù)發(fā)事件過程觀測截止時間τ=10,失效事件時間Di是當(dāng)Xi1=1時由均值為10的指數(shù)分布生成的,當(dāng)Xi1≠1時由均值為100/Zi的指數(shù)分布生成的,刪失時間為Yi,Yi=min(Ci,Di).對于樣本量n,回歸參數(shù)α和β設(shè)為

(α,β)=((-3,3)T,(-0.5,0.5)T), (α,β)=((-3,3)T,(-0.3,0.1)T), (α,β)=((-3,3)T,(-0.2,0.4)T).情形1基準(zhǔn)強度函數(shù)為λ0(t)=t/10.

情形2基準(zhǔn)強度函數(shù)為λ0(t)=log(1+t).

4.2 數(shù)值模擬結(jié)果

圖1 n=100時的擬合曲線

圖2 n=400時的擬合曲線

表2 n=400時的模擬效果

5 總結(jié)

本文拓寬了現(xiàn)有文獻中模型的應(yīng)用條件,提出了一種既不需要脆弱性分布的參數(shù)假設(shè),也不需要重復(fù)事件過程的泊松假設(shè)的半?yún)?shù)估計模型.該模型協(xié)變效應(yīng)分別為改變時間指數(shù)尺度效應(yīng)和改變風(fēng)險的乘法變化效應(yīng).本文對于模型的可識別性,假設(shè)Λ0(τ)=1,同時估計uz.使用廣義估計方程的方法和隨機過程理論對模型參數(shù)進行估計,證明估計參數(shù)的相合性和漸近正態(tài)性.利用MATLAB軟件對估計參數(shù)進行數(shù)值模擬檢驗,驗證了估計參數(shù)的漸近正態(tài)性.數(shù)值模擬結(jié)果表明,本文提出的估計方法表現(xiàn)良好,估計值與真值偏差較小,幾乎無偏.