基于混合地理加權(quán)Fay-Herriot模型的小域估計

李騰,魏傳華,于力超

( 中央民族大學(xué)理學(xué)院統(tǒng)計學(xué)系,北京100081)

1.引言

近二十年來,小域估計(Small area estimation)作為抽樣調(diào)查領(lǐng)域的一種方法,得到越來越多的重視.理論上,小域估計模型的設(shè)定以及有關(guān)的估計和檢驗等統(tǒng)計推斷問題得到了深入的研究.應(yīng)用上,小域估計方法已經(jīng)在經(jīng)濟、社會、環(huán)境、衛(wèi)生和流行病等多個領(lǐng)域的實際問題中得到應(yīng)用.有關(guān)小域估計的詳細(xì)介紹可參考著作[14].

小域估計方法主要分為基于設(shè)計的方法和基于模型兩類方法,其中基于模型的方法能夠充分利用到輔助信息而得到重視.一般來說,基于模型的小域估計方法所使用的模型一般分為區(qū)域?qū)哟文Ｐ秃蛦卧獙哟文Ｐ蛢深?其中區(qū)域?qū)哟文Ｐ褪褂幂^多的是Fay和Herriot[7]提出的能將直接估計和未知參數(shù)聯(lián)系起來的一類隨機效應(yīng)模型,文獻(xiàn)稱之為Fay-Herriot模型,關(guān)于該模型的理論研究可參考文[1,9,16,21].

目前關(guān)于Fay-Herriot模型的研究和使用上大都假定了各個區(qū)域之間是不相關(guān)的,同時使用的是線性混合效應(yīng)模型.然而，實際問題中因變量與自變量之間的關(guān)系不一定是線性形式,錯誤的模型設(shè)定會導(dǎo)致推斷偏離實際情況,為了解決這一問題,多種非參數(shù)與半?yún)?shù)混合效應(yīng)模型被用到小域估計中,具體內(nèi)容可參考文[10,17].另外在很多類實際問題中,抽樣區(qū)域之間并不是獨立的,而是具有空間效應(yīng),從而多類空間混合效應(yīng)模型被提出用以小域估計,有關(guān)文獻(xiàn)可參考文[3,11,15].

空間效應(yīng)一般分為空間自相關(guān)性和空間非平穩(wěn)性兩類,以往的研究大都是基于空間自相關(guān)性來建立空間混合效應(yīng)模型.Chandra等[4]基于地理加權(quán)回歸方法提出了一類空間變系數(shù)混合效應(yīng)模型.地理加權(quán)回歸模型是由Brunsdon等[2]提出來的一類空間變系數(shù)模型,模型系數(shù)作為地理位置比如經(jīng)度緯度的光滑函數(shù),可以直接刻畫因變量和自變量之間的關(guān)系隨著地理位置的變化而變化,關(guān)于地理加權(quán)回歸的詳細(xì)介紹可參考著作[8].實際數(shù)據(jù)分析中，更為普遍的情形是一部分系數(shù)是地理位置的光滑函數(shù),另外一部分是常數(shù),這種半?yún)?shù)空間變系數(shù)模型被稱為混合地理加權(quán)回歸模型,關(guān)于這類模型的研究可參考文[13,18―20].為了更好的刻畫抽樣區(qū)域之間的空間非平穩(wěn)性,本文將在傳統(tǒng)Fay-Herriot模型的基礎(chǔ)上提出一類半?yún)?shù)空間變系數(shù)混合效應(yīng)模型,并研究相應(yīng)的小域估計量及其性質(zhì).

本文剩余部分做如下安排:第二節(jié)主要介紹半?yún)?shù)空間變系數(shù)混合效應(yīng)Fay-Herriot模型的估計;第三節(jié)將進(jìn)行小域估計量的均方誤差估計;第四節(jié)對所提方法進(jìn)行數(shù)值模擬;總結(jié)和展望將在最后一節(jié)給出.

2.模型的估計

