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

段冬艷， 劉 欣

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

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

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

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

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

1 模型和估計

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

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

(1)

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

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

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

(2)

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

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

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

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

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

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

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

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

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

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

下面的等價性定理給出了設(shè)計為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è)計ξ*是個體參數(shù)的A-最優(yōu)設(shè)計，并且在設(shè)計點ξ*處達(dá)到上確界。

以及

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

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

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

φ(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è)計ξ*是個體參數(shù)的A-最優(yōu)設(shè)計，并且在設(shè)計點ξ*處達(dá)到上確界。

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

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

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

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

(4)

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

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

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

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