關(guān)于分層線性模型樣本容量問題的研究

關(guān)于分層線性模型樣本容量問題的研究

1研究背景很多社會研究都涉及分層數(shù)據(jù)結(jié)構(gòu)，例如，經(jīng)濟(jì)學(xué)家探求在多個國家中經(jīng)濟(jì)政策是如何影響居民的消費(fèi)行為，研究采集的觀測數(shù)據(jù)不僅包括以國家為層次的經(jīng)濟(jì)指標(biāo)，還包括以家庭為單位的信息，因此整個觀測的數(shù)據(jù)結(jié)構(gòu)是分層的。此時，同屬一個層次的個體之間的相關(guān)性會大于來自不同層次的個體之間的相關(guān)性，整個觀測樣本就不再具有獨立同分布性質(zhì)，如果繼續(xù)使用經(jīng)典的線性回歸模型，就會得到有偏的參數(shù)估計和錯誤的統(tǒng)計推斷結(jié)果。近年來，隨著分層線性模型統(tǒng)計理論的發(fā)展，一套完整的應(yīng)用于分層結(jié)構(gòu)數(shù)據(jù)的統(tǒng)計推斷方法已經(jīng)建立起來，并且能得到有效的參數(shù)估計。分層線性模型(hierarchicallinearmodels)的稱謂最早由Lindley和Smith(1972)[1]提出。這個模型在不同的研究領(lǐng)域有不同的稱呼，在社會學(xué)研究中，它經(jīng)常被稱為多層線性模型(multilevellinearmodel)；在生物統(tǒng)計研究中常用的名字是混合效應(yīng)模型(mixed-effectsmodels)和隨機(jī)效應(yīng)模型(random-effectsmodels)；計量經(jīng)濟(jì)學(xué)文獻(xiàn)稱之為隨機(jī)系數(shù)回歸模型(random-coefficientregressionmodels)等。模型的一般形式為：目前比較常用的分層線性模型的估計方法有極大似然法(ML)，約束極大似然法(REML)(Littell，1996)[2]和完全貝葉斯法(FullBayes)。理論研究表明，大樣本情況下ML和REML得到的參數(shù)估計都是一致最優(yōu)估計量，但是在樣本較小的情況下，REML在估計方差協(xié)方差成分和T時，考慮了固定效應(yīng)系數(shù)γ的不確定性，因此REML比ML得到了更可靠的和T的估計值。FullBayes考慮了所有輔助參數(shù)的不確定性，因此理論上得到的參數(shù)估計較ML和REML更可靠，但是此估計方法需要運(yùn)用MCMC算法，很多情況下，并非能得到一個收斂的Markov鏈，當(dāng)分層線性模型形式較復(fù)雜時（如待定參數(shù)向量增多、層數(shù)增加等），F(xiàn)ullBayes方法相當(dāng)復(fù)雜。因此很多統(tǒng)計軟件采用ML和REML估計分層線性模型，本研究采用REML方法估計。由于ML和REML方法估計的前提假設(shè)都是樣本量要足夠大，因此樣本量較小時，這兩種估計都是有偏的，由此得到的參數(shù)置信區(qū)間和假設(shè)檢驗都是不可信的。因此，分層線性模型樣本量問題的研究是一個重大課題，近十幾年來，很多的學(xué)者都致力于這方面的研究。Bryk和Raudenbush(1992)[3]很早就指出，可以憑借OLS回歸的經(jīng)驗法則：增加一個解釋變量至少需要增加10個觀測樣本，聯(lián)系到分層線性模型的估計，增加一個層2的結(jié)果變量（層1模型中待定的隨機(jī)參數(shù)）至少需要增加10個觀測樣本，這個準(zhǔn)則只是OLS回歸經(jīng)驗法則的平移，并沒有清楚回答模型要得到較好估計所需的最小樣本量問題。Kim(1990)[4]在研究分層線性模型斜率參數(shù)估計的時候，發(fā)現(xiàn)當(dāng)組數(shù)較少而組內(nèi)的觀測值相對很大時，固定效應(yīng)和隨機(jī)效應(yīng)的估計會產(chǎn)生很大的偏差，但是Kim的研究在相同的條件下，只進(jìn)行了50次模擬計算，因此結(jié)果不太可靠。Mok(1995)[5]的研究也得出了類似Kim的結(jié)論，并且他進(jìn)一步指出參數(shù)估計的偏差和效率更大程度取決于樣本中的組數(shù)大小。Busing(1993)[6]指出當(dāng)組數(shù)超過300時，隨機(jī)效應(yīng)的估計才是無偏的。Kreft(1996)[7]運(yùn)用模擬技術(shù)探討了分層線性模型的勢，建議30/30準(zhǔn)則，即30個組，每組30個觀測值可以得到較可靠的估計。Hox(1998)[8]在Keft(1996)之后更詳盡闡述表明，50組，每組20個觀測值可以使交互效應(yīng)得到較可靠的假設(shè)檢驗；100組，每組10個觀測值可以較可靠進(jìn)行隨機(jī)效應(yīng)的假設(shè)檢驗。另外Bliese(1998)[9]明確指出組內(nèi)相關(guān)系數(shù)(ICC)與樣本量相互聯(lián)系。Snijders和Bosker(1999)[10]闡明對于較大的ICC值，增加組數(shù)也不會得到一致的標(biāo)準(zhǔn)差的估計。Mass和Hox(2005)[11]針對不同的層1和層2樣本量進(jìn)行了模擬研究，結(jié)果表明層2的樣本數(shù)小于50，會導(dǎo)致有偏的協(xié)方差成分估計?？v觀以上研究，其方法主要是運(yùn)用數(shù)據(jù)模擬的方法，針對不同層1、層2樣本量的組合，估計分層線性模型，比較固定效應(yīng)、隨機(jī)效應(yīng)和方差協(xié)方差成分的估計值。在計算機(jī)發(fā)展日新月異的今天，進(jìn)行模擬計算簡單易行，應(yīng)用廣泛。但是對于此類問題，以往研究存在許多問題：①大部分研究考慮的是參數(shù)估計的準(zhǔn)確性，只有很少的一部分重視估計的標(biāo)準(zhǔn)差的準(zhǔn)確性。②許多研究在討論估計參數(shù)假設(shè)檢驗或構(gòu)造置信區(qū)間時，想當(dāng)然的運(yùn)用了大樣本的近似理論，即參數(shù)估計近似服從正態(tài)分布，但是實際上，在樣本量較小時，分層線性模型的參數(shù)估計值，尤其是方差協(xié)方差成分的估計值并非服從正態(tài)分布(Raudenbush(1984))。2研究方法本研究仍采用數(shù)據(jù)模擬的方法，用R語言編譯計算程序。另外，研究采用簡單的兩層模型，每層模型各有一個解釋變量：在數(shù)據(jù)模擬過程中，研究分別改變組數(shù)，組內(nèi)個體數(shù)和組間相關(guān)系數(shù)ICC，在樣本變化的條件下，構(gòu)造參數(shù)估計的置信區(qū)間，計算置信區(qū)間覆蓋真實參數(shù)的覆蓋率來考查估計值的分布情況。①組內(nèi)個體數(shù)取值為5，30，50，ICC取值為1，2，3，考察組數(shù)從5至100改變情況下，估計值的分布如何變化；②組數(shù)取值為30，50，100，ICC取值為1，2，3，考察組內(nèi)個體數(shù)從5至100改變情況下，估計值的分布如何變化。在構(gòu)造置信區(qū)間時，以往研究直接將方差用近似標(biāo)準(zhǔn)差代替，構(gòu)造正態(tài)分布的置信區(qū)間，而本研究用另外一種方法計算參數(shù)估計的標(biāo)準(zhǔn)差。R程序lme4程序包中的lmer命令可以提供分層線性模型的計算，它提供固定效應(yīng)系數(shù)和方差協(xié)方差成分的估計值，以及固定效應(yīng)系數(shù)估計的標(biāo)準(zhǔn)差，但是不提供方差協(xié)方差成分估計的標(biāo)準(zhǔn)差，lmer命令的編寫者DougBates(2008)表示，第一，方差協(xié)方差成分估計的標(biāo)準(zhǔn)差的計算是非常復(fù)雜的，很多情況下，估計值已經(jīng)在參數(shù)空間以外，即便是估計出來，也是無意義的，因此有些軟件（如SAS）提供標(biāo)準(zhǔn)差也并不可靠。第二，當(dāng)參數(shù)分布差不多是對稱的情況下，提供參數(shù)估計的標(biāo)準(zhǔn)差才是有意義的，因為可以由此構(gòu)造置信區(qū)間。HarveyGoldstein(1999)提到，樣本不是很大的情況下，如果隨機(jī)擾動項不再服從正態(tài)分布，分層線性模型的固定效應(yīng)和隨機(jī)效應(yīng)的估計仍然是一致的，但是其估計的標(biāo)準(zhǔn)差不能用來構(gòu)建置信區(qū)間和進(jìn)行顯著性檢驗。一種替代的方法是運(yùn)用密集計算中的Jackknife和Bootstrap計算估計的標(biāo)準(zhǔn)差，由Jackknife和Bootstrap理論可知估計的標(biāo)準(zhǔn)差可以利用正態(tài)分布構(gòu)建置信區(qū)間和進(jìn)行顯著性檢驗。3研究步驟3.1模擬數(shù)據(jù)的生成(4)討論覆蓋真實值情況：如果以上的置信區(qū)間確實覆蓋了真實值，則取值為1；否則取0。(5)計算覆蓋率：重復(fù)第一步至第四步1000次，統(tǒng)計覆蓋次數(shù)，計算覆蓋率。(6)在不同樣本量情況下，考察覆蓋率隨樣本量如何變化。固定n(=30，50，100)，p(=1，2，3)，N從5到100變化；固定N(=5，30，50)，p(=1，2，3)，n從5到100變化。3.2.2運(yùn)用參數(shù)Bootstrap方法和Jackknife方法相比，參數(shù)Bootstrap方法步驟(2)中有差異，它將(2)分成以下幾步。構(gòu)造置信區(qū)間和覆蓋率的計算和Jackknife的過程一樣。4研究結(jié)果4.1固定效應(yīng)和方差協(xié)方差成分的置信區(qū)間覆蓋情況由于不論是Jackknife還是Bootstrap方法，運(yùn)算量都是非常大的，因此我們先考察在組成樣本三個條件N，n，p取不同值的組合下置信區(qū)間的覆蓋情況。4.1.1Jackknife方法表1的結(jié)果表明，當(dāng)N=5時，計算的覆蓋率都離95%較遠(yuǎn)，這時若n的取值較大，如n=100，覆蓋率會相對提高，但是方差協(xié)方差成分中的覆蓋率還是很低，這說明分層線性模型的估計很大程度上依賴于N的個數(shù)。如果僅僅考慮參數(shù)點估計，當(dāng)N的取值很小，即使n取值很大，的估計值也不能收斂，因此存在某個整體較大的樣本量（較大的n和較小的N），方差協(xié)方差成分的估計是相當(dāng)不可靠的。另外從表1中，還可以看出組內(nèi)相關(guān)系數(shù)(ICC)的改變對置信區(qū)間的覆蓋率沒有太大影響。當(dāng)N達(dá)到30，n達(dá)到30時，固定效應(yīng)和的估計是比較可靠的，其覆蓋率基本上達(dá)到93%，只有三個值小于93%；但是此時對于而言，大部分的覆蓋率仍然偏低，只有個別能達(dá)到93%以上。表1N，n，p不同取值的組合下，Jackknife和Bootstrap方法計算參數(shù)95%置信區(qū)間覆蓋率(%)注：“J”表示Jackknife法得到的覆蓋率；“B”表示Bootstrap法得到的覆蓋率。4.1.2Bootstrap方法和Jackknife相比，Bootstrap方法得出的結(jié)論很相似，但是對于N較小的時候，可以清楚地看到ICC的值越大，其95%的置信區(qū)間的覆蓋率越??；另外對于σ的估計和Jackknife不一樣，當(dāng)N=5時，只要提高組內(nèi)的樣本量n，例如當(dāng)n=100時，其95%的置信區(qū)間的覆蓋率的均值為94.5%；當(dāng)n=30，N=5,95%的置信區(qū)間的覆蓋率的均值也達(dá)到93%，這說明σ的估計依賴的是整個樣本容量，若組數(shù)稍小時，可以通過增加組內(nèi)的個體數(shù)，提高估計的可靠性。從以上的分析可以看出，Jackknife和Bootstrap方法只是在樣本組數(shù)較小時有差異，而且其差異主要表現(xiàn)在σ的估計上，這種差異來源于這兩種方法計算時再抽樣的過程不同，當(dāng)然這種ParametricBootstrap的方法在樣本組數(shù)較小時優(yōu)于Jackknife的方法，但當(dāng)樣本數(shù)增加時區(qū)別不大，而ParametricBootstrap的計算量遠(yuǎn)遠(yuǎn)超過Jackknife的計算，因此本研究下面的分析，當(dāng)n和N連續(xù)的變化時，我們考慮Jackknife的計算方法。4.2N從5到100變化時參數(shù)估計置信區(qū)間的覆蓋率現(xiàn)在考慮n取30，50，100；ICC取1，2，3的不同組合，N從5到100變化。因為ICC的變化對覆蓋率的影響不大，因此以下的覆蓋率均為ICC三種取值下的平均覆蓋率。另外，為了考察覆蓋率的可接受性，我們利用隨機(jī)模擬的方法建立相應(yīng)樣本量下的覆蓋率的置信區(qū)間（利用隨機(jī)數(shù)均值覆蓋率的置信區(qū)間），構(gòu)成了相應(yīng)的置信帶。4.2.1固定效應(yīng)以為例，其他三個固定效應(yīng)的變化情況相似。從圖1清楚地看到，當(dāng)N取值很小時（如N=5，10），固定效應(yīng)的覆蓋率是很低的；當(dāng)N達(dá)到35時，三種情況下的覆蓋率都進(jìn)入了置信帶。4.2.2方差協(xié)方差考慮σ和的變化參照。①σ的估計情況；當(dāng)N很小時（如N=5,10），標(biāo)準(zhǔn)差σ覆蓋率遠(yuǎn)遠(yuǎn)低于置信帶的下界；隨著n的增加，有更多的覆蓋率進(jìn)入置信帶，平均當(dāng)N達(dá)到50時，三種情況下的覆蓋率完全進(jìn)入置信帶。②的估計情況：的覆蓋率在n的三種取值情況下，都比固定效應(yīng)γ和標(biāo)準(zhǔn)差σ相對應(yīng)的覆蓋率低，因此N從5變化至100時，只有較少的部分進(jìn)入置信帶，對于n的三種取值，N達(dá)到70時，覆蓋率幾乎完全進(jìn)入置信帶。對于較小的N，隨著n的增加，參數(shù)覆蓋率更靠近置信帶。因此對于τ的估計，可以通過n的增加得以改善。4.3n從5到100變化時參數(shù)估計置信區(qū)間的覆蓋率現(xiàn)在考慮N取5，30，50；ICC取1，2，3的不同組合，N從5到100變化。因為ICC的變化對覆蓋率的影響不大，因此以下的覆蓋率均為ICC三種取值下的平均覆蓋率值。另外，為了考察覆蓋率的可接受性，我們利用隨機(jī)模擬的方法建立相應(yīng)樣本量下的覆蓋率的置信區(qū)間（利用隨機(jī)數(shù)均值覆蓋率的置信區(qū)間），構(gòu)成了相應(yīng)的置信帶。4.3.1固定效應(yīng)以為例，其他三個固定效應(yīng)的變化情況相似。圖4（見下頁）清楚地表明，N為5時，不論n如何增加，參數(shù)置信區(qū)間的覆蓋率都遠(yuǎn)離置信帶；當(dāng)N取30時，只有兩點在置信帶外；當(dāng)N取50時，整條覆蓋率連線全部進(jìn)入了置信帶。圖1n的不同取值下的95%的置信區(qū)間覆蓋率圖2n的不同取值下σ的95%的置信區(qū)間覆蓋率圖3n的不同取值下的95%的置信區(qū)間覆蓋率圖4N取5，30，50時，的95%的置信區(qū)間覆蓋率4.3.2方差協(xié)方差考慮σ和，的變化參照。①σ的估計情況。圖5顯示出，N為5時，不論n如何增加，參數(shù)σ置信區(qū)間的覆蓋率都遠(yuǎn)離置信帶；當(dāng)N取30時，大部分的覆蓋率都進(jìn)入了置信帶；當(dāng)N取50時，除一點外，整條覆蓋率連線全部進(jìn)入了置信帶。②的估計情況。圖6表明的估計和固定效應(yīng)參數(shù)與σ有很大差別，當(dāng)N取50時，還有很多覆蓋率在置信帶以外，這表明對于方差協(xié)方差成分的估計N為50是不夠的。結(jié)合圖3，我們可以看到對于方差協(xié)方差成分的估計，N達(dá)到70或者更多才可靠。另外，圖6還顯示了一個有趣的現(xiàn)象，當(dāng)N取30和50時，較小的n（比如n＜20）的覆蓋率穩(wěn)定地進(jìn)入了置信帶，而較大的n（例如：30＜n＜90）的覆蓋率反而在置信帶以外，這表明這時較小的n的參數(shù)估計比較大的n更可靠。5研究結(jié)論從以上的深入分析中，我們可以看到分層線性模型估計的可靠性與樣本中的組數(shù)，組內(nèi)個體樣本數(shù)的取值有密切關(guān)系，總結(jié)如下：(1)樣本中的組數(shù)N。分層線性模型估計的可靠性很大程度上依賴于N的取值，隨著N的增加，固定效應(yīng)，方差協(xié)方差成分的參數(shù)估計的95%的置信區(qū)間的真實覆蓋率能夠進(jìn)入置信帶。當(dāng)N較小時（如N=5），