一類線性混合模型的廣義估計(jì)方程方法應(yīng)用

余俊飛,趙慧秀

(南京理工大學(xué) 理學(xué)院,江蘇 南京 210094)

0 引言

含有兩個(gè)方差分量的線性混合模型是常用的縱向數(shù)據(jù)統(tǒng)計(jì)分析模型，它在線性模型中引入隨機(jī)效應(yīng)，建模了試驗(yàn)個(gè)體的不均勻性，并且刻畫了同一個(gè)體內(nèi)部觀測(cè)值之間的相關(guān)性.該模型在生物學(xué)、醫(yī)學(xué)、經(jīng)濟(jì)學(xué)等領(lǐng)域具有廣泛的應(yīng)用，因此對(duì)于該類模型參數(shù)估計(jì)的研究具有重要意義.

線性混合模型又稱為隨機(jī)效應(yīng)模型，已有文獻(xiàn)中有許多關(guān)于模型參數(shù)估計(jì)的研究.Liard and Ware[1]提出兩階段隨機(jī)效應(yīng)模型，討論了結(jié)合經(jīng)驗(yàn)貝葉斯和極大似然估計(jì)的重復(fù)測(cè)量模型擬合方法.Dempster and Selwyn[2]闡述了含有兩個(gè)方差分量線性混合模型的統(tǒng)計(jì)計(jì)算方法，并比較了基于限制極大似然估計(jì)的EM算法和牛頓迭代法的效率.王松桂[3]概述了線性混合模型在參數(shù)估計(jì)、假設(shè)檢驗(yàn)等方面的重要成果，并且提出了模型參數(shù)的一種譜分解估計(jì).吳密霞等[4]在一組簡(jiǎn)單條件下證明了線性混合模型固定效應(yīng)和方差分量可以同時(shí)達(dá)到最優(yōu)估計(jì).Gumedze and Dunne[5]綜述了線性混合模型不同組成部分的參數(shù)估計(jì)和推斷，重點(diǎn)討論了方差參數(shù)的估計(jì)以及隨機(jī)效應(yīng)的推理過程.以上方法多是基于數(shù)據(jù)分布的假設(shè)，對(duì)于縱向數(shù)據(jù)的統(tǒng)計(jì)分析，Liang and Zeger[6-7]提出的廣義估計(jì)方程具有廣泛的應(yīng)用.該方法不假設(shè)響應(yīng)變量的分布，通過給定的工作相關(guān)陣替代真實(shí)的個(gè)體內(nèi)部相關(guān)矩陣，其中工作相關(guān)陣的選擇不影響回歸參數(shù)及其方差估計(jì)的相合性[8].

本文結(jié)合廣義估計(jì)方程處理一類線性混合模型的參數(shù)估計(jì)問題，該類模型中含有兩個(gè)方差分量，其中隨機(jī)效應(yīng)的方差用來刻畫個(gè)體內(nèi)部的相關(guān)性.在此基礎(chǔ)上可以計(jì)算出方差分量的矩估計(jì)，并應(yīng)用廣義估計(jì)方程估計(jì)模型參數(shù).文中以一組大鼠繁殖數(shù)據(jù)進(jìn)行了實(shí)例分析.

1 基本模型和參數(shù)估計(jì)

1.1 線性混合模型

考慮縱向數(shù)據(jù)下的線性混和模型，假設(shè)有m個(gè)被試個(gè)體,xij和yij分別是第i(i=1,2,…，m)個(gè)個(gè)體在第j(j=1,2,…，ni)次觀測(cè)的協(xié)變量和響應(yīng)變量.Yi是ni×1的個(gè)體響應(yīng)向量，Xi為已知的ni×p協(xié)變量矩陣，Zi是已知的ni×q維設(shè)計(jì)矩陣，含有兩個(gè)方差分量的線性混合模型表示為

Yi=Xiβ+Zibi+εi，

(1)

E(Yi)=Xiβ，

(2)

則模型中個(gè)體的響應(yīng)變量服從多元正態(tài)分布.這樣的模型在縱向數(shù)據(jù)分析中經(jīng)常遇到，通過最大化對(duì)數(shù)似然函數(shù)可以求出模型的極大似然估計(jì).

1.2 廣義估計(jì)方程

var(yij)=φv(μij)，

其中：φ是散布參數(shù)；v是已知的方差函數(shù).

基于以上假設(shè)，Weddueburn[9]提出了如下的擬似然估計(jì)方程

(3)

其中：

對(duì)于個(gè)體內(nèi)部具有相關(guān)性的縱向數(shù)據(jù)，Vi的結(jié)構(gòu)是非對(duì)角矩陣.Liang and Zeger在此基礎(chǔ)上提出了廣義估計(jì)方程方法.

定義一個(gè)對(duì)角矩陣

(4)

(5)

方程(5)不同工作相關(guān)陣下相關(guān)系數(shù)的估計(jì)方法相同.根據(jù)α估計(jì)的不同方式，可以迭代求出相應(yīng)的參數(shù)估計(jì)值.

則

(6)

1.3 模型參數(shù)估計(jì)

含有兩個(gè)方差分量線性混合模型的邊際均值和協(xié)方差矩陣計(jì)算公式為

E(Yi)=E[E(Yi|bi)]，

Vi=cov[E(Yi|bi)]+E[cov(Yi|bi)],

(7)

由此可以得到式(2)的相應(yīng)結(jié)果.

首先對(duì)隨機(jī)截距模型進(jìn)行分析，即模型中只含有單個(gè)隨機(jī)效應(yīng)，其協(xié)方差矩陣為