假定目標(biāo)總體U由m個小區(qū)域構(gòu)成,且各小區(qū)域間都不重疊,每個區(qū)域可看做一個單元.若將第i個區(qū)域(i=1,2···m)的目標(biāo)參數(shù)記為Yi,xi和zi為輔助變量.其中,xi的系數(shù)β是固定的,zi的系數(shù)α(ui,vi) 隨地理位置(ui,vi) 變化而變化(u和v分別表示研究區(qū)域i的經(jīng)度和緯度),bi是隨機效應(yīng),yi為Yi的直接估計量,ei為抽樣誤差.根據(jù)上面的符號,區(qū)域?qū)哟蔚陌雲(yún)?shù)空間變系數(shù)混合效應(yīng)模型可記為如下形式

其中β為p×1維向量,α(ui,vi) 是q×1 維向量.常值系數(shù)β和系數(shù)函數(shù)α(ui,vi)對應(yīng)于固定效應(yīng).bi為隨機效應(yīng),且它的均值為0,方差為σ2b,隨機誤差項ei ～N(0,Ψi),其中Ψi已知.同時假設(shè)不同區(qū)域間的隨機誤差項是相互獨立的,Ψ= diag1≤i≤m(Ψi).為了將模型表達(dá)得更加清晰,將模型(2.1)記為如下矩陣形式

其中

y的協(xié)方差矩陣記為COV(b+e)=Σ,且有Σ=σ2bIm+Ψ.

現(xiàn)運用文[6,18]中的Profile最小二乘方法估對模型(2.1) 進(jìn)行估計.模型中常值系數(shù)假設(shè)已知,那么(2.1)可記為

進(jìn)一步模型可表示為以下的矩陣形式

Chandra等[4]所研究的空間變系數(shù)混合效應(yīng)模型便是模型(2.3)這種形式,同時空間變系數(shù)混合效應(yīng)模型中,Chandra等用局部常數(shù)方法來估計系數(shù)函數(shù).本部分將采用局部線性方法對模型(2.3)中系數(shù)函數(shù)進(jìn)行估計.設(shè)研究區(qū)域中某一點為(u0,v0),將(u0,v0) 到第i區(qū)域的距離設(shè)為d0i.利用Brunsdon等[2]在(u0,v0) 點構(gòu)造一組權(quán),使wi(u0,v0) =k(d0i/h),其中k(·) 為核函數(shù),h為光滑參數(shù)(或窗寬),我們將在下文介紹.

根據(jù)Taylor展開公式,在(u0,v0) 鄰域內(nèi)有

其中α(u)(u0,v0) 和α(v)(u0,v0) 分別表示α(ui,vi) 關(guān)于經(jīng)度u和緯度v的偏導(dǎo)數(shù)在(u0,v0)處的值.于是α(u,v) 在(u0,v0) 處的估計可通過使得

達(dá)到最小得到.由加權(quán)最小二乘法可得

其中Γ(u0,v0) = (Z,U(u0,v0)Z,V(u0,v0)Z),Z= (zT1,zT2,··· ,zTm)T,U(u0,v0) = Diag(u1?u0,u2?u0,··· ,um?u0),V(u0,v0)=Diag(v1?v0,v2?v0,··· ,vm?v0),W(u0,v0)=Diag(w1(u0,v0),w2(u0,v0),··· ,wm(u0,v0)).如果取(u0,v0) 為(ui,vi) ,則由(2.5)式得α(u,v) 在各樣本點的估計為

從而M的估計可定義為

其中

將式(2.4)中M用其估計值代替,整理可得如下線性混合效應(yīng)模型

整理得

從而定義β的廣義profile最小二乘估計為

由混合線性模型的理論可知b的經(jīng)驗最佳線性無偏估計量(EBLUP)為

那么第i個區(qū)域的隨機效應(yīng)的EBLUP表示為

其中l(wèi)Ti為1×n維向量(0,0,··· ,0,1,0,··· ,0) ,在第i個位置取值為1,其他位置取值為0.

