使用数据增强分析具有单个协变量的捕获-再捕获模型

总结

我考虑使用影响检测概率的单个协变量分析捕获-再捕获模型。联合似然的贝叶斯分析是使用一种灵活的数据增强方案进行的，该方案有助于使用马尔可夫链蒙特卡罗方法进行分析，并在免费软件中简单明了地实现。该方法用于草地田鼠的研究(草甸田鼠)其中记录了连续协变量（体重）的辅助数据，并认为检测概率与体重有关。在第二个示例中，该模型应用于使用双观测器协议的空中水禽调查。观测的基本单位是个体鸟类的集群，集群的大小（离散协变量）被用作检测概率的协变量。

1.简介

在动物种群研究中，个体间具有异质检测概率的捕获-再捕获模型被广泛使用。一类异质性模型是将单个检测概率视为由某些概率分布描述的潜在参数。近年来，此类模型受到了相当大的关注，已经开发出许多特定模型，这些模型对应于单个检测概率参数的替代概率分布(诺里斯和波洛克，1996年;库尔和阿格雷斯蒂，1999年;Dorazio和Royle，2003年;Royle和Nichols，2003年;Royle，2006年). 这些模型的解释和实用性受到了以下方面的质疑链接（2003）他表示，为了最实际的目的，N个在混合分布类别中无法识别。缓解这一问题的一种可能方法是尝试明确建模异质性，例如，识别被认为可以解释检测概率变化的单个协变量。明显的协变量包括年龄、性别、大小或群体大小。

捕获-再捕获模型中的单个协变量问题已经引起了相当大的关注（参见波洛克，2002年;威廉姆斯、尼科尔斯和康罗伊，2002年，第300-302页，供审查）。已经提出了三种基本策略来对个体协变量进行建模。建议的一种方法波洛克、海因斯和尼科尔斯（1984）需要将个体分层为有限数量的离散类，从而产生K（K）具有地层人口规模的地层{N个_k}^K（K）_k=1。在这种方法下N个_k参数是推理的对象。这种方法适用于离散协变量，尤其是那些具有少量值的协变量。一个缺点是，参数维数随协变量的基数或连续协变量的层数而增加。通过指定较少的地层来避免维度增加，会导致协变量效应的近似值较差，并且随着额外协变量的考虑，问题变得更加复杂。第二种广泛使用的策略是Huggins–Alho程序(哈金斯，1989年;Alho，1990年). 这是基于从所谓的“条件似然”（捕获数据的条件似然）导出的单个可检测性估计值。在这种方法下，N个是一个导出的参数，其估计基于广义Horvitz–Thompson估计（HTE），并采用了经典的渐近推断方法。第三种方法，即所谓的“联合似然”方法，规定了遭遇数据和协变量的联合分布x个例如，通过为协变量规定概率分布克(x个). 这种方法与经典广义线性模型中缺失协变量的当代处理方法一致。请参见Ibrahim等人（2005年）进行全面审查。波洛克（2002）提供了不同方法的良好概述。

在本文中，我考虑了个体协变量模型的贝叶斯分析。分析基于遭遇历史观测值和个体协变量的联合分布。马尔可夫链蒙特卡罗（MCMC）的贝叶斯分析是通过对未知样本量的多项式模型使用有效的数据增强过程来辅助的(Royle、Dorazio和Link，2007年). 数据增强产生了一种灵活的模型重新参数化，允许进行简单的贝叶斯分析。提供了两个示例。首先，对草地田鼠进行活体诱捕研究，其中体重被认为可以解释探测概率的异质性。第二个例子涉及水禽调查数据，其中评估了离散协变量（集群大小）的影响。

2.联合可能性和数据增强

假设一个封闭的种群N个采样时间：T型场合，生成n个独特的个人，其中个人我 = 1, 2, … , n个被捕获{年_我}^n个_我=1次。为了清楚起见，我们假设检测概率在T型但以下发展直接扩展到捕获历史方面的公式（中提供了一个示例第4节). 假设每个人身上都测量了一个辅助变量，比如x个_我这一变量被认为会影响个体的可检测性。我们假设捕获是具有参数的独立同分布（i.i.d.）伯努利试验第页(x个_我; θ₁) ≡第页_我例如，

