異方差分層線性模型的A-最優設計

陳佳穎，劉 欣

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

分層線性模型也稱混合效應模型或隨機系數模型，是同時考慮觀測對象的群體效應和個體效應的兩階段模型，通常用于對多個個體重復測量的數據進行建模。其理論分析日臻完善，主要有Vonesh等[1]和Verbeke等[2]的研究。

隨著分層線性模型在經濟學、藥物動力學、心理學等方面的廣泛應用，關于分層線性模型的最優試驗設計問題得到了持續關注。在觀測誤差為同方差的場合已有不少研究結果。在群體參數方面：Fedorov等[3]建立了刻畫D-最優設計的等價性定理，隨后Schmelter[4-5]、Debusho等[6-7]在這方面做了一系列研究。在個體參數方面：Prus等[8]研究了D-最優設計和線性最優設計，隨后Prus[9-10]討論了G-最優設計和極小極大設計，He等[11]考慮了R-最優設計。

實際中，觀測誤差為異方差的情況經常出現，針對異方差模型的分析也是統計研究的重點之一。對于異方差分層線性模型的研究工作目前較少，Cheng等[12]研究了群體參數的D-、G-、A-、I-最優設計問題，獲得了隨機系數一次回歸模型的最優設計。Liu等[13]考慮了個體參數的D-最優設計問題，給出了D-最優設計的充分必要條件。

本文研究分層線性模型中的最優設計問題，與大多數現有文獻不同，本文考慮的是異方差模型。討論個體參數的最優設計，并選擇研究目前結論尚少的A-最優設計。同時，通過計算效率將最優設計和等權重設計進行了比較。

1 異方差分層線性模型和預備知識

假設在設計區域χ上，對n個個體進行重復觀測，其中對個體i觀測mi次，并用yij表示第i個個體的第j次觀測值。本文考慮的異方差分層線性模型如下：

yij=fT(xij)βi+ε(xij)，j=1,2,…,mi,

i=1,2,…,n

(1)

式中：xij∈χ表示試驗條件；f=(f1,…,fp)T表示已知的回歸函數向量；隨機向量βi=(βi1,…,βip)T為個體i的個體參數；ε(xij)為隨機誤差。此外，模型滿足：

(1)E(βi)=β=(β1,…,βp)T，Cov(βi)=σ2D,i=1,2,…,n；

(2)E(ε(xij))=0,var(ε(xij))=σ2/λ(xij),j=1,2,…,mi,i=1,2,…,n；

(3)不同個體的觀測值、個體參數向量和觀測誤差都不相關。

其中，總體平均值β未知，σ2已知，λ(x)是定義在設計區域χ上的正實值函數。

當λ(x)恒等于1時，模型(1)就是Prus等[8]研究的同方差分層線性模型。

(2)

(3)

式中：In為n×n單位矩陣；Jn為所有元素等于1的n×n矩陣；?表示矩陣的Kronecker積。

2 個體參數的A-最優設計

本節以個體參數的精準預測為目標討論最優設計問題。

參照Prus等[8]的做法，本文考慮Kiefer[14]定義的近似設計為

(4)

式中：支撐點x1,x2,…,xk給出觀測的條件，而權重w1,w2,…,wk表示在這些點上進行的總觀測次數的相對比例。對于近似設計ξ，其在模型(1)下的標準化信息矩陣可定義為

(5)

相應地，對于近似設計ξ，式(3)中定義的均方誤差矩陣可表示為

(6)

式中：Δ=mD。特別地，當D是非奇異的時，式(6)可簡化為

(7)

(8)

當D是非奇異時，A-準則函數簡化為

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

(9)

記Ξ為設計空間中所有具有非奇異信息矩陣的近似設計的全體。如果設計ξ*使得式(8)達到最小，即

(10)

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

下文的等價性定理給出了一個設計為A-最優的充分必要條件。

定理1對任一設計ξ，定義其敏感函數為

(11)

則設計ξ*是個體參數的A-最優設計當且僅當

(12)

此外，上確界在ξ*的支撐點處達到。

且

因此，A準則函數φpred,A(ξ)是凸的。

令δx是設計點為x的單點設計，根據文獻[15]可知：

進而得到，設計ξ*∈Ξ是個體參數的A-最優設計

即ξ*∈Ξ是個體參數的A-最優設計當且僅當

推論1設D是非奇異的，定義設計ξ的敏感函數為

(13)

設計ξ*是個體參數的A-最優設計當且僅當

(14)

此外，上確界在ξ*的支撐點處達到。

例1隨機截距一次回歸模型的A-最優設計。

考慮χ=[0,1]上的隨機截距一次回歸模型：

yij=βi1+β2xj+ε(xj)

(15)

式中：βi1是均值為β1、方差為σ2d1的隨機截距。對此模型，p=2,f(x)=(1,x)T且矩陣D=diag(d1,0)。

Chang[16]在研究加權多項式回歸模型的D-最優設計問題時考慮了如式(16)所示的方差權重函數。

(16)

在這一權重函數下，Cheng等[12]求出了隨機斜率一次回歸模型的群體參數的D-最優設計。

下面在此權重函數下求模型(15)的個體參數的A-最優設計。首先，根據Cheng等[12]研究中的引理3.1可以確定A-最優設計有兩個支撐點，分別是x1=0,x2=1，即A-最優設計具有如式(17)所示的形式。

(17)

其中w∈[0,1]為待定的權重參數。形如式(17)的設計所對應的A-最優準則函數為

對準則函數求導并令其等于0，得到方程：

(18)

圖1 等權重設計ξeff相對于A-最優設計的效率圖Fig.1 Efficiency of equireplicated design ξeff relative to A-optimal design

例2隨機斜率一次回歸模型的A-最優設計。

考慮χ=[0,1]上的隨機斜率一次回歸模型：

yij=β1+βi2xj+ε(xj)

(19)

式中:只有斜率βi2是隨機的，其均值為β2，方差為σ2d2，因此矩陣D=diag(0,d2)。

Cheng等[12]也給出了權重函數為

λ2(x)=x2+1

(20)

時，模型(19)的群體參數的幾類最優設計。

現在來求其個體參數的A-最優設計，求解方法同例1，只需考慮形如式(17)的兩點設計。

此時A-準則函數為

圖2 隨機斜率一次回歸模型的A-最優設計的權重函數(異方差結構為λ2(x))Fig.2 Weight function of A-optimal design for the straight line regression with random slope (heteroskedasticity structure is λ2(x))

圖3 等權重設計ξeff相對于A-最優設計的效率圖Fig.3 Efficiency of equireplicated design ξeff relative to A-optimal design