在對(u0,v0)處的模型系數(shù)函數(shù)值α(u0,v0)進(jìn)行估計時候,每個樣本點(ui,vi) 都對應(yīng)一個權(quán)重wi(u0,v0).一般來說,距離(u0,v0) 近的觀測值對(u0,v0) 的回歸函數(shù)(或參數(shù))估計影響程度大,距離遠(yuǎn)的觀測值對其影響程度小.因此,設(shè)d0i為(u0,v0)到(ui,vi) 的距離,我們應(yīng)將較大的權(quán)重賦予距離近的點,較小的權(quán)重賦予距離遠(yuǎn)的點,實際應(yīng)用中可以采用Gauss核函數(shù)作為權(quán)函數(shù)

其中h稱為光滑參數(shù)或窗寬,反映了擬合曲線的光滑性.h的大小通過“去一觀測主題”交叉證實法來確定.

現(xiàn)在可得到小域估計中第i個未抽區(qū)域的目標(biāo)參數(shù)估計量,即地理加權(quán)經(jīng)驗最佳線性無偏估計量.記為.該區(qū)域目標(biāo)參數(shù)可根據(jù)模型參數(shù)估計量和本身的輔助變量進(jìn)行估計,的具體表達(dá)式如下:

由于參數(shù)β,α(ui,vi)和bi的估計需要σ2b和Ψ作為輔助信息,但協(xié)方差參數(shù)θ=σ2b是未知的,我們可以類似于文[4],利用極大似然估計(ML)方法或約束的極大似然估計(REML)方法來對θ=σ2b進(jìn)行估計,REML方法是ML方法的修正,它考慮了ML 方法估計固定效應(yīng)造成自由度損失的問題.值得注意的是這兩種方法其實都是profile似然估計方法,也就是將模型中的非參數(shù)部分消去,從而只剩下參數(shù)分量.由于估計方法不同,兩種方法得到的估計值也是不同的.

ML算法

1) 計算各區(qū)域與樣本點間的距離d;

2) 計算各區(qū)域?qū)?yīng)的權(quán)重矩陣W(ui,vi);

3) 通過廣義交叉驗證確定光滑參數(shù)h;

4) 計算對數(shù)似然函數(shù)

將上式中的M用(2.7)中?M代替,可得如下的對數(shù)似然函數(shù)

5) 對l(β,α,θ)關(guān)于θ求一階偏導(dǎo),得s(β,α,θ)

其中Σ(θ)=?Σ/?θ=Im,Σ(θ)=?Σ?1/?θ=?Σ?1Σ?1;

6) 對-l(β,α,θ)關(guān)于θ求二階偏導(dǎo),得I(θ)

7) 假定θ初始值為0,通過得分算法對下式進(jìn)行迭代,這里的a指第a次迭代過程.

8) 當(dāng)?shù)諗繒r，我們得到θ的估計值

REML算法

1) 計算各區(qū)域與樣本點間的距離d;

2) 計算各點對應(yīng)的權(quán)重矩陣W(ui,vi);

3) 通過廣義交叉驗證確定光滑參數(shù)h;

4) 計算極大似然函數(shù)

5) 對lR(β,α,θ)關(guān)于θ求一階偏導(dǎo),得sR(β,α,θ):

6) 并對?lR(β,α,θ)關(guān)于θ求二階偏導(dǎo),得IR(θ),

7) 假定θ初始值為0.通過得分算法對下式進(jìn)行迭代,這里的a指第a次迭代過程,

8) 當(dāng)?shù)諗繒r,我們得到θ的約束極大似然估計值

3.小域估計量的均方誤差估計

我們采用文[12]對混合效應(yīng)模型均方誤差的求解方法對小域估計量的均方誤差進(jìn)行求解.得到如下結(jié)論:

定理3.1小域估計量是抽樣集y的線性組合.兩者有如下關(guān)系:A、B和D分別為

證

定理3.2小域估計量的均方誤差其中

g3i(θ) =的漸進(jìn)協(xié)方差,其中A、B含義同(3.1)式.

證和式(2.13)類似,當(dāng)協(xié)方差變量θ=σb2已知時,我們可得到該估計量的地理加權(quán)最佳線性無偏估計量記t(θ) .

假設(shè)β真值已知,通過β真值估計的α(ui,vi) 和bi來表示該估計量設(shè)為t?(θ).

假設(shè)β，α(ui,vi)和bi真值均已知,那么該統(tǒng)計量為

