上海交通大學(xué)生命科學(xué)技術(shù)學(xué)院生物信息與生物統(tǒng)計(jì)系(200240) 秦 飛 俞章盛

【提 要】 目的 建立更為一般化的面板計(jì)數(shù)數(shù)據(jù)模型，既包含具有時(shí)變效應(yīng)的協(xié)變量又包含具有非參數(shù)效應(yīng)的協(xié)變量。方法 使用偽似然樣條估計(jì)方法及全似然樣條估計(jì)方法，并通過(guò)數(shù)值模擬定性和定量地比較這兩種方法。結(jié)果 兩估計(jì)方法在樣本量為200時(shí)都具有很小的偏差和均方誤差(MSE)；全似然樣條估計(jì)方法的偏差和MSE更小(比偽似然樣條估計(jì)方法的MSE小約35%)，但計(jì)算時(shí)間遠(yuǎn)遠(yuǎn)超過(guò)偽似然樣條估計(jì)方法(約是其200倍)；兩種估計(jì)方法用于小兒哮喘數(shù)據(jù)后得到了相似的估計(jì)結(jié)果，揭示了白細(xì)胞介素IL-9的時(shí)變效應(yīng)及白細(xì)胞介素IL-5的非參數(shù)效應(yīng)。結(jié)論 綜合MSE及計(jì)算時(shí)間，在此模型下本研究更推薦使用偽似然樣條估計(jì)方法。

面板計(jì)數(shù)數(shù)據(jù)經(jīng)常出現(xiàn)在電子醫(yī)療病例、臨床試驗(yàn)、流行病學(xué)研究中。對(duì)于面板計(jì)數(shù)數(shù)據(jù)，研究者僅調(diào)查到在兩次相鄰的觀測(cè)時(shí)間點(diǎn)間復(fù)發(fā)性事件發(fā)生的次數(shù)，但未調(diào)查到復(fù)發(fā)性事件發(fā)生的具體時(shí)刻。一個(gè)典型的例子可見(jiàn)于小兒哮喘研究[1]。在此研究中，研究人員對(duì)105個(gè)小兒哮喘病患者連續(xù)跟蹤了5年，通過(guò)電話隨訪獲得截至到此刻小兒哮喘復(fù)發(fā)的次數(shù)，但未調(diào)查到復(fù)發(fā)的具體時(shí)刻，所以收集到的數(shù)據(jù)為面板計(jì)數(shù)數(shù)據(jù)。

在對(duì)面板計(jì)數(shù)數(shù)據(jù)的建模方面，Zhao等[2-3]先后建立了具有時(shí)變效應(yīng)及具有非參數(shù)效應(yīng)協(xié)變量的面板計(jì)數(shù)數(shù)據(jù)模型，但尚未有研究在一個(gè)面板計(jì)數(shù)數(shù)據(jù)模型中同時(shí)考慮到這兩者，而這種情況在實(shí)際應(yīng)用中是可能出現(xiàn)的，比如Cai等[4]就在Cox模型中同時(shí)包括了具有時(shí)變效應(yīng)及具有非參數(shù)效應(yīng)的協(xié)變量。在估計(jì)方面，Zhao等[2-3]僅采用了偽似然函數(shù)下的樣條估計(jì)方法，但沒(méi)有采用全似然函數(shù)下的樣條估計(jì)方法。因此，本文將建立同時(shí)含有時(shí)變效應(yīng)及非參數(shù)效應(yīng)協(xié)變量的面板計(jì)數(shù)數(shù)據(jù)模型，采用偽似然樣條估計(jì)及全似然樣條估計(jì)的方法并對(duì)其進(jìn)行對(duì)比，最后將此方法應(yīng)用于小兒哮喘研究中[1]。

模型及估計(jì)

1.模型

對(duì)于面板計(jì)數(shù)數(shù)據(jù)，研究者可以觀測(cè)到n個(gè)獨(dú)立的樣本Ui=(Ki，Ti，N(i)，Xi，Zi)，i=1，2，…，n，其中Ki為觀測(cè)次數(shù)，Ti={Ti，1，Ti，2，…，Ti，Ki}為觀測(cè)時(shí)間點(diǎn)，N(i)={N(Ti，1)，N(Ti，2)，…，N(Ti，Ki)}為面板計(jì)數(shù)值，Xi=(Xi，1，Xi，2，…，Xi，p)T為具有時(shí)變效應(yīng)的協(xié)變量，Zi=(Zi，1，Zi，2，…，Zi，q)T為具有非參數(shù)效應(yīng)的協(xié)變量。類似于其他非參數(shù)模型，本文假定某個(gè)協(xié)變量是具有時(shí)變效應(yīng)還是非參數(shù)效應(yīng)是提前設(shè)定好的。本文建立以下非參數(shù)模型：

