含缺失協(xié)變量的混合效應(yīng)模型的多元單邊檢驗(yàn)

齊培艷，仵旭靚，段西發(fā)

（太原科技大學(xué) 應(yīng)用科學(xué)學(xué)院，山西 太原，030024）

0 引言

人類免疫缺陷病毒（Human Immunodeficiency Virus，HIV）病毒動力學(xué)研究是近年來艾滋病研究的熱門，非線性混合效應(yīng)（NLME）模型［1-2］、半?yún)?shù)非線性混合效應(yīng)（SNLME）模型［3］被廣泛應(yīng)用于HIV病毒動力學(xué)、藥代動力學(xué)分析和消長研究。這些模型通?；谏蓴?shù)據(jù)的底層機(jī)制，其參數(shù)具有一定意義的物理解釋，如抗HIV治療期間病毒衰減率或隨機(jī)效應(yīng)的方差等參數(shù)應(yīng)該是非負(fù)的，要滿足這些條件就需要對參數(shù)進(jìn)行自然限制或約束。因此，應(yīng)該在包含限制的前提下對參數(shù)進(jìn)行單邊或序約束檢驗(yàn)，在假設(shè)檢驗(yàn)中加入自然參數(shù)限制或約束會產(chǎn)生比相應(yīng)的無約束檢驗(yàn)更有效的結(jié)果［4］。

約束統(tǒng)計(jì)推斷一直以來都受到極大的關(guān)注，Silvapulle和Sen［5］對其研究進(jìn)展進(jìn)行了綜述。然而，大多數(shù)文獻(xiàn)關(guān)注的是方差分量和均值的檢驗(yàn)［6-7］，對回歸模型的單邊檢驗(yàn)的研究較為有限。Brice等［8］對潛變量模型（LVMs）進(jìn)行了修正的 Wald檢驗(yàn)；Verbeke和 Molenberghs［9］討論了線性混合效應(yīng)（LME）模型的約束推理；Qeadan等［10］證明了一類線性混合模型的廣義似然比檢驗(yàn)（LRT）和Wald檢驗(yàn)在檢驗(yàn)方差分量時(shí)的等價(jià)性；Zhou等［11］從HIV病毒動力學(xué)模型出發(fā)，提出了具有截尾響應(yīng)的NLME模型的基于近似似然的多元單邊檢驗(yàn)，仿真結(jié)果表明所提出的檢驗(yàn)比相應(yīng)的雙邊檢驗(yàn)或無約束檢驗(yàn)更有效；Wang［12］對NLME模型中的多元單邊檢驗(yàn)發(fā)展了似然比檢驗(yàn)（LRT）和Wald檢驗(yàn)，并證明了它們的有效性要強(qiáng)于傳統(tǒng)的雙邊檢驗(yàn)。本文主要考慮SNLME模型中固定參數(shù)的單邊假設(shè)檢驗(yàn)。

在實(shí)踐中，數(shù)據(jù)往往會在應(yīng)答或協(xié)變量中出現(xiàn)缺失，使用完全情況（CC）方法丟棄所有的不完全觀測數(shù)據(jù)，可能會使結(jié)果出現(xiàn)偏差［13］。為了減少這種誤差，Wang等［14］和Zhou等［15］分別討論了基于多重插補(bǔ)和經(jīng)驗(yàn)似然的具有缺失協(xié)變量的NLME模型的多元單邊檢驗(yàn)。本文針對具有缺失協(xié)變量的半?yún)?shù)非線性混合效應(yīng)模型中固定參數(shù)的多元單邊檢驗(yàn)問題，提出了基于多重插補(bǔ)的似然比檢驗(yàn)（LRT）和Wald檢驗(yàn)，并通過實(shí)例分析對HIV病毒動力學(xué)研究中雙相指數(shù)病毒衰變模型的病毒衰減率是否嚴(yán)格為正進(jìn)行了檢驗(yàn)。這里假設(shè)缺失數(shù)據(jù)是隨機(jī)缺失的，即缺失可能與觀察到的數(shù)據(jù)有關(guān)，而不是與缺失的數(shù)據(jù)有關(guān)［16］。

