基于傾向得分多層模型的非概率抽樣統計推斷

劉 展

（湖北大學a.數學與統計學學院；b.應用數學湖北省重點實驗室，武漢430062）

0 引言

網絡調查由于其具有比電話調查更短的數據收集周期和更低的費用而急劇增長，目前已經在市場研究的調查數據收集中占主導地位。網絡調查可視為一種單一的模式，但其抽取樣本的方法有多種，包括基于概率抽樣的網絡調查與基于非概率抽樣的網絡調查，其中最值得注意的是基于非概率抽樣的候選者數據庫網絡調查。這里網絡調查的候選者數據庫，簡稱網絡候選者數據庫，指自愿完成網絡調查的上網人群，如果在后續的網絡調查中選擇他們作為調查對象，他們將配合完成調查[1]。候選者數據庫網絡調查是從網絡候選者數據庫中抽取樣本進行網絡調查，其獲得的樣本為網絡候選者數據庫的調查樣本?？梢钥吹?，候選者數據庫網絡調查本質上屬于非概率抽樣調查，其樣本屬于非概率樣本，入樣概率未知，采用傳統的抽樣推斷理論進行推斷存在一定的困難。因此，非概率樣本的統計推斷問題，特別是候選者數據庫網絡調查的推斷問題，就成為網絡調查發展中迫切需要解決的問題。

傾向得分調整方法是一種被廣泛應用于非概率抽樣推斷的方法。Duncan和Stasny（2001）[2]利用傾向得分方法降低了電話調查的覆蓋偏差。Schonlau等（2004）[3]與Lee（2006）[4]指出傾向得分調整方法可能有助于減少抽樣框覆蓋不全、無回答、非概率抽樣等產生的偏差。Brick（2013）[5]提出使用傾向得分加權來解決單元無回答的缺失問題。Schonlau 等（2009）[6]將性別、年齡、教育、收入等作為X，是否上網作為Y建立模型，構建非概率樣本的傾向得分權數，最終得到的加權估計接近于總體參數真值。Valliant和Dever（2011）[7]結合網絡候選者數據庫和參考樣本，建立logistic回歸模型來估計傾向得分，并基于傾向得分對非概率樣本構造了不同的權數來估計總體。然而，以上建立傾向得分模型估計傾向得分來推斷總體均是在單一層次數據背景下進行的，即假設網絡候選者數據庫的數據結構是單一層次的，然而在一些情況下網絡候選者數據庫的數據結構可能是多層次或嵌套結構。多層次或嵌套結構就是在上一級分層的基礎上進行再分層，形成一層套一層的嵌套結構或多層次結構，單一層次結構就是一個層次，不再逐層細分。例如，學生往往在一個班級內可以被分成不同的類，反過來這些班級又在學校內可以被分成不同的群。如果忽略數據的嵌套結構，仍然采用單一層次結構下的方法進行分析，可能會導致誤導性或不準確的結果。因此，當網絡候選者數據庫數據具有嵌套結構時，需要考慮針對嵌套結構數據的傾向得分估計方法。Li等（2013）[8]在觀察性研究中針對多層（嵌套）結構數據，探究了建立多層模型估計傾向得分來估計平均處理效應，但是并未探討總體推斷的問題。關于多層模型在醫療護理、衛生政策、教育等諸多領域中的應用研究文獻越來越多，同時傾向得分分析也日益流行，然而對多層次數據結構的傾向得分分析還沒有得到深入的研究。基于此，本文針對具有嵌套結構數據的網絡候選者數據庫和參考樣本，探索構建傾向得分多層模型來推斷總體的方法。為了便于討論，本文主要考慮兩個層次的嵌套結構數據。

1 多層模型

