基于混合地理加權Fay-Herriot模型的小域估計

李騰,魏傳華,于力超

( 中央民族大學理學院統計學系,北京100081)

1.引言

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

小域估計方法主要分為基于設計的方法和基于模型兩類方法,其中基于模型的方法能夠充分利用到輔助信息而得到重視.一般來說,基于模型的小域估計方法所使用的模型一般分為區域層次模型和單元層次模型兩類,其中區域層次模型使用較多的是Fay和Herriot[7]提出的能將直接估計和未知參數聯系起來的一類隨機效應模型,文獻稱之為Fay-Herriot模型,關于該模型的理論研究可參考文[1,9,16,21].

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

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

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

2.模型的估計

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

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

其中

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

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

進一步模型可表示為以下的矩陣形式

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

根據Taylor展開公式,在(u0,v0) 鄰域內有

其中α(u)(u0,v0) 和α(v)(u0,v0) 分別表示α(ui,vi) 關于經度u和緯度v的偏導數在(u0,v0)處的值.于是α(u,v) 在(u0,v0) 處的估計可通過使得

達到最小得到.由加權最小二乘法可得

其中Γ(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用其估計值代替,整理可得如下線性混合效應模型

整理得

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

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

那么第i個區域的隨機效應的EBLUP表示為

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

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

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

現在可得到小域估計中第i個未抽區域的目標參數估計量,即地理加權經驗最佳線性無偏估計量.記為.該區域目標參數可根據模型參數估計量和本身的輔助變量進行估計,的具體表達式如下:

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

ML算法

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

2) 計算各區域對應的權重矩陣W(ui,vi);

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

4) 計算對數似然函數

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

5) 對l(β,α,θ)關于θ求一階偏導,得s(β,α,θ)

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

6) 對-l(β,α,θ)關于θ求二階偏導,得I(θ)

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

8) 當迭代收斂時，我們得到θ的估計值

REML算法

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

2) 計算各點對應的權重矩陣W(ui,vi);

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

4) 計算極大似然函數

5) 對lR(β,α,θ)關于θ求一階偏導,得sR(β,α,θ):

6) 并對?lR(β,α,θ)關于θ求二階偏導,得IR(θ),

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

8) 當迭代收斂時,我們得到θ的約束極大似然估計值

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

我們采用文[12]對混合效應模型均方誤差的求解方法對小域估計量的均方誤差進行求解.得到如下結論:

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

證

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

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

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

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

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

根據小域估計量的MSE定義式得到

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

式(3.5)中

其中

則

式(3.5)中

因此有

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

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

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

結合式(3.11)有

且

4.數值模擬

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

模擬中假定模型形式為

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

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

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

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

其中,

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

表1 參數β的估計結果

表2 小域估計量的指標對比

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

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

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

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

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

5.總結

為了更好的刻畫區域之間的空間效應,本文提出了一類基于混合地理加權的Fay-Herriot模型用以小域估計.論文給出了目標參數的估計,并研究了其均方誤差的估計問題.論文結果推廣了小域估計的現有模型.