本文第1節(jié)對SNLME模型及多元單邊檢驗(yàn)進(jìn)行了介紹；第2節(jié)討論了具有缺失協(xié)變量的SNLME模型的多重插補(bǔ)方法；第3節(jié)，基于多重插補(bǔ)的多個(gè)估計(jì)提出了兩種檢驗(yàn)：似然比檢驗(yàn)（LRT）和Wald檢驗(yàn)；第4節(jié)通過HIV實(shí)例分析對提出的方法進(jìn)行了有效性檢驗(yàn)；第5節(jié)對本文進(jìn)行了討論和總結(jié)。

1 SNLME模型及多元單邊檢驗(yàn)

1.1 SNLME模型

其中μp和ξiq分別是固定系數(shù)和隨機(jī)系數(shù)的未知向量?；趆i(t)的假設(shè)，可以將ξiq看作獨(dú)立同分布零均值隨機(jī)向量的實(shí)現(xiàn)?？紤]基于百分位結(jié)點(diǎn)的自然三次樣條基，通過Akaike信息準(zhǔn)則（AIC）或貝葉斯信息準(zhǔn)則（BIC）來確定結(jié)點(diǎn)數(shù)（即最佳模型）。

用wp(t)和hiq(t)代替w(t)和hi(t)，則SNLME模型可以用下面的參數(shù)NLME模型來近似：

1.2 多元單邊檢驗(yàn)

2 具有缺失協(xié)變量的SNLME模型的多重插補(bǔ)

在HIV病毒動力學(xué)研究中，CD4細(xì)胞計(jì)數(shù)（協(xié)變量）往往不能完全觀察到，這就會導(dǎo)致協(xié)變量存在測量誤差和缺失，因此需要對協(xié)變量過程進(jìn)行建模。本文考慮以下LME模型來描述協(xié)變量過程：

2.1 貝葉斯框架下的多重插補(bǔ)

設(shè)θ=(α，ψ)≡(α，β，δ，R，A，D)是模型（3）和（5）未知參數(shù)的集合，用a和b表示協(xié)變量模型（5）中的隨機(jī)效應(yīng)ai和響應(yīng)模型（3）中的bi的集合。設(shè)f(·|·)和π(·)分別為條件密度函數(shù)和先驗(yàn)密度函數(shù)，并設(shè)(X|Y)表示給定Y的X的條件分布，則未觀測到的“真”協(xié)變量值=Uijα+Vijai完全由平均參數(shù)α和隨機(jī)效應(yīng)ai決定，的多重插補(bǔ)可 (α，ai)從的多重插補(bǔ)中獲得，可由以下給定觀測數(shù)據(jù)的(α，a)的后驗(yàn)預(yù)測分布生成：

一般上式積分是高維的，計(jì)算十分復(fù)雜，因此，不能直接根據(jù)觀測數(shù)據(jù)計(jì)算(α，a)的后驗(yàn)預(yù)測分布。作為替代方法，馬爾可夫鏈蒙特卡羅方法（MCMC）方法可用于從表達(dá)式（6）中取樣，且不需要計(jì)算積分，同時(shí)可以利用Gibbs抽樣和Metropolis-Hasting（M-H）算法來實(shí)現(xiàn)MCMC方法。

為了避免計(jì)算表達(dá)式（6）中的積分，從(α，a)的后驗(yàn)預(yù)測分布中生成所需的(α，a)的估算，這里使用數(shù)據(jù)增強(qiáng)算法引入潛在變量(b，ψ)。具體來說，從初始值開始，通過迭代以下兩步來創(chuàng)建 (α，a)的估算：對于t=1，2，…，