E[N(i)(t)]=Λ0(t)·exp[(β1(t)T·Xi+β2(Zi)]

(1)

其中，Λ0(t)為單調(diào)非減的非負(fù)函數(shù)，β1(t)=(β1，1(t)，β1，2(t)，…，β1，p(t))T為時(shí)變系數(shù)，β2(Zi)=β2，1(Zi，1)+β2，2(Zi，2)+…+β2，q(Zi，q)為非參數(shù)效應(yīng)。本文僅呈現(xiàn)當(dāng)p=q=1時(shí)的情形，本文的估計(jì)方法很容易拓展到p，q>1的情形。

2.樣條函數(shù)

(2)

(3)

(4)

3.偽似然估計(jì)

正如Zhang等[5-6，14-15]所述，首先假設(shè)潛在的計(jì)數(shù)過(guò)程{N(t)：t≥0}為非齊次的泊松過(guò)程，然后忽略每個(gè)個(gè)體的計(jì)數(shù)數(shù)據(jù){N(Ti，1)，N(Ti，2)，…，N(Ti，Ki)}之間的相關(guān)性從而得到的偽似然估計(jì)量通常具有漸近正態(tài)性和相合性。參照Z(yǔ)hang等[5-6，14-15]的思路，本文建立模型(1)的偽似然函數(shù)如下：

L(1)=P[N(T1，1)=N1，1，N(T1，2)=N1，2，…，N(T1，K1)=N1，K1，…，N(Tn，1)=Nn，1，…]

對(duì)L(1)取對(duì)數(shù)并忽略無(wú)關(guān)項(xiàng)，然后將樣條近似式(2)～(4)代入，得到以下對(duì)數(shù)偽似然函數(shù)：

(5)

4.全似然估計(jì)

如Zhang等[5-6，14-15]所述，全似然估計(jì)量通常要比偽似然估計(jì)量估得更準(zhǔn)但同時(shí)計(jì)算量更大。本文也探討了全似然估計(jì)在模型(1)下的表現(xiàn)。全似然估計(jì)量是通過(guò)首先假設(shè)潛在的計(jì)數(shù)過(guò)程為非齊次泊松過(guò)程，然后利用面板計(jì)數(shù)數(shù)據(jù)增量間的獨(dú)立性而構(gòu)造。本文建立模型(1)的全似然估計(jì)函數(shù)如下：

exp[-[Λ0(Ti，j)·exp(β1(Ti，j)·Xi+β2(Zi))-Λ0(Ti，j-1)·

exp(β1(Ti，j-1)·Xi+β2(Zi))]]}

其中，Λ0(Ti，0)≡0，Ni，0≡0。對(duì)L(2)取對(duì)數(shù)忽略并無(wú)關(guān)項(xiàng)，然后將樣條近似式(2)～(4)代入，得到以下對(duì)數(shù)全似然函數(shù)：

(6)

同樣地，極大化式(6)便可得到參數(shù)的全似然估計(jì)。這一過(guò)程仍可通過(guò)相同的R函數(shù)constrOptim()來(lái)實(shí)現(xiàn)。

模擬研究

N(Ti，j)-N(Ti，j-1)～Po{Λ0(Ti，j)·exp[β1(Ti，j)·Xi+β2(Zi)]-Λ0(Ti，j-1)·exp[β1(Ti，j-1)·Xi+β2(Zi)]}

其中j=1，…，Ki，Ti，0≡0，N(Ti，0)≡0，Λ0(0)≡0。本文對(duì)真實(shí)函數(shù)設(shè)置以下兩種情形：

情形1：Λ0(t)=t+1，β1(t)=1.5·sin(0.05·πt)，β2(z)=sin(πz)；

情形2：Λ0(t)=t+1，β1(t)=0.15·t，β2(z)=sin(2πz)·I(z≤0.5)+0.5·sin(2πz)，其中I(·)為示性函數(shù)?？梢钥吹角樾?要比情形1更為復(fù)雜。

本文設(shè)置樣本量為50及200，產(chǎn)生500次蒙特卡洛數(shù)據(jù)。參照Lu等[5-6]的做法，本文使用三次樣條，樣條的內(nèi)部節(jié)點(diǎn)數(shù)設(shè)置為6，采用分位數(shù)的方法放置節(jié)點(diǎn)，即選擇所有不同觀測(cè)時(shí)間點(diǎn)的k/(m+1)分位點(diǎn)(k=0，1，…，m+1)為這m+2個(gè)節(jié)點(diǎn)的放置位置。

