多響應(yīng)線性混合效應(yīng)模型的A- 最優(yōu)設(shè)計(jì)

段冬艷， 劉 欣

(東華大學(xué) 理學(xué)院， 上海 201620)

混合效應(yīng)模型作為分析縱向數(shù)據(jù)、重復(fù)測量數(shù)據(jù)和面板數(shù)據(jù)等的常用工具，在諸多領(lǐng)域得到了應(yīng)用。從畜牧養(yǎng)殖到群體藥代動(dòng)力學(xué)、從空間統(tǒng)計(jì)到個(gè)體化醫(yī)學(xué)都可以找到成功應(yīng)用的案例[1-3]。隨著計(jì)算機(jī)技術(shù)的發(fā)展，對混合效應(yīng)模型的統(tǒng)計(jì)分析和應(yīng)用都取得了長足的進(jìn)展。由文獻(xiàn)[4-6]可知，目前混合效應(yīng)模型的性質(zhì)已經(jīng)得到了很好研究。

在許多實(shí)際問題中，響應(yīng)向量往往同時(shí)具有多個(gè)組成部分，并且這些組成部分可能相互關(guān)聯(lián)，研究者須對這些部分進(jìn)行聯(lián)合建模，多響應(yīng)混合效應(yīng)模型是分析這類數(shù)據(jù)的重要工具。例如：Dahm等[7]采用多響應(yīng)混合效應(yīng)模型分析動(dòng)物育種試驗(yàn)；Sun等[8]利用兩響應(yīng)混合效應(yīng)模型研究1973—1977年在烏普薩拉市出生的獨(dú)生子女的生長曲線；Verbeke等[6]采用多響應(yīng)線性混合效應(yīng)模型分析收縮壓和舒張壓的試驗(yàn)數(shù)據(jù)；Jensen等[9]采用多響應(yīng)線性混合效應(yīng)模型分析空氣質(zhì)量與辦公室工作績效之間關(guān)系的試驗(yàn)數(shù)據(jù)。

混合效應(yīng)模型的分析結(jié)果與獲取數(shù)據(jù)的試驗(yàn)方案有著密切的聯(lián)系，如試驗(yàn)中變量的觀測時(shí)間和藥物的劑量等試驗(yàn)條件的設(shè)置都會(huì)直接影響最終統(tǒng)計(jì)推斷的效力。如何有效安排試驗(yàn)是試驗(yàn)設(shè)計(jì)的目的。最優(yōu)試驗(yàn)設(shè)計(jì)方法作為試驗(yàn)設(shè)計(jì)理論的一個(gè)重要分支，正被越來越多地應(yīng)用于混合效應(yīng)模型的試驗(yàn)設(shè)計(jì)研究中。在單響應(yīng)混合效應(yīng)模型的最優(yōu)設(shè)計(jì)方面：Entholzner等[10]研究得出在混合模型下的最優(yōu)和最有效的設(shè)計(jì)；Schmelter[11-12]考慮了隨機(jī)系數(shù)模型中單群恒等設(shè)計(jì)的最優(yōu)性和多群恒等設(shè)計(jì)的最優(yōu)性問題；文獻(xiàn)[13-15]分別對隨機(jī)截距模型、隨機(jī)斜率模型和三次回歸隨機(jī)系數(shù)模型的優(yōu)化設(shè)計(jì)進(jìn)行了研究；文獻(xiàn)[16-17]討論了具有隨機(jī)截距項(xiàng)的線性和二次回歸模型的V- 最優(yōu)設(shè)計(jì)和D- 最優(yōu)設(shè)計(jì)；文獻(xiàn)[18]研究了分層模型中個(gè)體參數(shù)的最優(yōu)設(shè)計(jì)。

本文研究多響應(yīng)線性混合效應(yīng)模型的A- 最優(yōu)設(shè)計(jì)。A- 最優(yōu)設(shè)計(jì)準(zhǔn)則是最常用的設(shè)計(jì)準(zhǔn)則之一，其統(tǒng)計(jì)意義在于使未知參數(shù)估計(jì)量的各分量方差之和最小，因此A- 最優(yōu)設(shè)計(jì)能最大程度地提高參數(shù)估計(jì)的精度[19]。

1 模型和估計(jì)

假設(shè)試驗(yàn)共有n個(gè)被試個(gè)體，對第i個(gè)個(gè)體進(jìn)行mi次觀測。將第i個(gè)個(gè)體的第j次試驗(yàn)點(diǎn)xij處的觀測記為yij，其具有r個(gè)分量，可根據(jù)式(1)計(jì)算得出。

yij=F(xij)βi+εij,i=1, …,n,j=1, …,mi

(1)

式中：試驗(yàn)點(diǎn)xij∈χ，χ為設(shè)計(jì)區(qū)域；F為r×p維回歸函數(shù)矩陣；εij為均值為0、協(xié)方差矩陣為Σ的觀測誤差向量；βi為個(gè)體參數(shù)。βi的均值E(βi)=β；協(xié)方差矩陣Cov(βi)=D，D為半正定矩陣，其秩為q。假設(shè)隨機(jī)效應(yīng)與觀測誤差相互獨(dú)立，不同次試驗(yàn)的觀測誤差之間相互獨(dú)立。

通過分離固定效應(yīng)和隨機(jī)效應(yīng)，可將模型(1)改寫為

yij=F(xij)βi+F(xij)bi+εij,
i=1, …,n,j=1, …,mi

(2)