Jni是元素全為1的ni×ni矩陣.從Vi的結(jié)構(gòu)可以看出該模型刻畫了個(gè)體內(nèi)部觀測(cè)值之間的等相關(guān)性，這與廣義估計(jì)方程中假設(shè)的可交換相關(guān)矩陣是等價(jià)的.在廣義估計(jì)方程框架下可以采用矩估計(jì)和擬加權(quán)最小二乘等方法估計(jì)等相關(guān)矩陣下的相關(guān)系數(shù)，并通過對(duì)應(yīng)關(guān)系求出模型的方差分量.令

N=n1+n2+…+nm，rij=

(Yij-μij)/v(μij)1/2，

則α的矩估計(jì)為

(8)

根據(jù)最小二乘中估計(jì)方差的思想，可以采用如下公式估計(jì)散布參數(shù)φ.

(9)

根據(jù)方程(5)(8)(9)進(jìn)行迭代,可以求出回歸參數(shù)β、相關(guān)系數(shù)α和散布參數(shù)φ.顯然β是對(duì)應(yīng)模型(1)中固定效應(yīng)參數(shù)的估計(jì)，而方差分量估計(jì)的對(duì)應(yīng)結(jié)果為

(10)

對(duì)含有q個(gè)隨機(jī)效應(yīng)參數(shù)的隨機(jī)系數(shù)模型進(jìn)行分析，這時(shí)模型協(xié)方差矩陣的形式為

令

(11)

(12)

然后通過方程(3)、(11)、(12)進(jìn)行迭代，即可求出模型未知參數(shù)估計(jì)值.

2 實(shí)例分析

本文研究的大鼠繁殖數(shù)據(jù)可在文獻(xiàn)[2]中獲得，主要目的是評(píng)估試驗(yàn)劑量對(duì)大鼠繁殖能力的影響.試驗(yàn)中30只母鼠被隨機(jī)分配到3個(gè)不同的處理組，分別是對(duì)照組、高劑量組、低劑量組，每組均有10只母鼠.母鼠生產(chǎn)后記錄每只幼鼠的體重，最后得到27窩幼鼠數(shù)據(jù)，并以幼鼠的體重作為評(píng)估試驗(yàn)效應(yīng)的指標(biāo).

每只母鼠的情況具有一定的差異，一般認(rèn)為來自同一窩的幼鼠之間存在等相關(guān)性，運(yùn)用含有兩個(gè)方差分量的線性混合模型來分析該數(shù)據(jù)是合適的.模型的響應(yīng)變量是幼鼠的體重(weight)，相應(yīng)的協(xié)變量為所在窩中母鼠是否在高劑量組(high)、母鼠是否在低劑量組(low)、幼鼠性別(sex)、所在窩中幼鼠只數(shù)(high).其中幼鼠平均體重和所在窩中幼鼠只數(shù)的關(guān)系如圖1所示.

圖1 幼鼠平均體重與幼鼠只數(shù)的散點(diǎn)圖

基于以上數(shù)據(jù)特征建立如下線性混合模型

(13)

應(yīng)用廣義估計(jì)方程方法得到模型參數(shù)，結(jié)果如表1所列.為評(píng)估該方法的估計(jì)效果，表中包含了極大似然估計(jì)的結(jié)果用于對(duì)比分析.

從表1結(jié)果來看，兩種方法的參數(shù)估計(jì)結(jié)果大體上相同，結(jié)果表明藥物劑量、幼鼠性別、每窩幼鼠只數(shù)對(duì)幼鼠的體重有明顯影響，這些情況與數(shù)據(jù)的直觀分析結(jié)果是一致的.表1結(jié)果顯示廣義估計(jì)方程方法回歸參數(shù)估計(jì)的標(biāo)準(zhǔn)差更小，該方法不假設(shè)響應(yīng)變量的聯(lián)合分布，因此能夠得到穩(wěn)健的估計(jì)結(jié)果.

隨機(jī)模擬能夠評(píng)估方法的有效性，以表1中應(yīng)用極大似然估計(jì)得到的結(jié)果為真值，即

表1 實(shí)例數(shù)據(jù)參數(shù)估計(jì)結(jié)果

應(yīng)用統(tǒng)計(jì)軟件依據(jù)公式(13)生成模擬數(shù)據(jù)，首先隨機(jī)生成一組數(shù)據(jù)，應(yīng)用兩種方法估計(jì)相關(guān)參數(shù).而后生成100組數(shù)據(jù)，計(jì)算各參數(shù)估計(jì)的平均值，其結(jié)果如表2所列.

表2的模擬數(shù)據(jù)分析結(jié)果表明兩種方法固定效應(yīng)的參數(shù)估計(jì)與真值差異不大.廣義估計(jì)方程固定效應(yīng)估計(jì)的標(biāo)準(zhǔn)差較小，對(duì)于方差分量，在已知分布情形下，極大似然估計(jì)更接近真值.

表2 模擬數(shù)據(jù)參數(shù)估計(jì)結(jié)果

3 結(jié)語

文中運(yùn)用隨機(jī)截距模型擬合了大鼠繁殖數(shù)據(jù)，采用廣義估計(jì)方程方法得到的參數(shù)估計(jì)取得了較好的效果，這說明提出的方法是可行的.實(shí)際應(yīng)用中，隨機(jī)系數(shù)模型常用來建模時(shí)間相關(guān)的縱向數(shù)據(jù)，該模型同樣可以應(yīng)用廣義估計(jì)方程方法進(jìn)行分析.