圖1展示了情形1下樣本量為50時(shí)對(duì)這三個(gè)函數(shù)Λ0(t)，β1(t)，β2(z)500次估計(jì)的均值曲線及2.5%，97.5%分位數(shù)曲線。從圖1可以看出，這兩種方法的估計(jì)結(jié)果都存在一定的偏離，也如前文所預(yù)料，這三個(gè)函數(shù)的全似然估計(jì)的均值曲線更接近于真實(shí)曲線，且置信區(qū)間更窄。當(dāng)樣本量增大為200時(shí)，所有的均值曲線相對(duì)樣本量為50時(shí)都更接近真實(shí)曲線且置信區(qū)間更窄，全似然估計(jì)方法此時(shí)仍舊估得更準(zhǔn)且置信區(qū)間更窄，但同時(shí)也注意到此時(shí)兩方法的均值曲線都幾乎跟真實(shí)曲線重合(圖2)。

圖1 情形1，樣本量為50時(shí)的估計(jì)結(jié)果

圖2 情形1，樣本量為200時(shí)的估計(jì)結(jié)果

本文也對(duì)這兩種估計(jì)量進(jìn)行了定量比較。從表1可以看出，當(dāng)樣本量為50時(shí)，全似然估計(jì)下這三個(gè)函數(shù)的估計(jì)偏差、均方誤差均小于偽似然估計(jì)，但偽似然估計(jì)的計(jì)算時(shí)間遠(yuǎn)遠(yuǎn)小于全似然估計(jì)。當(dāng)樣本量增大為200時(shí)，這兩種估計(jì)下函數(shù)的估計(jì)偏差、均方誤差都大大減小，全似然估計(jì)方法仍估得更準(zhǔn)，但同時(shí)計(jì)算時(shí)間方面?zhèn)嗡迫还烙?jì)方法仍遠(yuǎn)具優(yōu)勢(shì)，并且也注意到此時(shí)兩估計(jì)方法的偏差幾乎可以忽略。情形2的模擬結(jié)果和結(jié)論類似情形1，由于篇幅有限，本文沒(méi)有展示。由于在接下來(lái)的實(shí)際應(yīng)用中樣本量為105，本文還做了上述兩情形下在樣本量為100時(shí)的模擬，此時(shí)的估計(jì)偏差、均方誤差及計(jì)算時(shí)間介于樣本量為50時(shí)的結(jié)果和樣本量為200時(shí)的結(jié)果之間，得到的結(jié)論也類似(由于篇幅有限，未展示)。總之，就估計(jì)的準(zhǔn)確性和穩(wěn)定性而言，全似然估計(jì)方法優(yōu)于偽似然估計(jì)方法，前者的TMSE比后者小約35%。但就計(jì)算時(shí)間而言，偽似然估計(jì)占絕對(duì)優(yōu)勢(shì)(全似然估計(jì)計(jì)算時(shí)間約為偽似然的200倍)，再加上當(dāng)樣本量為200時(shí)，偽似然估計(jì)量已非常接近真實(shí)函數(shù)，所以在模型(1)下，本文推薦使用偽似然估計(jì)方法。

表1 情形1下兩種估計(jì)量的定量比較結(jié)果

實(shí)例應(yīng)用

本部分把前文介紹的方法應(yīng)用到小兒哮喘研究中，該研究的總體描述可參考文獻(xiàn)[1]。該研究旨在探究免疫因子特征和小兒哮喘之間的關(guān)系，納入了105名小兒哮喘患者，平均入組年齡為10.9月，50.5%的患者是女性，9.5%的患者母親在懷孕時(shí)抽煙。

白細(xì)胞介素IL-9是CD4+輔助細(xì)胞分泌的一種細(xì)胞因子，對(duì)哮喘小鼠模型的遺傳學(xué)研究表明，該細(xì)胞因子是支氣管高反應(yīng)性發(fā)病的決定性因素[16]。根據(jù)臨床經(jīng)驗(yàn)，有很多因素都會(huì)影響到小兒哮喘而且它們的影響效果會(huì)隨著時(shí)間的變化而變化，因此本文設(shè)定IL-9具有時(shí)變效應(yīng)。此外，另一種白細(xì)胞介素IL-5一直與過(guò)敏性鼻炎、哮喘等多種變應(yīng)性疾病相關(guān)[17]。根據(jù)Zhao等[3]的研究結(jié)果，IL-5具有明顯的非參數(shù)效應(yīng)。因此，在本文的面板計(jì)數(shù)數(shù)據(jù)模型中包括性別及是否吸煙這兩個(gè)協(xié)變量后，本文用偽似然估計(jì)和全似然估計(jì)這兩種方法來(lái)估計(jì)IL-9的時(shí)變效應(yīng)及IL-5的非參數(shù)效應(yīng)。