多層模型（Multilevel Modeling，簡記為MLM）是包含固定效應和（或）隨機效應的回歸模型，是一種用于評價嵌套結構數據的統計分析方法[9]。多層模型是通過解釋單元間的依賴關系、適當調整標準誤差、在所有層次上分解方差，來改進嵌套數據中個體效應的估計。此外，多層模型允許跨層次的相互作用，這種相互作用解釋了在一個層次測量的變量是如何影響在另一層次上所發生的關聯的[9]。換句話說，對嵌套數據建立多層模型允許研究者去調查結果的變動有多少與群內和群間個體之間的差異有關[10]。下面討論具體的多層模型方程。

對第j個群中第i個人觀測的連續因變量Y的多層模型如下：

式（1）中β0、β1是固定效應，u0j、u1j為隨機效應，eij為隨機誤差項。其中固定效應是未知的常數（參數），定義了預測變量X1ij與因變量Yij之間的關系；隨機效應是隨機變量，允許回歸模型中的系數根據研究的群或個體（對象）隨機變動，因此多層模型也稱為隨機系數模型。既然隨機效應是隨機變量，相應的就需要對隨機效應定義概率分布，這就意味著多層模型比一般的回歸模型有更多的參數需要估計。隨機效應與隨機誤差項通常的假定分布為：

其中D為隨機效應的方差協方差矩陣，為u0j的方差，為u1j的方差，σ01為u0j與u1j的協方差，σ2為eij的方差。

事實上，多層模型（1）也可以寫成另外一種形式，如式（3）所示，式（3）是更為常見的一種形式，很多軟件比如HLM、MLwiN等使用的都是這種形式。

其中β0j、β1j為隨機系數，并且β0j為第一層次方程的截距，β1j為第一層次方程的斜率，將第一層次與第二層次的方程結合即可得到式（1），第二層次中關于β0j的方程稱為截距方程，關于β1j的方程稱為斜率方程。式（3）清晰地定義了多層模型中第二層次所測量的第一層次預測變量，可以將第二層次的每個隨機系數方程視為一個僅有截距項的回歸模型，并且可以通過在這個模型中增加第二層次預測變量來解釋隨機效應的方差，比如在截距方程中增加性別Gender作為第二層次預測變量，則β0j=β0+β2Genderj+u0j，β2為未知常數（固定效應）。

式（3）只是一個簡單的二層次模型，在具體建模時可根據具體的問題在第一層次方程中添加預測變量如X2ij，…，Xpij等，在第二層次方程中也可以添加一些預測變量。總之，在多層模型中預測變量可以被添加到第一層次和第二層次的模型方程中，這些變量的截距或斜率可以跨層次的變動，即截距或斜率可以為隨機變量，是否限制或允許第一層次與第二層次預測變量跨層次變動是多層模型特殊的一個重要方面[9,11,12]。如果為三層次模型，則第二層次方程中的一些系數也是隨機的，可以同樣寫出隨機系數的第三層次方程，依次可類推到多個層次的模型。

對于多層模型的估計，可以從矩陣形式來考慮。在隨機效應的條件下，第 j個對象或群的觀測向量分布為 Yj|uj～N(Xjβ+Zjuj，Rj)，其中Xj為預測變量設計矩陣，β為固定效應向量，Zj為更小一些的預測變量設計矩陣，uj為隨機效應向量，Rj為隨機誤差的方差協方差矩陣。所有隨機效應平均下Y的邊際分布為Vj)，其中Xjβ為第j個對象或群的觀測Y的期望向量，D為隨機效應的方差協方差矩陣對象或群的觀測Y的邊際方差協方差矩陣。在此基礎上，采用某種形式的極大似然估計方法來估計模型參數。假設對象或群是相互獨立的，首先寫出似然函數形式為的邊際概率密度函數，β、θ為待估參數，對有些模型來說似然函數可能需要近似。然后，使用迭代數學方法（比如Newton-Raphson）找到使得（近似的）似然函數值達到最大的模型參數估計。

2 傾向得分多層模型及總體估計

