基于混合效應模型的群體數目估計

白永娟，李好奇

（長江師范學院 數學與統計學院，重慶 涪陵 408100）

0 引言

捕獲再捕獲研究是一種估算群體數目的方法，可以用來更好地了解潛在群體數目的動態變化。所討論的群體，根據是否存在出生、死亡、遷入、遷出可以分為開放群體和封閉群體。本文主要討論有多個觀測機構的開放群體數目估計問題。對于開放群體數據，有一些特征需要注意。第一，數據是捕獲再捕獲數據，屬于有偏抽樣，即僅僅被捕獲到至少一次的個體被觀測；第二，開放群體的個體數目是隨時間變化的，即不同時間的群體數目在不斷變化；第三，每個個體被捕獲概率是不同的，還要考慮個體協變量對捕獲概率的影響，以及無法觀測的個體異質性存在。

對于封閉群體數目估計，有很多文獻提出估計群體數目的方法。比如泊松對數線性模型(Poisson log-linear model)[1-3]、多項式模型[4]、樣本覆蓋方法[5]。對于開放群體，也有一些方法被提出來，如文獻[6-9]。這些方法都沒有擴展到多重列表問題。最近Lin等[10]提出了半參數方法來估計開放群體多重列表問題，但是沒有考慮協變量特征如性別、年齡等對捕獲概率的影響。考慮協變量特征的開放群體多重列表估計問題很少有文獻進行研究。

本文提出廣義混合線性回歸模型來估計多重列表的開放群體數目，同時考慮個體協變量特征對捕獲概率的影響。所提出的模型允許不可觀測個體異質性存在。由于捕獲再捕獲數據的有偏抽樣屬性，使得廣義混合線性回歸模型的標準估計方法不可用，本文提出基于條件似然的估計方法，可以得到相關參數的極大似然估計，進而估計出群體數目。得到的估計量都將證明相合性，漸進正態性。

1 模型和估計方法

1.1 符號及模型

把整個捕獲時間劃分為等長度的小時間區間t=1，2，…，T，假設有d個捕獲機構。在每個時間區間t,有nt個個體至少被捕獲一次，記錄詳細的個體特征和對應捕獲機構。令ytij表示個體i被機構j在時間t被捕獲的示性函數，被捕獲取值為1，否則為0，Xti表示對應的協變量。令yti=(yti1，…，ytid)′和則觀測數據 (Xti，yti)僅在δti=1的時候被觀測。在時間區間t內的群體個數表示為νt。本文目的是對任意給定時間段t，基于觀測數據(Xti，yti)，估計出未知群體數目νt。假設ptij表示在時間段t內個體i被機構j捕獲的概率，考慮下面的模型：

其中j=1，…，d和i=1，…，νt,βtj反映的是機構j隨時間改變的捕獲能力，因為群體數目會隨時間改變。ai是隨機效應，反映對象的特殊響應趨勢，例如基于已知協變量特征判斷個體的被捕獲概率很低，但由于隨機效應存在，實際被捕獲概率很高。此外個體對多個機構響應的相關性可以通過隨機效應ai來表示。本文假設ai是均值為零方差為σ2的正態隨機變量。

1.2 模型估計

本文給出全似然函數[11]:

其中f(nt)表示從νt個體中捕獲到nt個個體的二項概率，f(Xti|δti=1)是Xti的條件密度函數，f(yti|Xti，δti=1)是yti的條件密度函數，則：

其中pt表示在第t個時間段內個體平均被捕獲的概率。用ft(·)表示Xti的密度函數。可以證明f(Xti|δti=1)=f(δti=1|Xti)ft(Xti)/pt，因此：

令qti表示概率f(Xti|δti=1)，則根據式

（4）可以得到：

現在考慮f(δti=1|Xti)，表示至少被捕獲一次的概率，可以被寫為：

其中ptij(x，a)是ptij中Xti，ai分別用x和a代替。從式（5）和式（6）可以看到pt是βtj，αj，σ2和qti的函數。對于yti的條件密度函數有：

把式（3）、式（5）至式（7）代入式（2），可以得到對數似然函數：

計算 log{L(β，α，σ2，ν)} 關于αd，t=1，…，T，i=1，…，nt的導數并令導數為0，即可得到得分方程：

其中λt是拉普拉斯乘子，qti具有限制條件而的展開形式是ν的函數，在[n，∞)上是凹函數，

tt具 有 連 續 二 階 導 數 ，在νt=nt/p?t處 的 一 階 導 數 為-log{1-p?t},其中p?t是pt的估計量。

其中：

1.3 漸近性質