其中，在本例中，参数θ₁是向量θ₁= (α₀, α₁).

我现在要解决的问题是估算人口总数，N个，其中n个对观察到的个体进行采样，并将模型参数包括方程（1）中检测函数的参数。这样做的主要困难在于样本是一个有偏见的样本，即样本中出现的个体往往具有较高的第页_我因此，个体出现在样本中的概率也是未知模型参数的函数。与丢失数据的传统处理类似，我们通过指定协变量的概率分布来进行处理x个，说吧克(x个 | θ₂). 对于许多连续的协变量，假设克为平均值μ和标准偏差σ的正态分布。符号θ₂将表示的参数克不管他们的人数。

为了在该模型下进行推断，我开发了以下观测值的可能性年和x个从联合分配年和x个以个人出现在样本中的事件为条件，即年> 0. 抑制对参数的依赖性年和x个是

这里我们使用了传统的“括号表示法”，例如[u个]表示的概率密度函数（pdf）u个, [u个 | v（v）]表示的条件pdfu个鉴于v（v）等等。观测值年是整数k = 1, 2, … , T型因此[年,x个 | 年 = 0]对于所观察到的相遇频率而言是相同的0。因此，条件似然为

假设捕获是i.i.d.Bernoulli试验，因此[年 | x个]是带索引的二项式T型，则捕获的边际概率为

哪里θ= (θ₁, θ₂)，纸盒（0；T型,第页(x个; θ₁))表示带索引的二项式pdfT型和参数第页(x个)根据方程式（1）。因此，条件似然的具体形式是

可以使用传统方法将其最大化，以获得θ的最大似然估计（MLE）₁和θ₂。MLE可用于获得所谓的“条件MLE”N个即。，⁠.

无条件似然是方程（3）和n个鉴于N个，它是带索引的二项式N个和参数1−（f）(0;θ). 这就产生了

有条件和无条件MLEN个是渐近等价的(萨纳塔南，1972年)，两者都被广泛应用于不同的人口规模估计环境中。

2.1通过数据增强进行贝叶斯估计

对上一节中描述的无条件或条件可能性的分析相对简单，主要困难在于方程（2）中积分的计算。在大多数情况下，这将以数字方式完成。对于具有第页单个协变量，方程式（2）替换为第页-维度模拟。随后，通过采用个体协变量模型的贝叶斯公式，避免了评估似然所需的集成，在该模型下，使用MCMC的传统方法表征后验分布。对无条件似然的分析引入了一个额外的复杂性，即参数空间的维数本身就是一个参数（即未观察到的数量x个的，“零观测值”的数量等于N个−n个,N个未知）。这个问题激发了以下描述的数据增强方案Royle等人（2007）用于在指数未知的多项式模型中建模和估计，N个在随后的分析中，利用此数据增强方案来固定部分观测变量的维数x个和年.

基本策略是用一个固定的已知数字来增加观察到的数据集M（M）−n个，对所有零捕获历史进行建模，并对增强数据集（大小M（M）)作为使用未知但可估计的零膨胀参数的完整数据模型的零膨胀版本。本质上，估算人口规模的经典多项式模型被重新参数化为如下形式的场地占用模型：MacKenzie等人（2002年）事实上，封闭人口规模估算模型和场地占用模型之间的这种二重性在很大程度上激发了数据增强的想法，但也可以基于贝叶斯论证进行更正式的论证(Royle等人，2007年). 具体而言，他们指出，数据增强相当于为N个关于整数到M（M）因此，对于M（M）≫N个，这为N个这通常被认为是非信息性的，或表达了关于N个只要M（M）足够大。贝叶斯构造为选择M（M）-它应该足够大，以便对后文摘要没有影响，否则它的使用（作为一种信息性的前文）应该像惯例一样需要一些理由。采用数据增强来分析单个协变量模型的主要好处是，它使用MCMC方法生成了简单高效的贝叶斯实现。