2.2 算法的實(shí)現(xiàn)

對于上一小節(jié)提出的貝葉斯框架下的多重插補(bǔ)方法，需要指定參數(shù)θ的先驗(yàn)分布。首先假設(shè)參數(shù) α，β，δ，R，A，D 是獨(dú)立的，即 π(θ)≡π(α，β，δ，R，A，D)=π(α)π(β)π(δ)π(R)π(A)π(D)，其中平均參數(shù)α和β服從正態(tài)先驗(yàn)分布，而精度參數(shù)δ，R，A和D服從伽馬（Gamma）或威沙特（Wishart）先驗(yàn)分布，如下所示：

其中，超參數(shù) γ，λ，p，q，r，Σα，Σβ，ΣR，ΣA，ΣD已知，Gamma 分布 G(γ，λ)和 Wishart分布 W(ΣR，p)的均值分別為γλ和pΣR。在實(shí)踐中，先驗(yàn)信息可以從以前的研究或參考文獻(xiàn)中獲得。對于某些具有可靠先驗(yàn)信息的參數(shù)，可以使用方差較小的強(qiáng)先驗(yàn)；對于其他沒有足夠先驗(yàn)信息的參數(shù)，則可以采用非信息性先驗(yàn)分布（差異較大）。

在確定了觀測數(shù)據(jù)的模型和模型中未知參數(shù)的先驗(yàn)分布后，可以使用數(shù)據(jù)增強(qiáng)算法和Gibbs抽樣與M-H算法模擬多個(gè)獨(dú)立的條件分布(α，a)。具體來說，使用Gibbs抽樣，從條件分布f(α|D，a，b，ψ)，f(ai|yi，zi，bi，θ)，f(bi|yi，zi，ai，θ)，(i=1，…，n)和 f(ψ|D，α，a，b)中迭代樣本。

用 β-表示參數(shù) θ中除去 β 外其他的參數(shù)，δ-2，A-1，R-1，D-1同理，則精度參數(shù) δ-2，A-1，R-1，D-1的全條件分布為：

3 具有缺失協(xié)變量的SNLME模型的多元單邊假設(shè)檢驗(yàn)

如第2節(jié)所討論的，將SNLME模型通過3次樣條基函數(shù)近似的參數(shù)NLME模型生成適當(dāng)?shù)亩嘀夭逖a(bǔ)之后，我們獲得了m個(gè)完整的數(shù)據(jù)集，接下來將這m個(gè)結(jié)果進(jìn)行組合便可產(chǎn)生總體推斷。對于假設(shè)檢驗(yàn)問題這里提供了兩種方法來組合多個(gè)完整數(shù)據(jù)分析結(jié)果：1.可以組合多個(gè)完整數(shù)據(jù)檢驗(yàn)統(tǒng)計(jì)量，即LRT統(tǒng)計(jì)量，然后提出總體檢驗(yàn)；2.可以組合完整數(shù)據(jù)的參數(shù)估計(jì)，然后基于組合估計(jì)提出單個(gè)檢驗(yàn)統(tǒng)計(jì)量，即Wald型檢驗(yàn)統(tǒng)計(jì)量。

4 HIV實(shí)例分析

本節(jié)采用SNLME模型來建立HIV病毒動力學(xué)模型，經(jīng)驗(yàn)多項(xiàng)式LME模型來建立協(xié)變量過程。將所提出的基于多重插補(bǔ)的檢驗(yàn)TLR和TW應(yīng)用于真實(shí)數(shù)據(jù)的多元單邊檢驗(yàn)問題中，并與樸素CC法進(jìn)行比較。