假設網絡候選者數據庫的協變量數據具有嵌套結構，比如協變量包括了班級、學校等變量，此時網絡候選者數據庫單元為不同學校不同班級里的學生。從該網絡候選者數據庫中隨機抽取一個樣本進行調查，得到一個網絡候選者數據庫的調查樣本；同時假設獲得了另一個具有同樣嵌套結構數據的參考樣本，比如采用多階段隨機抽樣得到了另一個概率樣本。設調查樣本中第j個群第i個單元的入樣概率 πij=π(Xij)，Xij=(Xij1，Xij2，…，Xijp)′是協向量，p是協變量數目，則有：

其中，V是網絡候選者數據庫；SV是從網絡候選者數據庫中隨機抽取的調查樣本；P(i j ∈ SV|Xij，V)是網絡候選者數據庫中第j個群第i個單元進入網絡調查樣本的概率，常常已知；P(ij∈V|Xij)是第j個群第i個單元自愿加入網絡候選者數據庫的概率，往往未知需要估計。設示性變量D表示單元是否加入網絡候選者數據庫，D=1記為單元加入網絡候選者數據庫，D=0記為單元未加入網絡候選者數據庫，那么第 j個群第i個單元加入網絡候選者數據庫的概率就是p(Xij)=P(Dij=1|Xij)=P(ij∈V|Xij)。

根據傾向得分的定義：給定輔助變量X的條件下單元接受處理（參與或回答）的條件概率，可知p(Xij)=P(ij∈V|Xij)本質上就是傾向得分。由于p(Xij)為一個概率，其范圍在0到1之間，取值有限，直接建立多層模型不太合適，因此考慮使用聯結函數，如最常用的Logit聯結，即logit(p(Xij))=log(p(Xij)(1-p(Xij)))，將概率范圍(0，1)轉換為(-∞，∞)，再建立多層模型。在建立多層模型時，有很多的模型選擇，比如第一層次的預測變量可以假定在群層次（第二層次）上是固定的或變動的；第二層次的預測變量可以僅包含在截距方程中，僅影響截距，也可以既包含在截距方程中，也包含在斜率方程中，同時影響截距和斜率?？傊?，在模型假定上進行不同層次預測變量的添減以及截距、斜率隨機性的設置，就會產生許多不同的多層模型。下面分別進行討論。

（1）具有固定斜率隨機截距且不包含第二層次預測變量的傾向得分多層模型。

這里logit(p(Xij))表示傾向得分的Logit函數；β0為斜率方程中的斜率（固定）；β1為第一層次方程中的斜率（固定）。該模型假定不同群中X1與logit(p(Xij))回歸方程的斜率相同截距不同，截距不同是由群的差異所引起的，但并未考慮群層次上具體變量對截距的影響。此模型通過固定斜率，假定X1對logit(p(Xij))的影響在不同群中是一樣的。

（2）具有固定斜率隨機截距且包含第二層次預測變量的傾向得分多層模型。

其中Wj為第二層次上的預測變量，γ01為固定斜率。此模型假定第二層次預測變量W只對X1與logit(p(Xij))回歸方程的截距產生影響，換句話說不同群中X1與logit(p(Xij))回歸方程的截距不同斜率相同，截距受到預測變量W的影響。

（3）具有隨機斜率固定截距且不包含第二層次預測變量的傾向得分多層模型。

此模型假定不同群中X1與logit(p(Xij))回歸方程的截距相同斜率不同，斜率不同是由群的差異所引起的，但并未考慮群層次上具體變量對斜率的影響。

（4）具有隨機斜率固定截距且包含第二層次預測變量的傾向得分多層模型。

其中Wj為第二層次上的預測變量，γ11為固定斜率。此模型假定不同群中X1與logit(p(Xij))回歸方程的斜率不同截距相同，斜率受到預測變量W的影響，即第二層次預測變量W對第一層次預測變量與因變量之間的關系產生影響。