对于包含单个协变量的模型，数据包括n个观察到的检测频率年₁,年₂, … , 年_n个和相关协变量值x个₁,x个₂, … , x个_n个对于通过数据增强进行的估计，引入零伪操作(年_n个+1= 0,年_n个+2= 0, … ,年_M（M）=0），缺少协变量值{x个_我}^M（M）_我=n个+1和一组潜在指标变量{z（z）_我}^M（M）_我=1观察到的(z（z）_我=1）对于我 = 1, 2, … , n个并且未被观察到我 = n个+ 1, … , M（M）那么，假设z（z）_我～伯努利（ψ）。参数ψ可以解释为增加总体的一个元素（大小M（M）)是抽样人口中的一员N个)推断的对象是人口数量，⁠.

扩充数据的模型由以下三个组件组成：

• (1)

  z（z）_我～伯努利（ψ）；

• (2)

  [年_我 | 第页(x个_我)]=二项式(T型,z（z）_我第页(x个_我))使用logit(第页_我) =α₀+α₁x个_我;

• (3)

  [x个_我]=正常（μ_x个, σ²_x个).

使用传统的MCMC方法进行估计和推断很简单。以下各节给出了两个示例。MCMC算法要求对每个“缺失”进行采样x个_我从其完全条件分布来看，即[x个_我 | 年_我= 0]∝{(1 −第页(x个_我))^T型[x个_我]}. 可以使用Metropolis–Hastings算法等通用方法对此进行采样。否则，MCMC算法很容易实现。使用现有软件可以避免这些细节。随后的分析在WinBUGS软件中进行(Gilks、Thomas和Spiegelhalter，1994年). 在随后的分析中提供了此处描述的基本协变量模型的一些扩展。

3.应用：连续协变量

在这个例子中，草地田鼠的种群(草甸田鼠)1981年，在马里兰州的劳雷尔连续五天进行了采样。采样是使用以玉米为诱饵的10×10网格陷阱进行的。这里考虑的数据是56名成年男性的遭遇史，从6月27日开始连续五天的五个样本中获得。在第一次捕获时测量了每个个体的协变量体重（克），有一些先验的观点认为体重对检测概率有积极影响。在接下来的分析中，将体重协变量标准化为平均值为零和单位方差。有关该研究的更多详细信息，请参阅尼科尔斯、波洛克和海恩斯（1984）数据被广泛使用Williams等人（2002年）其他一些作者也进行了分析，包括沃尔特（1990）;Lee、Lee和Gee（2003）;杨和超（2005）; 和张、刘、游（2005）.

对于微型客车数据，遭遇频率分布为n个_k=（12，8，9，12，15）对于k = 分别为1、2、3、4、5。根据简单零模型（“M（M）₀“）是⁠，未捕获个体数的MLE发生在参数空间的边界上n个₀= 0. 考虑到该采样设备（网格中布置的诱饵陷阱）的高检测概率，这是一个合理的结果。个人至少被捕获一次的概率为1−（1−0.60）⁵= 0.99.

我将模型扩展到包括个体协变量“体重”。我们假设小鼠群体中的体重正态分布，平均μ_米和方差σ²_米.用M（M）−n个 = 100次观察年 = 0，扩展数据的模型如中所述第2节。在125和50个零的基础上增加分析次数，结果与蒙特卡罗误差一致。为了进行分析，需要结构参数的先验分布，为此，我采用了传统的默认先验值，表面上，这些先验值表示的模型参数先验信息很少。对于平均参数，使用平均值为0且方差为1000的正态先验值，而对于方差参数的倒数（τ=1/σ²)使用形状和比例均等于0.001的伽玛先验值。

MCMC算法的实现很简单。在这里，我将重点介绍WinBUGS软件的实现。这只需要模型结构的“伪代码”表示。也就是说，数据、潜在变量和参数的概率分布规范，以及参数之间的任何函数链接。WinBUGS模型规范Microtus公司Web附录A中提供了数据，附录中给出了从R中执行WinBUGS的R脚本。

后向分布N个基于50000个蒙特卡洛从后验分布中提取的数据，总结如下表1。基于估计分布的95%后向间隔为[56，69]。体重对检测概率的影响似乎非常显著，α₁高于0。模型参数的各种其他后验总结见表1.

3.1 Huggins–Alho程序

Huggins–Alho程序(哈金斯，1989年;Alho，1990年)是分析单个协变量模型的事实标准。这是一个启发式程序，通过与HTE的类比证明是合理的。我将此估计值应用于Microtus公司数据集，并生成⁠，这与完全基于模型的结果一致。将Huggins–Alho程序扩展到更一般的情况需要制定模型，其中N个通过调节移除n个.