這些數(shù)據(jù)來自包含45名HIV感染患者的抗HIV治療研究，他們在研究開始時(shí)接受了一種有效的抗逆轉(zhuǎn)錄病毒方案的治療，在開始治療后第0，2，7，10，14，28天和第6，8，12周測量病毒載量（RNA），在整個(gè)研究中也記錄了CD4細(xì)胞數(shù)目，每位患者的觀察次數(shù)從3到9不等，一共觀察了296次。其中，缺失數(shù)據(jù)出現(xiàn)在隨時(shí)間變化的協(xié)變量CD4中，缺失率大概10%左右。圖1顯示了45名病人治療期間病毒載量的軌跡，可以看到，不同患者之間的軌跡不同，但病毒載量的總體趨勢都隨著時(shí)間的推移在下降，這表明抗HIV的治療可能是有效的。為了避免估計(jì)的不穩(wěn)定性（過大或過小估計(jì)），本文標(biāo)準(zhǔn)化了CD4計(jì)數(shù)并重新定標(biāo)原始時(shí)間（以天為單位），以使新的時(shí)間刻度在0和1之間。

圖1 抗HIV治療期間病人的病毒載量軌跡Fig.1 Viral load trajectory of patients during anti-HIV treatment

考慮以下HIV病毒動力學(xué)模型來描述病毒載量：

采用AIC或BIC準(zhǔn)則選擇最佳模型。使用自然三次樣條基與百分位數(shù)結(jié)點(diǎn)的線性組合來逼近非參數(shù)函數(shù)w(t)和hi(t)，令ψ0(t)=?0(t)≡1，取與（2）式中相同的三次樣條基函數(shù)，其中q≤p。采用AIC/BIC準(zhǔn)則由表1可確定模型（2）中p=3，q=1即：

表1 非參數(shù)函數(shù)w(t)和hi(t)的模型選擇Table 1 Model selection for nonparametric functions w(t)and hi(t)

模 型 中 固 定 效 應(yīng) 為 β=(P1，P2，λ1，λ2，β1，β2，β3)T，隨 機(jī) 效 應(yīng) 為 bi=(b1i，b2i，b3i，b4i)T，假 設(shè)eiji.i.d～N(0，σ2I)，bi～N(0，D)，隨機(jī)誤差和隨機(jī)效應(yīng)相互獨(dú)立。表示真實(shí)的（但不可觀察的）CD4計(jì)數(shù)。眾所周知的是，CD4計(jì)數(shù)的測量通常有很大的誤差。因此，這里假設(shè)初始病毒衰減率與真實(shí)的CD4值有關(guān)，而不是與觀察到的CD4值有關(guān)。

對于協(xié)變量CD4，本文采用經(jīng)驗(yàn)多項(xiàng)式LME模型模擬其軌跡：

在HIV病毒動力學(xué)研究中，研究人員經(jīng)常感興趣的是治療的有效性，它可以通過真實(shí)衰減率是否嚴(yán)格為正來表明結(jié)果。模型中λ1，λ2分別表示生產(chǎn)性感染細(xì)胞，長期性感染細(xì)胞和潛伏性感染細(xì)胞的病毒衰減率，由于第一階段為快速指數(shù)衰變階段，而第二階段衰變緩慢［17］，因此，在建立模型時(shí)，可以假設(shè) P1＞P2，λ1＞λ2，以使模型具有可辨識性。

對于HIV病毒動力學(xué)模型（10），考慮多元單邊檢驗(yàn)問題（4），假設(shè)D是對角矩陣，則可以得到所有參數(shù)的集合為 θ=(P1，P2，λ1，λ2，β1，β2，β3，δ2，d11，d22，d33，d44)T，因此，矩陣R有以下形式：

其中矩陣R的秩r=2。

這樣，針對病毒衰減率是否嚴(yán)格為正的檢驗(yàn)問題可以寫成如下形式：

接下來，使用組合LRT檢驗(yàn)和總體Wald檢驗(yàn)對參數(shù)進(jìn)行檢驗(yàn)，并得出結(jié)論。

