異方差分層線性模型的A-最優(yōu)設(shè)計(jì)

陳佳穎，劉 欣

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

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

隨著分層線性模型在經(jīng)濟(jì)學(xué)、藥物動(dòng)力學(xué)、心理學(xué)等方面的廣泛應(yīng)用，關(guān)于分層線性模型的最優(yōu)試驗(yàn)設(shè)計(jì)問(wèn)題得到了持續(xù)關(guān)注。在觀測(cè)誤差為同方差的場(chǎng)合已有不少研究結(jié)果。在群體參數(shù)方面：Fedorov等[3]建立了刻畫(huà)D-最優(yōu)設(shè)計(jì)的等價(jià)性定理，隨后Schmelter[4-5]、Debusho等[6-7]在這方面做了一系列研究。在個(gè)體參數(shù)方面：Prus等[8]研究了D-最優(yōu)設(shè)計(jì)和線性最優(yōu)設(shè)計(jì)，隨后Prus[9-10]討論了G-最優(yōu)設(shè)計(jì)和極小極大設(shè)計(jì)，He等[11]考慮了R-最優(yōu)設(shè)計(jì)。

實(shí)際中，觀測(cè)誤差為異方差的情況經(jīng)常出現(xiàn)，針對(duì)異方差模型的分析也是統(tǒng)計(jì)研究的重點(diǎn)之一。對(duì)于異方差分層線性模型的研究工作目前較少，Cheng等[12]研究了群體參數(shù)的D-、G-、A-、I-最優(yōu)設(shè)計(jì)問(wèn)題，獲得了隨機(jī)系數(shù)一次回歸模型的最優(yōu)設(shè)計(jì)。Liu等[13]考慮了個(gè)體參數(shù)的D-最優(yōu)設(shè)計(jì)問(wèn)題，給出了D-最優(yōu)設(shè)計(jì)的充分必要條件。

本文研究分層線性模型中的最優(yōu)設(shè)計(jì)問(wèn)題，與大多數(shù)現(xiàn)有文獻(xiàn)不同，本文考慮的是異方差模型。討論個(gè)體參數(shù)的最優(yōu)設(shè)計(jì)，并選擇研究目前結(jié)論尚少的A-最優(yōu)設(shè)計(jì)。同時(shí)，通過(guò)計(jì)算效率將最優(yōu)設(shè)計(jì)和等權(quán)重設(shè)計(jì)進(jìn)行了比較。

1 異方差分層線性模型和預(yù)備知識(shí)

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

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

i=1,2,…,n

(1)

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

(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)不同個(gè)體的觀測(cè)值、個(gè)體參數(shù)向量和觀測(cè)誤差都不相關(guān)。

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

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

(2)

(3)

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

2 個(gè)體參數(shù)的A-最優(yōu)設(shè)計(jì)

本節(jié)以個(gè)體參數(shù)的精準(zhǔn)預(yù)測(cè)為目標(biāo)討論最優(yōu)設(shè)計(jì)問(wèn)題。

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

(4)

式中：支撐點(diǎn)x1,x2,…,xk給出觀測(cè)的條件，而權(quán)重w1,w2,…,wk表示在這些點(diǎn)上進(jìn)行的總觀測(cè)次數(shù)的相對(duì)比例。對(duì)于近似設(shè)計(jì)ξ，其在模型(1)下的標(biāo)準(zhǔn)化信息矩陣可定義為

(5)

相應(yīng)地，對(duì)于近似設(shè)計(jì)ξ，式(3)中定義的均方誤差矩陣可表示為

(6)

式中：Δ=mD。特別地，當(dāng)D是非奇異的時(shí)，式(6)可簡(jiǎn)化為

(7)

(8)

當(dāng)D是非奇異時(shí)，A-準(zhǔn)則函數(shù)簡(jiǎn)化為

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

(9)

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

(10)

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

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

定理1對(duì)任一設(shè)計(jì)ξ，定義其敏感函數(shù)為

(11)

則設(shè)計(jì)ξ*是個(gè)體參數(shù)的A-最優(yōu)設(shè)計(jì)當(dāng)且僅當(dāng)

(12)

此外，上確界在ξ*的支撐點(diǎn)處達(dá)到。

且

因此，A準(zhǔn)則函數(shù)φpred,A(ξ)是凸的。

令δx是設(shè)計(jì)點(diǎn)為x的單點(diǎn)設(shè)計(jì)，根據(jù)文獻(xiàn)[15]可知：

進(jìn)而得到，設(shè)計(jì)ξ*∈Ξ是個(gè)體參數(shù)的A-最優(yōu)設(shè)計(jì)

即ξ*∈Ξ是個(gè)體參數(shù)的A-最優(yōu)設(shè)計(jì)當(dāng)且僅當(dāng)

推論1設(shè)D是非奇異的，定義設(shè)計(jì)ξ的敏感函數(shù)為

(13)

設(shè)計(jì)ξ*是個(gè)體參數(shù)的A-最優(yōu)設(shè)計(jì)當(dāng)且僅當(dāng)

(14)

此外，上確界在ξ*的支撐點(diǎn)處達(dá)到。

例1隨機(jī)截距一次回歸模型的A-最優(yōu)設(shè)計(jì)。

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

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

(15)

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

Chang[16]在研究加權(quán)多項(xiàng)式回歸模型的D-最優(yōu)設(shè)計(jì)問(wèn)題時(shí)考慮了如式(16)所示的方差權(quán)重函數(shù)。

(16)

在這一權(quán)重函數(shù)下，Cheng等[12]求出了隨機(jī)斜率一次回歸模型的群體參數(shù)的D-最優(yōu)設(shè)計(jì)。

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

(17)

其中w∈[0,1]為待定的權(quán)重參數(shù)。形如式(17)的設(shè)計(jì)所對(duì)應(yīng)的A-最優(yōu)準(zhǔn)則函數(shù)為

對(duì)準(zhǔn)則函數(shù)求導(dǎo)并令其等于0，得到方程：

(18)

圖1 等權(quán)重設(shè)計(jì)ξeff相對(duì)于A-最優(yōu)設(shè)計(jì)的效率圖Fig.1 Efficiency of equireplicated design ξeff relative to A-optimal design

例2隨機(jī)斜率一次回歸模型的A-最優(yōu)設(shè)計(jì)。

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

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

(19)

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

Cheng等[12]也給出了權(quán)重函數(shù)為

λ2(x)=x2+1

(20)

時(shí)，模型(19)的群體參數(shù)的幾類最優(yōu)設(shè)計(jì)。

現(xiàn)在來(lái)求其個(gè)體參數(shù)的A-最優(yōu)設(shè)計(jì)，求解方法同例1，只需考慮形如式(17)的兩點(diǎn)設(shè)計(jì)。

此時(shí)A-準(zhǔn)則函數(shù)為

圖2 隨機(jī)斜率一次回歸模型的A-最優(yōu)設(shè)計(jì)的權(quán)重函數(shù)(異方差結(jié)構(gòu)為λ2(x))Fig.2 Weight function of A-optimal design for the straight line regression with random slope (heteroskedasticity structure is λ2(x))

圖3 等權(quán)重設(shè)計(jì)ξeff相對(duì)于A-最優(yōu)設(shè)計(jì)的效率圖Fig.3 Efficiency of equireplicated design ξeff relative to A-optimal design