3.2模型扩展

数据增强框架可以直接扩展，以适应更通用的模型。例如Microtus公司诱捕器是以玉米为诱饵的，这就增加了通常被称为行为反应的可能性，在这种情况下，我们可能期望在初始捕获后检测概率会增加。我在这里考虑经典的行为反应模型M（M）_b条以及建议的扩展杨和超（2005）允许短期行为反应（他们称之为短暂反应）。这类更广泛的模型具有以下形式

协变量x个_1,ij公司是上一次捕获的指示器。也就是说，x个_1,ij公司=1，如果在采样前捕获个体j个该模型包含类似马尔可夫自回归的结构年_我,j个−1这是对被俘的短期短暂反应。

在WinBUGS中描述模型相当简单（请参阅Web附录A）。结果如所示表2，我们发现似乎没有非常强烈的行为反应。有一些证据表明，这是一种轻微的短暂反应，从这个意义上说，对于该参数，大多数后部质量都在0以上。

我考虑进一步扩展，以包括体重协变量：

协变量x个_1,ij公司是以前捕获的指示器，如前所述，并且x个_2,我是体重协变量。我假设小鼠群体的体重是平均μ的正态分布_x个和方差σ²_x个即。，x个_2,我～正常（μ_x个, σ²_x个). 与之前的模型一样，MCMC对模型的增广数据分析非常简单。结果总结如下表3，其中我们看到（第2列和第3列）个体协变量体重的巨大影响，并且在存在体重协变量的情况下，持续的行为反应非常强烈。另一方面，短暂的反应很弱。因此，删除了模型的这个组件，并重新安装了得到的更简单的模型（中的第4-8列表3). 我们看到两个α的后验分布₁（行为效应）和α_三（体重系数）集中在0以上。因此，结果表明，一旦捕获个体，检测概率会大幅增加，而且体重较大的个体捕获概率也较高。还应注意，当模型中包含这两种效应时N个确实从零开始偏移。这是一种直观的效果，因为较高的明显发现概率是由于先前捕获的个体的捕获率较高，以及较重个体被捕获的倾向。

体重协变量的估计总体平均值和标准差表明，在观察到的个体上测量的协变量只代表一个稍微有偏差的样本。而体重的样本平均值和标准偏差分别为41.79和11.96E类[x个]和Var[x个]（在回传中给出的μ和σ估计值的后验平均值后表3)通过求解获得(x个−41.79）/11.96=1.122，得出总体值E类[x个]=39.45和Var((x个− 41.79)/11.96) = 1.122²得出SD的总体标准偏差(x个) = 13.42. 因此，如预期的那样，体重的人口分布比样本略低，变化更大。个人被捕获的概率秒对平均体重的标准偏差进行了评估秒 = （−2，−1，0，1，2），对于之前未捕获的假设个体，结果为（0.532，0.848，0.984，0.999，1.000）。体重低于平均值两个标准差的个体在以下研究中被捕获的可能性只有50%左右：T型 = 5个周期。因此第页是由第页关于协变量的体重。现在，考虑对一个个体进行相同的计算是之前捕获的。这将产生（0.766、0.964、0.998、0.999、1.000）。对此的解释如下：初始捕获后，低于平均体重的单个两个标准差的捕获概率大约加倍。由于平均捕获概率高（α₀)对于体重远高于平均值的个体来说，初始捕获几乎没有影响。

3.3选型

之前的结果似乎相当明确，支持体重的影响，以及在存在体重的情况下“永久”的行为反应。在这里，我解决了基于后验模型概率的形式化模型选择问题。每个模型的后验概率都是通过指定一组潜在指标变量来计算的，例如，每个模型效应对应一个变量w个_j个对于j个第个效应，并在每个效应之前施加一个伯努利w个_j个，比如有参数π_j个。这个概念是由Kuo和Mallick（1998）（另请参见康登[2005年，第3.2节]）。潜在指标变量包括：

具有先验分布：

这里我假设它们是相互独立的。通过将线性预测器指定为