根據(jù)小域估計量的MSE定義式得到

式中第一項為已知β時小域估計量的方差,記為g1i(θ);第二項為回歸系數(shù)時引起的方差部分,記為g2i(θ).

式(3.5)中

其中

則

式(3.5)中

因此有

因為θ=σ2b未知,我們需要用估計值代替真值進(jìn)行均方誤差求解.所以

由于式(3.8)很難直接處理,因此借助文[5]中Taylor展開的方法得到近似估計,記d(θ) =得

首先采用REML算法對小域估計量的均方誤差進(jìn)行估計.由Rao(1975)我們知

結(jié)合式(3.11)有

且

4.數(shù)值模擬

本部分利用數(shù)值模擬來考察所提模型的有效性.為了進(jìn)行比較,數(shù)值模擬實驗假定在協(xié)方差參數(shù)已知和協(xié)方差參數(shù)需要估計的情況下,分別對模型參數(shù)進(jìn)行估計.

模擬中假定模型形式為

其中i=1,...,m2,研究區(qū)域是邊長m ?1個單位長度的正方形,n=m2個樣本觀測值正好落在分布均勻的m×m格子點上,則每個觀測點的經(jīng)緯度坐標(biāo)(ui,vi) 為

其中mod(i ?1,m)為i ?1與m相除后的余數(shù),的整數(shù)值部分.

假定固定系數(shù)β= 3,固定系數(shù)所對應(yīng)的變量xi滿足xi ～N(1,1),變系數(shù)所對應(yīng)的變量zi滿足zi ～U(0,1).其中zi的系數(shù)與小域的經(jīng)緯度(ui,vi) 是相關(guān)的,滿足α(ui,vi)=ui+vi.此外我們假定隨機效應(yīng)bi服從正態(tài)分布N(0,0.72),隨機擾動項ei服從正態(tài)分布N(0,0.82),值得注意的是,隨機擾動項的條件我們是已知的.模擬中采用Gauss 核函數(shù)作為權(quán)函數(shù),窗寬h通過交叉驗證法來確定.

設(shè)定模擬重復(fù)次數(shù)為T=500,樣本容量m2分別為36,64和100.利用前面介紹的估計方法,基于以上數(shù)據(jù)分三種情況(協(xié)方差參數(shù)已知,ML方法估計協(xié)方差,REML 方法估計協(xié)方差)給出參數(shù)的估計.同時我們以符號“MEAN”“SD”分別表示固定系數(shù)進(jìn)行500 次模擬實驗的估計值均值、標(biāo)準(zhǔn)差.同時記yit和分別為樣本中區(qū)域i第t次模擬的真實值和估計值.選取指標(biāo)(平均絕對相對誤差)和(平均相對根均方誤差)來評估估計量的優(yōu)劣.指標(biāo)定義如下:

其中,

與此同時,我們給出三種情形下估計值在次模擬下的均方誤差均值.模擬結(jié)果見表1至表3:

表1 參數(shù)β的估計結(jié)果

表2 小域估計量的指標(biāo)對比

表3 小域估計量均方誤差均值的比較

我們將協(xié)方差已知、協(xié)方差未知(ML算法估計)、協(xié)方差未知(REML算法估計)情況下β的估計值記為?β.