討論n=mtin{nt}趨于無窮的時候，本文給出所提估計量的漸進分布。符號 →d表示“依分布收斂”。求出lN(θ)關于θ的導數，得到得分函數：

假定θ?=(β?，α?，σ?2)是得分方程U(θ)=0 的解。進一步，通過泰勒展開可以得到：

利用參數模型中極大似然估計量標準漸近理論[12]，在正則條件下：

其中I(θ)是參數θ的費希爾信息矩陣。根據delta方法，可以得到：

其中：

給定條件X1，…，Xnt，逼近式(12)的第二部分均值為0，第一部分對于X1，…，Xnt是可測的，且條件均值為0，最后一部分對于隨機變量nt是可測的，均值也是0。

其中bt定義如方程 (13)。 另外p?t=nt/ν?t，利用類似的方法可以得到：

1.4 方差估計

其中是1/π(Xti;θ0)的樣本方差其中是的樣本均值是pt的估計量。

這個方法在B≥100的時候效果較好。

2 模擬

利用數值例子來說明本文方法的效果。兩維協變量Xti不隨時間變化。Xti第一個成分服從標準正態分布，獨立于第二成分。第二部分以相等概率取值1和0。每一種設置進行500次重復模擬。

情形1：設置為T=5，d=4，每一期的群體個數分別為ν=200；情形2：設置為T=5，d=4，每一期的群體個數分別為ν=1000；情形3：設置為T=5，d=8，每一期的群體個數分別為ν=200。

表1給出了情形1下所提方法的結果，包括基于500次重復計算的偏差，標準差。從表1可以看出，在樣本量較小的情形下，估計結果效果良好，能很好地估計出群體數目。情形2相對于情形1，樣本量增加了，其他設置保持不變，從基于500次重復計算的結果來看，所提方法依舊有效。由于同情形1類似，故結果未列出。表2給出了情形3下所提方法的結果，包括基于500次重復計算的偏差，標準差。情形3是對于捕獲列表增加的情況，考察所提方法的效果。從表3展示的結果看，捕獲列表較多情形下，所提方法仍然效果較好。

表1 設置T=5,d=4,v=200下模擬結果

表2 設置T=5,d=8,v=200下的模擬結果

3 結論

對于多列表捕獲再捕獲問題，本文通過混合效應模型，對原始捕獲數據進行分析，得到每個時間段群體個數相對客觀的估計。同時，利用隨機效應評估每個個體對捕獲機構的反應。個體隨機效應彌補了個體協變量不能描述的個體反映。

捕獲再捕獲數據在多維列表情況下，為了分析的簡單，一般設定各個機構的捕獲是獨立進行的，即假設各個列表獨立。這個假設在一些情況下可以放松，假設各個捕獲列表之間具有相關性，這個可以作為下一階段研究的內容。

參考文獻：

[1]Fienberg S E.The Multiple Recapture Census for Closed Population and Incomplete 2k Contingency Tables[J].Biometrika,1975，(59).

[2]Cormack R M.Log-linear Models for Capture-recapture[J].Biomet?rics,1989,(45).

[3]International Working Group for Disease Monitoring and Forecasting.Capture Recapture and Multiple-Record Systems Estimation.I:Histo?ry and Theoretical Development[J].Am.J.Epidemiol,1995,(142).

[4]Cormack R M,Jupp P E.Inference for Poisson and Multinomial Mod?els for Capture-Recapture Experiments[J].Biometrika,1991,(78).

[5]Chao A,Lee S M.Estimating the Number of Classes via Sample Cover?age[J].J.Amer.Statist.Assoc,1992,(87).

[6]Huggins R M,Yip P S F.Estimation of the Size of an Open Population From Capture-Recapture Data Using Weighted Martingale Methods[J].Biometrics,1999,(55).

[7]Huggins R M,Yang H C,Chao A.Population Size Estimation Using Local Sample Coverage for Open Populations[J].J.Statist.Plann.Infer?ence,2003，(113).

[8]Yang H C,Huggins R M.The Estimation of the Size of the Open Popu?lation Using Local Estimating Equations[J].Statist Sinica,2003,(13).

[9]Yang H C,Huggins R M,Clark A S S.Estimation of the Size of an Open Population Using Local Estimating Equations II:A Partially Parametric Approach[J].Biometrics,2003,(59).

[10]Lin H,Yip P S,Chen F.Estimating the Population Size for a Multi?ple List Problem With an Open Population[J].Statistica Sinica,2009,(19).

[11]Chen K.Parametric and Semiparametric Models for Recapture and Removal Studies:A Likelihood Approach[J].J.R.Statist.Soc.B,2001,(63).

[12]Van der Vaart A W.Asymptotic Statistics[M].Cambridge:Cambridge University Press,1998.