在数据增强框架内对该扩展模型进行分析不会带来任何额外的困难。一个更广泛的技术考虑是，众所周知，后验模型概率对参数的先验值很敏感（例如。，艾特金，1991年;Link and Barker，2006年)在评估后验模型概率时，模糊先验通常不是无害或“无信息”的。Akaike信息标准的使用似乎在很大程度上避免了这个问题，因为它强加了一个特定的、可能是不希望的先验，这是样本量的函数(卡丹和拉扎尔，2004年). 这个问题的一个解决方案是，在一个模型下计算后验模型概率，其中参数的先验值固定在完整模型下的后验分布上(艾特金，1991年). 在本例中，先验值是回归系数四维后验分布的正态近似值（从中总结出表3派生）。使用此先验分布，从100000次燃烧后蒙特卡罗绘图中获得后验模型概率。

八种可能的模型由二进制序列索引(w个₁,w个₂,w个_三)，其中w个₁,w个₂、和w个_三表示持续行为反应、短暂反应和体重效应。模型（1,1,1）、（1,0,1）、（0,1,1）和（0,0.1）的后验概率分别为0.272、0.347、0.279和0.102。其他四个模型的后验概率为0。我们看到，基本模式倾向于具有体重的模型（所有四个具有正概率的模型都具有体重）。然而，后验概率在包含体重和一个或两个行为反应的三个模型中大致均匀分布。

3.4模型灵敏度

根据联合分布制定个体协变量模型年和x个生成一个明显基于模型的过程。具体而言，需要对检测概率和协变量之间关系的函数形式以及协变量的总体分布进行假设。后者不同于基于条件似然的经典方法，因此很自然会质疑这个额外的假设是否重要，因为结果对选择克(x个).链接（2003）在允许个体异质性的相关模型类别中（所谓的“模型M（M）_小时“），异质性分布的不同参数选择产生了关于N个人们自然会质疑这种模糊性是否会出现在当前的模型类别中。注意，在所考虑的模型中链接（2003）,第页被视为具有分布的个体特定的潜在参数克(第页). 在目前的模型类别中，我们有小时(第页_我) =α₀+α₁x个_我.所考虑的模型类别之间的主要差异链接（2003）这里考虑的单个协变量模型是，在后者中，我们观察到一些关于克(x个)凭借观察x个_我针对样本中出现的个人。与链接（2003）.

对选择的敏感性克在完整模型（包括行为效应和体重协变量）下对以下几种选择进行评估克包括前面考虑的正态、双指数、logistic、log-normal和t吨-具有两个自由度的分布。对于感兴趣的参数，在每个模型下的结果，N个，在Web中进行了总结表1给出了后验平均偏差。这些结果确实表明了N个选择克(x个). 然而，对后验平均值的影响N个个体的整体可检测性较高，在一定程度上减轻了这种影响。体重系数似乎对选择克.

4.示例：野生动物调查中的集群大小

在这个例子中，我考虑了绿头鸭的航测数据(Anas platyrhynchos公司)美国鱼类和野生动物管理局在2005年美国东北部和加拿大东部年度水禽种群调查期间收集的数据(Koneff等人，2008年). 固定翼测量中的样本单位为18英里直线段(史密斯，1995). 采用双观察者抽样方案，由两名观察者（前排座椅、后排座椅）进行抽样(马格努森、考利和格里格，1978年). 这产生了三个形式为{（0，1），（1，0），（1,1）}的可观察到的遭遇历史，表明该星团是只被后座观测者看到，还是只被前座观测器看到，还是两个观测者都看到。此类数据的模型相当于一项捕获-再捕获研究T型 = 2个周期和检测概率第页₁（对于前排观察员），以及第页₂（后视观察员）。预期前座观测器的检测概率高于后座观测者。本次调查共观察到162群鸟类。检测历史频率见表4观察到的簇大小在1到7之间，很少观察到的簇大小大于2。在下文中，假设种群簇大小分布为1+泊松（λ），即右移泊松分布，平均值为1+λ。聚类检测概率模型我是

观察员k = 1, 2.

数据增加通过添加大量（0，0）遭遇历史进行。在这个分析中，添加了125个零遭遇历史，这证明在以下意义上是足够的：N个集中在远离上限的地方（参见图1). 要在WinBUGS中实现模型，请注意，每个观察都是一个样本大小为1的四维多项式随机变量的实现。多项式对应于遭遇历史{（0，0），（0，1），（1，0）和（1，1）}的单元概率为{（1−第页_我1)(1 −第页_我2), (1 −第页_我1)第页_我2,第页_我1(1 −第页_我2),第页_我1第页_我2}用于群集我由x个_我个人。该多项式模型的零膨胀类似于第3节WinBUGS模型规范见Web附录A。

