楊新平，徐登國

（楚雄師范學院 數(shù)學與統(tǒng)計學院，云南 楚雄 675000）

0 引言

線性模型y=Xβ+e的GLS估計，作為一種較為經(jīng)典的估計方法，國內(nèi)外許多學者從各方面做了大量研究。Magnus（1978）[1]在較弱的條件下給出了一階似然估計的條件，并且證明了GLS估計及信息陣的一些性質(zhì)。Magnus等（1980）[2]又借助排除矩陣再次探討了這一問題。對于更一般的混和線性模型，固定效應參數(shù)的似然估計方法與線性模型參數(shù)的似然估計方法類似，只是方差參數(shù)的GLS估計更復雜，需采用Newton-Raphson迭代算法或者Fisher-Scoring迭代算法進行處理。二水平線性模型作為混合線性模型中的一種，在社會科學和行為科學中得到了廣泛的應用，它的參數(shù)估計也受到了普遍的關注。從理論上看，絕大多數(shù)研究都將二水平線性模型納入混和線性模型范疇，應用矩陣直接和的相關知識，將水平2的各組觀測值作為對角塊，給出了精練的估計表達式[3,4]。對于初學者來講，估計的公式較抽象，難理解。如果對水平2采用分組（塊）探討似然估計，雖然得到的公式復雜，但推導過程直觀性強，容易理解。對此問題的研究，Raudenbush和Bryk（2002）[5]進行了簡略陳述。該文針對這一問題，做了細致的推導，清楚地表達了相應結論，導出了一階似然條件和Fisher-Scoring迭代算法的公式和調(diào)整的參數(shù)估計方法，用JSP數(shù)據(jù)作為例子進行了闡述和說明。關于JSP數(shù)據(jù)，Browne（2011）和 Goldstein（2015）[6]曾用方差成分模型進行了研究，并給出IGLS估計。石磊等（2013）[4]也進行了更深入研究，給出相應估計。本文先建立兩水平對中模型，給出模型相應參數(shù)的IGLS估計，再對參數(shù)的估計進行調(diào)整，提高了各個參數(shù)的估計精度。

1 模型的似然估計

對于兩個水平層次數(shù)據(jù)，用l表示水平1的第l個觀測單位，j表示水平2的第j個分組，每組有nj個觀測值，j=1，2，···，J。建立兩水平模型時，先對水平1建模得到模型：

然后將αj進行隨機處理，并表示為水平2相應解釋變量的一般化模型：

其中，Wj是d×p階矩陣，其元素由數(shù)據(jù)中的元素構成，β是p×1向量，uj～N(0，Σ)，并且滿足：Eujrj=0，β和uj，rj獨立，Σ 是d×d階正定矩陣。將式（2）代入式（1）得到一般的模型表達式：

引理1：在模型（3）的假設條件及β，Vj二階可微的條件下，有如下一階ML條件：

其中：γ=(2π)-n/2。

將式（4）代入式（5）得：

對于極大似然值θ，β，必須滿足dΛ=0。但dθ≠0，dβ≠0。因此得到一階似然條件方程：

引理2：在模型（4）的假設及β，Vj二階可微的條件下，參數(shù)β和θ信息矩陣是：

求二階微分：

所以：

式（4）代入式（6）得：

所以信息矩陣為：

定理：模型（4）的假設及β，Vj二階可微的條件下，參數(shù)的IGLS估計可通過下列三個迭代方程得到：

其中：Dv(Σ(l))=vec(Σ(l))，給定初始值，根據(jù)v(Σ0)意義，用v(Σ0)的元素構成 Σ0，始值 Σ0滿足維列向量，階方陣，它們分別由引理1及引理2計算得到，即：

證明：對于方程（7），由引理1及引理2的證明容易得到向量及矩陣A((v(Σ))′，(σ2)|Y，U，Z)，分別計算它們在第l步的值，再用計算出的第l步的值分別替換Newton-Raphson迭代公式中右邊相應各項就得到方程（7）。對于方程（8），由于使用公式容易得到：寫成迭代表達式即證。對于方 程（8），由 引 理 1，有得正規(guī)方程滿足滿秩的條件，則寫成迭代表達式即證。

在迭代方程中，由方程（7）通過l步計算得到的Σ的更新估計并不能保證其正定性，因此，對于組間觀測數(shù)量差異較大的數(shù)據(jù)，建立二水平線性模型后，并用迭代廣義最小二乘估計方法得到的水平2協(xié)方差估計不一定正定。如果是負定矩陣，在進行bootstrap抽樣之前，需進行估計值的調(diào)整，使得協(xié)方差陣估計正定。