（5）具有隨機斜率隨機截距且均不包含第二層次預測變量的傾向得分多層模型。

此模型假定不同群中X1與logit(p(Xij))回歸方程的截距和斜率均不同，截距和斜率的不同是由群的差異所引起的，但并未考慮群層次上具體變量對截距和斜率的影響。

（6）具有隨機斜率與截距且均包含第二層次預測變量的傾向得分多層模型。

其中γ01為截距方程中預測變量固定的斜率，γ11為斜率方程中預測變量固定的斜率。此模型假定預測變量W對不同群中X1與logit(p(Xij))回歸方程的截距和斜率都產生影響。

（7）具有隨機斜率與截距且截距斜率其中之一包含第二層次預測變量的傾向得分多層模型。

此模型假定預測變量W對不同群中X1與logit(p(Xij))回歸方程的截距或斜率產生影響。以上只是考慮了第一層、第二層中一個預測變量的情形，可以根據具體的情況，添加多個預測變量。

在建立傾向得分多層模型之后，可以估計出logit(p(Xij))，得到傾向得分的估計帶入式（4）則得到 πij的估計值，即：

其中網絡候選者數據庫的調查樣本SV可分成H個群（第二層次），第j個群包含mj個個體（第一層次），j=1，…，H。如果總體規模N已知時，可將式（13）的分母換為N。該加權調整方法可稱為傾向得分逆加權方法。此外，對總體均值估計的方差可使用重抽樣方法，如Bootstrap、Jackknife等來進行估計。

3 模擬研究

現對本文所提出的總體估計方法的效果進行模擬研究。首先生成一個大小為100000（N=100000）的具有嵌套數據結構的有限總體，協變量（預測變量）包括個體層次（第一層次）上服從正態分布N(1，1)的變量X以及群層次（第二層次）上服從正態分布N(0，1)的變量W，群個數m為50，每個群包含的個體單元數都相等為Nm=2000，目標變量Y值由兩個協變量X、W以及多層線性模型Y=-0.5+0.5X-0.4W+0.3W×X+u1X+u0+ε來生成，其中u0、u1與ε分別由以下分布來產生：

在生成的有限總體中不放回簡單隨機抽取20個群（每個群2000個單元），然后在每個群中不放回簡單隨機抽取50個單元，得到樣本量為1000（nr=1000）的一個樣本，視其為參考樣本。然后，從20個群共40000個單元中去掉參考樣本單元后的剩余單元中，再分別從20個群中不放回簡單隨機抽取1000個單元，得到一個有20個群且每個群包含1000個單元的嵌套樣本，將其作為網絡候選者數據庫V，規模為M=20000。接著，采取分層隨機抽樣從網絡候選者數據庫V的20個群中不放回簡單隨機抽取500個單元，作為網絡候選者數據庫的調查樣本SV，樣本量為10000（nv=10000）。將上面的過程重復1000次，獲得1000組蒙特卡羅樣本，針對每一組樣本構建具有隨機斜率與截距且均包含第二層次預測變量的傾向得分多層模型以及單一層次結構下的Logistic傾向得分模型（忽略數據的嵌套結構）來估計傾向得分，并采用傾向得分逆加權方法來估計總體均值。這里的Logistic傾向得分模型為：

其中Xij=(Xij1，Xij2，…，Xijp)′為協向量，p為協變量個數，β為回歸參數。最后，計算1000組模擬樣本上總體估計的均值、方差、偏差與均方誤差。

表1 基于傾向得分多層模型的總體均值估計的模擬結果

表1為基于傾向得分多層模型和Logistic模型并使用傾向得分逆加權調整方法計算的總體估計的均值、方差、偏差與均方誤差。兩種方法估計的總體均值較為接近，且總體均值估計的方差基本相同。從偏差上看傾向得分多層模型的逆加權估計偏差相對來說更小，而兩種方法估計的總體均值的均方誤差一樣，說明對于具有嵌套數據結構的網絡候選者數據庫，基于傾向得分多層模型的總體估計相對于基于傾向得分Logistic模型的總體估計偏差更小，無偏程度更高，而兩個估計的效率基本一樣。