式中：bi=βi-β是個(gè)體參數(shù)與總體均值的差值。在模型(1)和(2)中，響應(yīng)向量yij的協(xié)方差矩陣為Cov(yij)=F(xij)DFT(xij)+Σ。同一個(gè)體的不同觀測之間是相關(guān)的，其協(xié)方差陣為Cov(yij,yik)=F(xij)DFT(xik),j≠k。不同個(gè)體之間的觀測是相互獨(dú)立的。將第i個(gè)個(gè)體的所有觀測記為Yi=(yTi1,yTi2, …,yTimi)T，其矩陣形式為

Yi=Fiβ+Fibi+εi

Y=Xβ+Zb+ε

式中：N=m1+m2+…+mn為總觀測次數(shù)。

式中：R=Cov(ε)=ΙN?Σ；G=Cov(b)=Ιn?D。

2 最優(yōu)設(shè)計(jì)

為準(zhǔn)確預(yù)測個(gè)體參數(shù)βi，本文期望通過試驗(yàn)設(shè)計(jì)使MSE矩陣在某種程度上達(dá)到最小。在考慮平衡設(shè)計(jì)的基礎(chǔ)上，假設(shè)試驗(yàn)設(shè)計(jì)有k個(gè)不同的試驗(yàn)點(diǎn)x1,x2, …,xk，重復(fù)數(shù)分別為n1,n2, …,nk，可將其記為

進(jìn)一步考慮Kiefer[22]定義的近似設(shè)計(jì)，即只要求

對任一近似設(shè)計(jì)ξ，其標(biāo)準(zhǔn)化信息矩陣定義為

式中：Δ=mD。當(dāng)D非奇異時(shí)，可將MSE矩陣簡化為

由于預(yù)測的精度可以通過MSE矩陣來刻畫，此處借鑒A- 最優(yōu)設(shè)計(jì)的思想，定義準(zhǔn)則函數(shù)為

(n-1)tr{Δ-Δ(M-1(ξ)+Δ)-1Δ})

令Ξ表示在χ上具有非奇異信息矩陣的所有近似設(shè)計(jì)的集合。

定義1如果設(shè)計(jì)ξ*，使得

成立，則稱ξ*為個(gè)體參數(shù)的A-最優(yōu)設(shè)計(jì)。

下面的等價(jià)性定理給出了設(shè)計(jì)為A-最優(yōu)的充要條件。

定理1記N(ξ,Δ)=Δ(M-1(ξ)+Δ)-1。令

φ(x,ξ)=tr{M-1(ξ)FT(x)Σ-1F(x)M-1(ξ)}+

(n-1)tr{N(ξ,Δ)M-1(ξ)FT(x)Σ-1F(x)M-1(ξ)NT(ξ,Δ)}

當(dāng)且僅當(dāng)

則設(shè)計(jì)ξ*是個(gè)體參數(shù)的A-最優(yōu)設(shè)計(jì)，并且在設(shè)計(jì)點(diǎn)ξ*處達(dá)到上確界。

以及

(n-1)tr{Δ(M-1(ξ)+Δ)-1M-1(ξ)(M(ξ)-

(n-1)tr{N(ξ,Δ)M-1(ξ)NT(ξ,Δ)}-

推論1當(dāng)D非奇異時(shí)，令

φ(x,ξ)=tr{M-1(ξ)FT(x)Σ-1F(x)M-1(ξ)}+

(n-1)tr{(M(ξ)+Δ-1)-1(FT(x)Σ-1F(x)+
Δ-1)(M(ξ)+Δ-1)-1}

(3)

當(dāng)且僅當(dāng)

則設(shè)計(jì)ξ*是個(gè)體參數(shù)的A-最優(yōu)設(shè)計(jì)，并且在設(shè)計(jì)點(diǎn)ξ*處達(dá)到上確界。

例1帶有隨機(jī)截距的線性回歸模型求解A-最優(yōu)設(shè)計(jì)

Kubokawa等[24]利用線性混合效應(yīng)模型分析神奈川縣各地至附近車站的距離和附近車站到東京車站距離對該區(qū)域地價(jià)影響。神奈川縣的許多居民常乘火車往返東京，所以居住地到附近車站和附近車站到東京車站的距離會(huì)影響地價(jià)，Kubokawa等[24]運(yùn)用如式(4)所示的模型對1996和2001年的地價(jià)數(shù)據(jù)進(jìn)行分析。

y1ij=β10+β11x1ij+β12x2ij+b1i+ε1ij

y2ij=β20+β21x1ij+β22x2ij+b2i+ε2ij

(4)

式中：y1ij和y2ij分別代表1996和2001年第i個(gè)小區(qū)域、第j個(gè)地塊的地價(jià)，i=1, …,n,j=1, …,m, (x1ij,x2ij)∈[-1, 1]2；x1ij表示該地塊到附近重要車站的(標(biāo)準(zhǔn)化)距離，x2ij表示該地塊附近車站到東京站的(標(biāo)準(zhǔn)化)距離。隨機(jī)效應(yīng)bi=(b1i,b2i)T的協(xié)方差陣為Σb，誤差向量εij=(ε1ij,ε2ij)T的協(xié)方差陣為Σ，并且bi和εij相互獨(dú)立。本例考慮模型(4)的最優(yōu)設(shè)計(jì)問題，利用定理1證明設(shè)計(jì)

是隨機(jī)效應(yīng)的A-最優(yōu)設(shè)計(jì)。

記H≡(Σ+mΣb)-1，由計(jì)算可得

注：例1給出的A- 最優(yōu)設(shè)計(jì)既與隨機(jī)效應(yīng)的方差協(xié)方差矩陣無關(guān)，也與誤差的方差協(xié)方差矩陣無關(guān)，這一點(diǎn)對其他模型不一定成立。