2 估計調(diào)整

3 實例

JSP數(shù)據(jù)來源于倫敦市48所小學728名學生，是數(shù)據(jù)集的一個子樣，在收集數(shù)據(jù)時，同一個學生考察在兩個時段的成績，變量名分別是：math.yr.3（記為yij），math.yr.1（記為x1ij），Gender.boy.1（記為x2ij），Social.class.manual.1（記為x3ij），School.ID。建立對中兩水平線性模型：

表1 JSP數(shù)據(jù)的二水平對中模型（10）的參數(shù)估計及調(diào)整后的參數(shù)估計匯總表

圖1 核密度估計，均值2.471(0.328),眾數(shù)3.332

模型（10）參數(shù)在Bootstrap鏈中的收斂性可通過由它的均值鏈的穩(wěn)定性來判斷，這里只給出-2loglikelihood、固定效應參數(shù)β1和各個隨機效應參數(shù)的均值鏈曲線（見圖2至圖7），其他均值鏈曲線略。

圖2 -2loglikelihood=4293.090的均值曲線

圖3 β1=0.590(0.026)的均值曲線

圖4 =2.471(3.884)的均值曲線

圖5 =-0.103(0.103)的均值曲線

圖6 =0.038(0.026)的均值曲線

圖7 =17.794(4.296)的均值曲線

從各自的均值鏈曲線可知，在5個重復集的情況，迭代次數(shù)達到80次以上，固定效應參數(shù)β1、-2loglikelihood值、隨機效應參數(shù)的穩(wěn)定性較好，而的穩(wěn)定性稍差的核密度估計接近正態(tài)分布密度（見圖8），用正態(tài)分布密度作近似分布得到的置信度為：

圖8 核密度估計，均值17.794(4.296),眾數(shù)17.454

95%的置信區(qū)間為（9.32884，26.16916），其他參數(shù)的漸近區(qū)間估計同樣計算。如果要提高穩(wěn)定性，減少模擬的誤差，可以通過增加重復集的數(shù)目和每個重復集的迭代次數(shù)實現(xiàn)，但會增加較大的運算量，程序運行需要大量時間。

由表1中的計算結果可知：經(jīng)過調(diào)整計算后的固定效應參數(shù)只是在社會背景上有變化，其他值變化不大，特別是固定效應的截距項均為31.4，說明若某一名學生在8歲時的數(shù)學績?nèi)绻苓_到平均成績25.97，在社會背景不變的條件下，那么在11歲時其期望的數(shù)學成績?yōu)?1.4分。-2loglikelihood由4235.920增加到4293.09值保持穩(wěn)定的值有一定的改變，固定效應參數(shù)的標準誤差減小，估計精度有所提高。

4 結論

（1）本文前半部分采用分組（塊）來探討似然估計，經(jīng)過詳細的推導，公式雖然相對復雜，但直觀性較強，對于初學者來說，較容易理解。兩水平線性模型在正態(tài)性假定及滿足引理1中三個條件的情況下，極大化對數(shù)似然函數(shù)首先得到固定效應的參數(shù)估計，再通過定理中的迭代方程（7）至方程（9）進行迭代計算，直至各個參數(shù)的序列收斂，就得到參數(shù)的IGLS估計。已經(jīng)證明：在正態(tài)性假定下，IGLS估計等價于極大似然估計ML[9]。IGLS估計是進一步實現(xiàn)Bootstap抽樣的基礎，通過自助抽樣，得到各個待估參數(shù)的自助樣本，求出各個參數(shù)的核密度估計，可完成進一步統(tǒng)計推斷。

（2）實例中用到的JSP數(shù)據(jù)，共有728個樣本數(shù)據(jù)，分成48個組，第46所學校有44個觀測值，分組觀測值個數(shù)最多，第42所學校有4個觀測值，分組觀測值個數(shù)最少，數(shù)據(jù)不平衡，用經(jīng)典的IGLS估計模型（10）的參數(shù)，得到的水平2的協(xié)方差估計負定，結論適用性差，通過用Bootstrap方法進行重抽樣求IGLS估計，估計值得到了改進，表現(xiàn)在獲得了正定的水平2協(xié)方差估計，固定效應的參數(shù)估計的標準差比IGLS的標準差小，同時還可以利用核密度估計進行統(tǒng)計推斷。因此，在經(jīng)典的統(tǒng)計方法的基礎上，采用自助抽樣或者貝葉斯方法改進參數(shù)估計精度的方法是值得推薦的。