4 實證分析

實證數據來自2014年美國綜合社會調查（General Social Survey，簡記為GSS），選取調查中的性別（gender）、年齡（age）、種族（race）、區域（region）、教育（education）、收入（income）和投票變量（vote12）共7個變量，2538個個案的一個數據集，其中區域變量為群的分類變量，共分9個群。因為數據中出現了“不知道（Don't know）”、“無回答（No answer）”、“不適用（Not applicable）”、“不合格（Ineligible）”等選擇結果，所以本文將這些結果的個案刪除。此外，為了簡化計算進一步對變量值進行了一些變換，獲得7個變量、191個單元的最終數據集，7個變量的類型與具體取值見表2。

表2 基于傾向得分多層模型推斷的實證變量

由于多層模型需要較多的數據，因此本文首先利用Bootstrap來構造一個偽總體。具體地，使用有放回簡單隨機抽樣方法，從191個單元中抽取一個大小為100的樣本，重復抽取100次，獲得10000個單元，將這10000個單元和191個單元匯總，從而形成一個大小為10191的偽總體（N=10191）。然后，先使用PPS抽樣從偽總體中抽取7個群即7個區域，再在每個區域中不放回簡單隨機抽取20個單元，得到一個大小為140的參考樣本（概率樣本，nr=140）。接著，從7個區域中去掉參考樣本單元之后，再分別從每個區域不放回簡單隨機抽取300個單元作為網絡候選者數據庫V（規模為M=2100）。進一步從網絡候選者數據庫的每個區域中不放回簡單隨機抽取一個樣本量為100的樣本，即網絡候選者數據庫的調查樣本SV（nv=700）。在偽總體固定的情況下，重復抽取參考樣本和網絡候選者數據庫的調查樣本1000次。

現對2012年總統選舉中投票人數所占的比例，即投票變量Yvote的均值進行估計，將年齡、性別、種族、教育、收入作為投票個體（第一層次）上的協變量，并且投票的個體嵌套于區域層次（第二層次）之中。在1000組參考樣本與網絡候選者數據庫的調查樣本上，建立具有固定斜率隨機截距且不包含第二層次預測變量的傾向得分多層模型以及Logistic傾向得分模型來估計傾向得分，并將傾向得分的逆作為權數進行加權調整來估計投票人數所占的比例。最后，計算1000組樣本上總體估計的均值、方差、偏差與均方誤差，見表3。

表3 基于傾向得分多層模型的總體均值估計的實證結果

由表3可知，建立傾向得分多層模型與Logistic模型并采用傾向得分逆加權方法計算的總體均值估計的均值都約為0.613，非常接近，并且兩種方法估計的總體均值的方差也相差較小。在偏差與均方誤差上，傾向得分多層模型逆加權估計的偏差絕對值比傾向得分Logistic模型逆加權估計的偏差絕對值要小，而且傾向得分多層模型逆加權估計的均方誤差也較小，表明從偏差和均方誤差兩個方面來看，對于嵌套結構的數據，傾向得分多層模型逆加權估計不僅無偏程度較高，而且較有效率，估計比較穩健。

5 結論

本文針對網絡候選者數據庫具有嵌套結構數據的情況，提出建立傾向得分多層模型來估計傾向得分，并利用傾向得分逆加權方法來對非概率樣本進行加權調整，從而推斷總體。模擬與實證研究表明，對于具有嵌套結構的數據，相對于基于傾向得分Logistic模型的逆加權估計來說，基于傾向得分多層模型的逆加權估計偏差更小，估計的效率更高。可見，傾向得分多層模型比傾向得分Logistic模型更適合具有嵌套結構數據的網絡候選者數據庫。