基于混合效應(yīng)模型的群體數(shù)目估計

白永娟，李好奇

（長江師范學(xué)院 數(shù)學(xué)與統(tǒng)計學(xué)院，重慶 涪陵 408100）

0 引言

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

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

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

1 模型和估計方法

1.1 符號及模型

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

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

1.2 模型估計

本文給出全似然函數(shù)[11]:

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

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

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

（4）可以得到：

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

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

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

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

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

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

其中：

1.3 漸近性質(zhì)

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

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

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

其中I(θ)是參數(shù)θ的費希爾信息矩陣。根據(jù)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 模擬

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

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

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

表1 設(shè)置T=5,d=4,v=200下模擬結(jié)果

表2 設(shè)置T=5,d=8,v=200下的模擬結(jié)果

3 結(jié)論

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

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

參考文獻：

[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.