為了解決協(xié)變量缺失問題，基于MCMC迭代，去掉初始的5萬次迭代之后，再從接下來的15萬次迭代中每500個(gè)MCMC抽樣中選取一個(gè)樣本，這樣就得到了300個(gè)用于統(tǒng)計(jì)推斷的樣本。然后基于第四節(jié)的介紹，為協(xié)變量CD4生成m=5個(gè)多重插補(bǔ)數(shù)據(jù)集，并對這5個(gè)數(shù)據(jù)集進(jìn)行參數(shù)估計(jì)。表2給出了基于估計(jì)結(jié)果得到的兩個(gè)感興趣參數(shù)的最大似然估計(jì)及均值，可以看到，正如之前假設(shè)的，患者第一階段的病毒衰減率遠(yuǎn)大于第二階段。

表2 五組多重插補(bǔ)結(jié)果的估計(jì)值及其平均值Table 2 Estimated and average values of five groups of multiple imputation

按照第四節(jié)給出的檢驗(yàn)統(tǒng)計(jì)量計(jì)算可以得到組合LRT檢驗(yàn)統(tǒng)計(jì)量為TLR=638.7，總體Wald檢驗(yàn)統(tǒng)計(jì)量為TW=621.1，且兩個(gè)檢驗(yàn)的P值均接近于0。因此，有強(qiáng)有力的證據(jù)拒絕H0，從而得出結(jié)論：在抗HIV治療期間，病毒載量顯著下降。而基于CC方法的檢驗(yàn)統(tǒng)計(jì)量為280.61，它的結(jié)果也拒絕了H0，但是比本文提出的檢驗(yàn)統(tǒng)計(jì)量TLR和TW要弱。

5 結(jié)論

多重插補(bǔ)方法［18］的主要優(yōu)點(diǎn)就是可以對輸入數(shù)據(jù)集使用現(xiàn)有的完全數(shù)據(jù)方法，并且可以對輸入數(shù)據(jù)集執(zhí)行各種完全數(shù)據(jù)假設(shè)檢驗(yàn)。本文基于多重插補(bǔ)方法對具有缺失協(xié)變量的半?yún)?shù)非線性混合效應(yīng)模型的多元單邊檢驗(yàn)問題進(jìn)行了討論，給出了兩個(gè)檢驗(yàn)統(tǒng)計(jì)量，似然比檢驗(yàn)（LRT）統(tǒng)計(jì)量和Wald檢驗(yàn)統(tǒng)計(jì)量，并推導(dǎo)了其漸近零分布。并且為了驗(yàn)證檢驗(yàn)統(tǒng)計(jì)量的性能，將給出的統(tǒng)計(jì)量應(yīng)用到真實(shí)的HIV治療研究中，通過計(jì)算我們得出了三個(gè)檢驗(yàn)統(tǒng)計(jì)量的值，實(shí)例分析表明基于多重插補(bǔ)的似然比檢驗(yàn)（LRT）統(tǒng)計(jì)量和Wald檢驗(yàn)統(tǒng)計(jì)量給出的結(jié)果相似，而且比樸素的CC方法更提供了更強(qiáng)有力的證據(jù)。

在實(shí)踐中，當(dāng)缺失數(shù)據(jù)是信息性的或不可忽略的時(shí)候，必須假設(shè)缺失數(shù)據(jù)機(jī)制模型，并將缺失數(shù)據(jù)模型合并到聯(lián)合模型中。同時(shí)，多元單邊檢驗(yàn)問題可以推廣到其他混合效應(yīng)模型，如廣義線性模型、非參數(shù)混合效應(yīng)模型、Logistic回歸模型、比例風(fēng)險(xiǎn)模型等。這些問題超出了本文的范圍，可以作為下一步的研究。

當(dāng)然本文使用的SNLME模型不只適用于抗HIV病毒的研究，也可用于其他病毒治療或藥物研究中，如對結(jié)核病患者藥物治療的建模，糖尿病患者血漿中亮氨酸動力學(xué)的群體模型等。