表1為參數(shù)β在協(xié)方差參數(shù)σ2b已知、用ML方法估計協(xié)方差參數(shù)、用REML 方法三種情況下的估計結(jié)果.由此表可看出: 1)整體上,500次模擬中參數(shù)β的估計值趨近于真值3,且在真值3附近上下波動.估計值的標(biāo)準(zhǔn)差波動幅度較小,且均在區(qū)間[0,0.30]內(nèi)浮動;2)當(dāng)協(xié)方差參數(shù)σ2b已知時,估計值?β0隨著樣本量的增大,越來越接近真值3,標(biāo)準(zhǔn)差也隨著樣本量的增大逐漸減小.當(dāng)協(xié)方差參數(shù)σ2b未知時,通過采用ML 方法或者REML 方法估計協(xié)方差參數(shù),我們同樣發(fā)現(xiàn)隨著樣本量的增加,在相同方法中所得到的參數(shù)β估計值逐漸趨近于真值3,同時估計值的標(biāo)準(zhǔn)差逐漸減小.由此說明,樣本量是影響參數(shù)估計的重要因素,樣本量越大,參數(shù)的估計越接近真實值,同時也越穩(wěn)定;3)當(dāng)樣本量相同時,協(xié)方差參數(shù)σ2b已知時得到的參數(shù)估計值最接近真值3,ML 方法估計協(xié)方差參數(shù)得到的參數(shù)估計值與REML方法估計協(xié)方差參數(shù)得到的參數(shù)估計值相比,更接近真值3.在標(biāo)準(zhǔn)差方面,REML 方法估計協(xié)方差參數(shù)所得到的標(biāo)準(zhǔn)差偏小一些,ML所對應(yīng)的標(biāo)準(zhǔn)差偏大一些.這說明REML方法得到的估計值更加穩(wěn)定.

表2為在協(xié)方差參數(shù)已知和協(xié)方差參數(shù)需要估計的三種情況下,所得小域估計量的平均絕對相對誤差和平均相對根均方誤差指標(biāo)比較.兩個指標(biāo)反映了小域估計值的擬合程度.樣本量相同時,協(xié)方差參數(shù)已知的ARE小于協(xié)方差參數(shù)未知時的ARE.同樣協(xié)方差參數(shù)已知得到的小于協(xié)方差參數(shù)未知時的.相同條件下,ML與REML所得到的指標(biāo)結(jié)果相差不大,REML方法得到的指標(biāo)相對更偏小一些.此外,隨著樣本量的增加,逐漸遞減,同時兩種估計方法所得到的指標(biāo)值越來越接近,說明兩種估計方法對協(xié)方差參數(shù)σ2b估計效果相似.

表3是協(xié)方差參數(shù)已知、協(xié)方差參數(shù)需估計這三種情況下小域估計量n個區(qū)域T次模擬中均方誤差均值的對比(確切的說,當(dāng)協(xié)方差參數(shù)σ2b已知時,得到的是小域估計量均方誤差.當(dāng)協(xié)方差參數(shù)未知需用ML算法和REML算法估計時,得到的是小域計量均方誤差的估計值.) 1)當(dāng)協(xié)方差參數(shù)σ2b已知時,均方誤差MSE(y)分為g1(y)和g2(y),g1(y)的取值反映了隨機效應(yīng)的估計影響,β的估計與g2(y)相關(guān).當(dāng)協(xié)方差參數(shù)σ2b需要估計時,g3(y) 反映的是協(xié)方差參數(shù)估計對總體均方誤差估計的影響.2)g1(y)取值最大,對MSE(y)影響最大,說明隨機效應(yīng)估計對均方誤差影響最大,這是由于隨機效應(yīng)設(shè)定時的標(biāo)準(zhǔn)差較大,存在一定的波動性,從而產(chǎn)生了較大的g1(y).g2(y)取值較小,說明β的估計對均方誤差的影響較小,即β的估計是無偏的.g3(y) 有一定的取值,說明ML算法和REML 算法對協(xié)方差參數(shù)進(jìn)行估計時有一定的偏差.3)協(xié)方差參數(shù)σ2b已知時,g1(y) 對均方誤差的影響占比較大,協(xié)方差參數(shù)σ2b未知時,影響均方誤差的主要是g1(y)和g3(y) 兩個部分.4)協(xié)方差參數(shù)σ2b未知時,ML 算法的均方誤差估計值與REML 的均方誤差估計值差別不大,說明兩種算法的估計沒有明顯差距.5)三種情形下,均方誤差的估計值隨著樣本量的增大逐漸減小.同時三部分的取值也逐漸減少.

5.總結(jié)

為了更好的刻畫區(qū)域之間的空間效應(yīng),本文提出了一類基于混合地理加權(quán)的Fay-Herriot模型用以小域估計.論文給出了目標(biāo)參數(shù)的估計,并研究了其均方誤差的估計問題.論文結(jié)果推廣了小域估計的現(xiàn)有模型.