本文使用三次樣條來(lái)估計(jì)模型中的未知函數(shù)；使用AIC準(zhǔn)則(AIC=-2·l(η)+2k)來(lái)選擇內(nèi)部節(jié)點(diǎn)個(gè)數(shù)，其中l(wèi)(η)為對(duì)數(shù)偽似然函數(shù)值或者對(duì)數(shù)全似然函數(shù)值，k為模型中參數(shù)的個(gè)數(shù)；和模擬研究類似，節(jié)點(diǎn)的位置選取依舊使用分位數(shù)的方法。對(duì)于偽似然估計(jì)，本研究選擇的節(jié)點(diǎn)數(shù)為7，對(duì)于全似然估計(jì)，選擇到的節(jié)點(diǎn)數(shù)為6，對(duì)時(shí)變效應(yīng)β1(t)和非參數(shù)效應(yīng)β2(z)的估計(jì)結(jié)果如圖3所示。從圖3可以看出，由于使用的內(nèi)部節(jié)點(diǎn)數(shù)更多，偽似然估計(jì)方法對(duì)β1(t)及β2(z)的估計(jì)曲線更震蕩些，但總體趨勢(shì)和全似然估計(jì)方法得到的曲線十分相似。這兩種估計(jì)結(jié)果都揭示了IL-9對(duì)小兒哮喘效應(yīng)的時(shí)變性，即在30個(gè)月齡之前，效應(yīng)隨著年齡的增大而增大，但其后效應(yīng)基本趨于平穩(wěn)。同之前Zhao等[3]的研究結(jié)果一致，IL-5的效應(yīng)仍舊是非參數(shù)的，效應(yīng)隨著IL-5值的增大而增大。本研究中模擬及實(shí)例應(yīng)用部分是在R軟件中實(shí)施的，實(shí)例應(yīng)用部分的代碼可通過(guò)以下ftp地址獲得：ftp：//public.sjtu.edu.cn/(用戶名yuzhangsheng，密碼public)。

圖3 小兒哮喘研究中對(duì)白細(xì)胞介素IL-9及IL-5的估計(jì)結(jié)果

小 結(jié)

本文建立了更為一般化的面板計(jì)數(shù)數(shù)據(jù)模型，其既包含具有時(shí)變效應(yīng)的協(xié)變量，也包含具有非參數(shù)效應(yīng)的協(xié)變量；本文使用了偽似然樣條估計(jì)方法和全似然樣條估計(jì)方法，并將這兩種估計(jì)方法進(jìn)行了對(duì)比；最后分析了小兒哮喘研究數(shù)據(jù)。

在模擬研究部分，發(fā)現(xiàn)當(dāng)樣本量為50和100時(shí)，這兩種估計(jì)方法都存在一定的偏差，但全似然樣條估計(jì)方法的偏差和均方誤差更小；同時(shí)由于全似然函數(shù)更為復(fù)雜，所以全似然樣條估計(jì)方法具有更長(zhǎng)的計(jì)算時(shí)間。當(dāng)樣本量增大到200時(shí)，兩種估計(jì)方法都具有可以忽略的偏差，且全似然樣條估計(jì)方法仍舊具有更小的偏差和均方誤差，但同時(shí)計(jì)算時(shí)間也長(zhǎng)得多。綜合均方誤差和計(jì)算時(shí)間，且考慮到當(dāng)樣本量為200時(shí)偽似然樣條估計(jì)方法已較準(zhǔn)確，本文推薦在類似的模型中使用偽似然樣條估計(jì)方法。

在實(shí)例應(yīng)用部分，這兩種估計(jì)方法得到了類似的估計(jì)結(jié)果，都揭示了白細(xì)胞介素IL-9的時(shí)變效應(yīng)及白細(xì)胞介素IL-5的非參數(shù)效應(yīng)，此結(jié)論也跟之前的研究[3]一致。探究在本文模型下兩種估計(jì)方法的理論性質(zhì)，研究更為高效的算法將是我們未來(lái)的研究方向。