多響應線性混合效應模型的A- 最優設計

段冬艷， 劉 欣

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

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

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

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

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

1 模型和估計

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

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

(1)

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

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

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

(2)

式中：bi=βi-β是個體參數與總體均值的差值。在模型(1)和(2)中，響應向量yij的協方差矩陣為Cov(yij)=F(xij)DFT(xij)+Σ。同一個體的不同觀測之間是相關的，其協方差陣為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為總觀測次數。

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

2 最優設計

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

進一步考慮Kiefer[22]定義的近似設計，即只要求

對任一近似設計ξ，其標準化信息矩陣定義為

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

由于預測的精度可以通過MSE矩陣來刻畫，此處借鑒A- 最優設計的思想，定義準則函數為

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

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

定義1如果設計ξ*，使得

成立，則稱ξ*為個體參數的A-最優設計。

下面的等價性定理給出了設計為A-最優的充要條件。

定理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(ξ,Δ)}

當且僅當

則設計ξ*是個體參數的A-最優設計，并且在設計點ξ*處達到上確界。

以及

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

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

推論1當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)

當且僅當

則設計ξ*是個體參數的A-最優設計，并且在設計點ξ*處達到上確界。

例1帶有隨機截距的線性回歸模型求解A-最優設計

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

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

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

(4)

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

是隨機效應的A-最優設計。

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

注：例1給出的A- 最優設計既與隨機效應的方差協方差矩陣無關，也與誤差的方差協方差矩陣無關，這一點對其他模型不一定成立。