兩水平線性模型的參數估計及調整

楊新平，徐登國

（楚雄師范學院 數學與統計學院，云南 楚雄 675000）

0 引言

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

1 模型的似然估計

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

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

其中，Wj是d×p階矩陣，其元素由數據中的元素構成，β是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二階可微的條件下，參數β和θ信息矩陣是：

求二階微分：

所以：

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

所以信息矩陣為：

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

其中：Dv(Σ(l))=vec(Σ(l))，給定初始值，根據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，有得正規方程滿足滿秩的條件，則寫成迭代表達式即證。

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

2 估計調整

3 實例

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

表1 JSP數據的二水平對中模型（10）的參數估計及調整后的參數估計匯總表

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

模型（10）參數在Bootstrap鏈中的收斂性可通過由它的均值鏈的穩定性來判斷，這里只給出-2loglikelihood、固定效應參數β1和各個隨機效應參數的均值鏈曲線（見圖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個重復集的情況，迭代次數達到80次以上，固定效應參數β1、-2loglikelihood值、隨機效應參數的穩定性較好，而的穩定性稍差的核密度估計接近正態分布密度（見圖8），用正態分布密度作近似分布得到的置信度為：

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

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

由表1中的計算結果可知：經過調整計算后的固定效應參數只是在社會背景上有變化，其他值變化不大，特別是固定效應的截距項均為31.4，說明若某一名學生在8歲時的數學績如果能達到平均成績25.97，在社會背景不變的條件下，那么在11歲時其期望的數學成績為31.4分。-2loglikelihood由4235.920增加到4293.09值保持穩定的值有一定的改變，固定效應參數的標準誤差減小，估計精度有所提高。

4 結論

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

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