感兴趣的数量是采样组的总数：⁠，以及总人数（回忆一下z（z）_我是一个潜在的指标，它指示着增加人口中的一个因素是否具有规模M（M）也是抽样人口的一个要素）。协变量值x个_我对于我 = n个+ 1,n个+ 2, … , M（M）是缺失值，这些值作为MCMC算法的一部分从其后验分布中采样。本研究的实际重点是获得测量区域的密度估计值，因此N个_印度需要按总样本面积进行缩放。就目前而言，我们忽略了这个细节。

基于从后验分布中提取的100000个蒙特卡罗模型参数的后验总结如下表5.簇大小系数β的后部质量集中在零附近。因此，集群大小似乎对检测概率没有影响。估计的平均检测概率按预期排序，第页₁=0.80和第页₂=0.51，分别用于前排和后排座椅观察者。最后给出了个体总数的后验分布，N个_印度，如所示图1和两者的后部总结N个_c（c）（集群数量）和N个_印度在中给出表5.

5.讨论

捕获-使用结构化异质性或“单个协变量”重新捕获模型，以解释许多问题中自然出现的检测概率变化。通过采用模型的贝叶斯公式，可以实现单个协变量模型的建模和推理的通用灵活框架。使用MCMC方法和数据增强方案可以实现简单的实现。

对于这里所考虑的模型，由数据增强引起的重新参数化产生了一个等效于站点占用模型的模型（例如。，MacKenzie等人，2002年)，一些位点缺少协变量值。由此产生的估计问题类似于广义线性模型中“不可忽略”的缺失问题。缺失是不可忽视的“……如果未能观察到一个值取决于本应观察到的值”(Ibrahim等人，2005年). 在当前背景下，协变量缺失的概率（即个体未被捕获）取决于通过检测概率的协变量值。

个体协变量模型的两种变体应用于动物调查数据。在Microtus公司例子，第页假设在样本中为常数，但在个体间因连续协变量而变化。在这种情况下，完整的数据可能性是零膨胀二项式。在使用双观测器采样协议的空中水禽调查示例中，第页假设样本（对应于观察者）会发生变化，并且作为对离散协变量的响应，集群的大小也会发生变化。在这种情况下，完整数据的似然是零膨胀多项式似然。

个体协变量模型的传统处理是基于条件似然的，在条件似然中估计包含概率，然后将其用于HTE中的种群规模。这种方法在概念上最没有吸引力的方面是N个（许多问题中的推理对象）是一个派生参数，明确表示为妨害参数的函数（通常不直接相关）。基于数据增强的模型的贝叶斯公式有几个优点，包括其通用扩展性。例如，如果没有记录一些被观察个体的协变量值，就不会有额外的并发症。其次，使用这种贝叶斯公式在小样本中进行推理是有效的。最后，实践者可以很容易地访问实现，因为重新参数化的模型可以直接在免费的软件WinBUGS中实现。

单个协变量模型的完全基于模型的公式与基于所谓的Huggins–Alho程序拟合此类模型的标准方法相反。该估计器的动机是与HTE进行类比，在实践中似乎受到青睐，主要是因为它是“基于设计的”，因此，人们可能希望它对模型假设（尤其是对单个协变量的分布选择）具有鲁棒性。然而，Little（2004年）对这一点进行了深入的讨论，指出HTE具有基于模型的合理性，并且HTE的性能可能对偏离该模型的情况非常敏感。因此，无论是否采用表面上无模型的程序，都应考虑对模型选择的敏感性。

6.补充资料

第3节和第4节中引用的WinBUGS型号规范以及第3.4节，可在“纸张信息”链接中找到，网址为生物计量学网站http://www.biometrics.tibs.org.

致谢

我感谢美国地质调查局Patuxent野生动物研究中心的J.D.Nichols博士允许我使用在连续协变量示例中分析的microus数据，感谢美国鱼类和野生动物局的Mark Koneff使用水禽调查数据。

工具书类