跳到主要内容
访问密钥 NCBI主页 MyNCBI主页 主要内容 主导航
Sankhya B（2008）。作者手稿；PMC 2022年11月1日提供。
以最终编辑形式发布为：
Sankhya B（2008）。2021年11月；83(2): 374–396.
2019年7月16日在线发布。 数字对象标识：2007年10月17日/13571-019-00198-7
预防性维修识别码：PMC8673559号
美国国立卫生研究院：美国国家卫生研究院1534857
PMID：34924732

具有非随机丢失的多元极值二元数据的共享空间模型

赵晓月,林章（Lin Zhang）、和Dipankar Bandyopadhyay公司
作者信息 版权和许可信息 PMC免责声明

摘要

牙周病（PD）的临床研究和试验在受试者的不同牙齿位置收集了大量数据。然而，它们呈现出一些统计复杂性。当我们的重点是了解极端PD进展的程度时，在具有对称（logit）链接的广义线性混合模型框架下进行标准分析可能不合适，因为二进制分裂（极端疾病与否）可能高度扭曲。此外，PD进展通常被假设为与空间相关，即近端牙齿可能具有与远端牙齿相似的PD状态。此外，观察到大量的缺失数据，这种缺失是非随机的，因为它反映了受试者的牙周健康状况。在本文中，我们通过一个共享（空间）潜在因素模型来解决上述所有问题，其中潜在因素通过广义极值回归对极端二进制响应进行联合建模，通过概率回归对非随机缺失牙齿进行联合建模。我们的方法是贝叶斯方法，推理框架由吉布斯-哈密顿蒙特卡罗技术支持。通过对PD真实数据集的仿真研究和应用，我们证明了我们的模型在模型拟合方面的潜在优势，并获得了与不考虑上述复杂性的备选方案相比的精确参数估计。

关键词：广义极值，哈密顿蒙特卡罗，潜在变量，非随机缺失，牙周病，空间

1 介绍

牙周病（Periontal disease，PD）或牙龈疾病仍然是全球口腔健康的主要负担，尤其是在美国。如果不进行治疗，它会通过牙齿周围的渐进性骨质流失和松动来损害生活质量，最终导致牙齿缺失。据充分记载，全世界约有10%至15%的人口易患严重的泛发性牙周炎(美国儿科学会，2005年). 用于评估PD状态的最重要的生物标记物是临床附着水平（CAL），以mm为单位，使用手动探针在所有28颗牙齿（在没有缺牙的受试者中）的每颗牙齿的六个预先指定的牙齿位置测量，不包括第三磨牙。图1提供了齿数、齿型和齿侧位置的象形描述。美国牙周病学协会(阿米蒂奇，1999年)根据PD的严重程度对牙齿进行分类，CAL≥3 mm表示PD中度至重度进展。

保存图片、插图等的外部文件。对象名为nihms-1534857-f0001.jpg

口腔内探测位置的图形表示：牙齿编号由T1-T7表示，有各种牙齿类型，如门牙（T1和T2）、犬牙（T3）、前磨牙（T4和T5）、臼齿（T6和T7）和牙侧位置（颊侧或口侧，与舌侧或舌侧）

在这篇论文中，我们有兴趣在一个记录居住在南卡罗来纳州沿海岛屿的讲古拉语的非洲裔美国人的牙病状况的激励性数据集中量化牙齿水平的二元PD状态（中重度PD与否）(Fernandes等人，2009年)，他们是2型糖尿病患者（以下简称GAAD数据），并研究其与各种受试者水平和牙齿水平协变量的关系。任意（可用）牙齿的二进制响应是通过根据≥3 mm阈值对相应牙齿（在预先指定的六个位置测量）的CAL值平均值进行二分来获得的，因此1表示“中度到重度”牙周炎组，0否则(Wiebe和Putnins，2000年). 在这种数据设置下，人们可能会尝试拟合广义线性混合模型（GLMM）(McCulloch和Searle，2004年)例如具有随机斜率或截距的logistic/probit回归模型。然而，对数据的仔细观察揭示了有趣的挑战。首先，在可用牙齿中，“中度至重度PD”类别的反应概率远低于“无至轻度”类别（17%）。在这种二进制响应的高度不平衡分布下，假设对称链接函数（例如logistic或probit）来拟合GLMM可能不充分，因为链接错误可能会导致回归参数估计中估计概率和偏差的均方误差增加，和平均响应概率(Czado和Santner，1992年；Roy和Dey，2014年). 其次，GAAD数据中有相当大比例的缺失数据（约26%）。无可争辩的是，PD是导致牙齿缺失的主要原因之一，与牙齿缺失较少的受试者相比，缺失牙齿较多的受试对象的牙周健康状况预计会恶化。因此，可以假设这种缺失是非随机的或信息性的(Reich和Bandyopadhyay，2010年)二进制PD状态，并保证其适当的建模。最后，PD进展可能表现出（潜在的）空间参照，即假设相邻牙齿的疾病状态与远端牙齿的疾病状况相似。理想情况下，适当的统计建模应考虑到上述三个限制，以避免有偏差的参数估计(Reich和Bandyopadhyay，2010年)。

现在有大量文献研究不平衡二进制响应的链接函数指定错误。其中绝大多数考虑了潜在变量框架(Chen等人，1999年)使用各种灵活的链接函数，例如广义偏移t吨-链接(Kim等人，2008年)，广义极值（GEV）链接(Coles等人，2001年)在本文中，我们重点探索GEV链接函数来建模空间相关的二进制响应。GEV连接及其无约束形状参数已被确定为处理严重不平衡独立的非常灵活的选择(Wang和Dey，2010年)以及空间相关(Li等人，2016a)回应。广义歪斜-t吨链接(Kim等人，2008年)优雅，但对形状参数的约束δ，以便0<δ≤1，大大降低了模型提供的可能偏斜范围(Wang和Dey，2010年). GEV提案没有这个问题。注意，除了观察到的数据外，缺失过程还可以表现出空间参考（GAAD数据集与前牙（如门牙和犬齿）相比，缺失后牙（即磨牙）的数量更多），从而为牙周整体健康评估提供信息。我们的分层贝叶斯方法通过联合通过共享空间因子对观测（多元）二进制数据过程（指示PD状态）和缺失数据过程进行建模(Reich和Bandyopadhyay，2010年)测量牙周健康。这两个过程之间共享的潜在牙周健康因素（随机效应项）通过观察到的受试者和牙齿水平协变量随牙齿而异，并使用条件自回归或CAR先验在空间上平滑(Banerjee等人，2014年). 在这种随机效应的条件下，观测数据过程和缺失过程是独立的。借助于本文采用的贝叶斯估计技术，可以利用观测数据过程和缺失数据过程中的所有可用信息来估计局部放电状态。

然而，通过标准马尔可夫链蒙特卡罗（MCMC）技术进行贝叶斯计算的一个重要瓶颈是模型参数的闭式后验分布不可用。广受欢迎的Metropolis-Hastings（M-H）算法的效率(Carlin和Louis，2008年)依赖于提案分布的质量，以及它快速探索整个似然区域的能力。当随机游走的方差很小时，典型的随机游走M-H在探索整个目标可能性时可能会很慢。另一方面，相对较大的方差会导致更高的拒绝率，链会在某个点上卡住，从而浪费计算时间。由于空间相关潜在变量注入的复杂关联结构，以及GEV链路引入的相关后验分布，M-H算法在这里不再适用。我们通过使用内吉布斯哈密顿蒙特卡罗（HMC）采样来规避这一点(Neal，2011年)从而加快马尔可夫链的移动和收敛。与随机行走的大都会相比，哈密顿动力学允许MCMC链沿着轨道进行大幅跳跃，从而使链在可能区域内移动更快。我们还通过Watanabe-Akaike信息标准（以下简称WAIC）评估了我们提出的GEV链路模型对其竞争对手（使用对称链路）的性能(渡边，2010)偏差信息标准或DIC(Spiegelhalter等人，2002年). 仿真研究和实际GAAD数据集的应用都说明了我们的GEV提案相对于备选方案的优势。

本文的其余内容如下。第2节开发空间共同的倾斜二元结果和缺失过程的建模框架，而第3节提出了贝叶斯推理框架，并使用MCMC对先验值、后验密度的形式和相关计算细节进行了必要的选择。第4节说明了联合模型及其对GAAD数据的应用，而第5节将我们的模型在不同丢失策略下的有限样本性能与对称probit链路模型进行了比较。最后，第6节提出了一些总结性意见。

2 关节空间模型

本节首先简要介绍GEV链接功能，然后继续进行模型开发。

2.1. GEV链接功能

GEV分布(麦克法登，1978年；Wang和Dey，2010年)是一系列连续概率分布，其分布函数GEV(μ,σ,ξ)由以下人员提供：

F类(x个∣μ,σ,ξ)={经验[−{1+ξ(x个−μσ)}+−1ξ],ξ≠0经验{−经验(−x个−μσ)},ξ=0}
(1)

哪里μ∈R（右）,σ∈R（右）+和ξ∈R（右）分别是位置、比例和形状参数，以及x个+=最大值（0，x个). 无约束形状参数ξ控制GEV密度等级的形状和尾部行为，足够灵活，以适应二进制响应中的广泛偏斜，并且可以根据响应曲线的偏斜进行识别和估计。特定值和范围ξ导致该类下的各种分布。例如，ξ>0、=0和<0分别产生Fréchet、指数衰减Gumbel和反向Weibull分布。

2.2. 空间相关性

接下来，我们介绍我们的数据设置。让年= (年1,年2, …,年n个)'表示an个×T型二进制响应矩阵，其中年我(t吨)表示齿的二进制CAL响应t吨,t吨= 1, …,T型受试者（=28）我,我= 1, …,n个.我们引入了一个潜在变量年我∗(t吨)，其中年我(t吨)=1，如果年我∗(t吨)≥0和年我(t吨)=0，如果年我∗(t吨)<0因此，年我(t吨)以伯努利分布，即。，年我(t吨)~伯尔尼(第页我(t吨))，其中第页我(t吨)=公共关系(年我∗(t吨)≥0)考虑到激励数据中的二进制响应是高度倾斜的，我们假设潜在变量的分布是倾斜的年我∗(t吨)如下：

年我∗(t吨)=μ我(t吨)+ϵ我(t吨)ϵ我(t吨)~GEV公司(0,ξ),
(2)

哪里μ我(t吨)确定潜在的位置参数年我∗(t吨)、和∈我(t吨)是GEV分布残差。以下内容Li等人（2016b），我们可以重写(2)作为GLMM，GEV链接为

第页我(t吨)=公共关系(年我(t吨)=1)=1−GEV公司(−μ我(t吨)；ξ)={1−经验{−(1−ξμ我(t吨))+−1ξ},ξ≠01−经验{−经验(μ我(t吨))},ξ=0}
(3)

其中GEV(x个；ξ)表示GEV分布的累积概率x个，带参数μ= 0,σ=1和形状ξ。的值μ和σ已修复可识别性问题。

在激发GAAD数据的背景下，牙齿水平的二元结果之间的空间相关性通过潜在向量表现出来μ我= (μ我(1), …,μ我(T型))'表示受试者我，因此μ我~MVN(x个我⊺β+ωα,∑)，其中MVN是具有平均值的多元正态密度x个我⊺β+ωα，的T型×T型正定方差-方差矩阵∑，使用x个我与参数向量相对应的主题级协变量矩阵β、和ω是带参数向量的齿级协变量矩阵α.多元潜在因素之间的空间关联μ我通过条件自回归（CAR）先验建模(贝萨格，1974年)的∑CAR模型假设协方差矩阵为μ我表单的σ2问(ρ)−1具有问(ρ) =D类–ρW，其中D类是一个T型×T型对角矩阵t吨第个对角线条目表示齿/位置的相邻数t吨、和W公司是邻接矩阵w个tt′=1，如果t吨和t吨'被视为相邻，否则为0。CAR模型包含两个参数；ρ∈[0,1]控制空间关联度，以及σ2（>0）控制变化幅度。在我们的数据分析中，与特定牙齿相邻的牙齿（位于同一个颌骨的两侧）被视为该牙齿的邻居。因此，邻域数∈{1,2}。

2.3. 非随机缺失数据

如前所述，在任何PD数据集中都可以观察到大量缺失的牙齿。这些缺失牙齿的位置不是随机的，并且与口腔该区域的牙周健康有关。为了在我们的模型中适应这种“信息性”，我们提出了一个“共享/联合随机效应”模型，将缺失牙齿的位置作为潜在过程的函数μ我(t吨). 共享参数模型为同时建模两个过程（丢失过程和（二进制）PD响应）提供了一个吸引人且直观的框架(Tsonaka等人，2009年). 让年我0(t吨)是（被观察到的）牙齿是否t吨主题缺少我，即。年我0(t吨)=1，如果是齿t吨缺少，否则为0。现在，引入一个潜在的连续变量年我0∗(t吨)，我们采用probit公式：

年我0(t吨)=我(年我0∗(t吨)>0),年我0∗(t吨)=一0+b0μ我(t吨)+ϵ我0(t吨),ϵ我0(t吨)~N个(0,1),
(4)

哪里一0和b0叙述潜在的过程μ我(t吨)牙齿缺失指示器。通过关节模型，μ我(t吨)包含在响应的存在性模型和二进制状态中，并且通过贝叶斯推理方案，存在性及其值都有助于μ我(t吨)因此β和α.何时b0=0，这对应于潜在μ我(t吨)以及牙齿缺失的位置，在这种情况下，牙齿的缺失不会影响估计β和α.

现在，在probit配方下年我0(t吨)=1，潜在变量年我0∗(t吨)遵循限制在零以上的截断正态分布；否则年我0∗(t吨)有界于零以下。我们可以将这个概率回归模型重写为公共关系(年我0(t吨) = 1) =Φ(一0+b0μ我(t吨))，其中Φ（·）是标准正态分布的累积密度函数。

三。 贝叶斯推断

在本节中，我们将介绍模型的贝叶斯推理框架。从模型参数的先验规范开始，我们提出了联合后验密度和相关条件，通过组合块Gibbs采样和内Gibbs HMC采样实现参数更新步骤。

3.1、。之前的选择

考虑由Θ= (β,α,σ2,ρ,ξ,一0,b0). 对于回归系数，我们设置β~MVN(一β,B类β)、和α~MVN(一α,B类α)，其中B类β和B类α都是对角矩阵。对于方差参数σ2在CAR规范中，我们可以选择非信息性的反式γIG(一σ,bσ)导致有条件的变戏法。空间关联参数的先验ρ是beta(一ρ,bρ)GEV形状参数ξ是正常的(一ξ,bξ). 统一的先验，其后验适当性在Wang和Dey（2010）也可以使用。最后，我们将非信息性平坦指定为Normal(μ0,σ02)对于系数一0和b0在probit链接函数中对非随机丢失进行建模。

3.2. 关节可能性和后验密度

结合联合（空间）GEV和缺失数据可能性（如第2节)在指定了先验值（没有特定的超参数值）的情况下，关节后部密度为：

π我=1n个MVN公司(μ我∣x个我⊺β+ωα；∑(σ2；ρ))×贝塔(ρ∣一ρ,bρ)×IG公司(σ2∣一σ,bσ)×MVN公司(α∣一α,B类α)×MVN公司(β∣一β,B类β)×正常(ξ∣一ξ,bξ)×π我=1n个πt吨=1T型[(1−GEV公司(−μ我(t吨)；ξ))年我(t吨)][(GEV公司(−μ我(t吨)))1−年我(t吨)]×π我=1n个πt吨=1T型q个我0(t吨)年我0(t吨)(1−q个我0(t吨))1−年我0(t吨)×N个(一0∣μ0,σ02)×N个(b0∣μ0,σ02)

哪里q个我0(t吨) =公共关系(年我0(t吨) = 1) =Φ(一0+b0μ我(t吨))，Beta（·，·）表示β密度，IG（·，.）表示反γ密度，MVN（·，…）表示多元正态密度，normal（·、…）表示正态密度。每个参数的完整条件分布如下：

  1. 回归系数β和α:
    β∣休息∝MVN公司(β∣一β,B类β)×π我=1n个MVN公司(μ我∣x个我⊺β+ωα；∑(σ2；ρ))~MVN公司(Gg（千兆克）,G公司),哪里G公司=(B类β−1+∑我=1n个x个我∑−1x个我⊺)−1,克=(一β⊺B类β−1+∑我=1n个μ我∑−1x个我⊺−α⊺ω⊺∑−1∑我=1n个x个我⊺)⊺
    α∣休息∝MVN公司(α∣一α,B类α)×π我=1n个MVN公司(μ我∣x个我⊺β+ωα,∑(σ2；ρ))~MVN公司(α∣小时,H（H）),哪里H（H）=(B类α−1+n个×ω⊺∑−1ω)−1,小时=(一α⊺B类α−1+∑我=1n个μ我⊺∑−1ω−∑我=1N个β⊺x个我∑−1ω)⊺.
  2. 空间差异σ2和关联参数ρ:
    σ2∣休息∝π我=1n个MVN公司(μ我∣x个我⊺β+ωα；∑(σ2；ρ))×IG公司(σ2∣一σ,bσ)~IG公司(σ2∣N个∗秒−第页2+一σ,bσ+∑我=1n个12[μ我−(x个我⊺β+ωα)]⊺(D类−ρW公司)[μ我−(x个我⊺β+ωα)])ρ∣休息∝贝塔(ρ∣一ρ,bρ)×π我=1n个MVN公司(μ我∣x个我⊺β+ωα,∑(ρ,σ2))
  3. 形状参数ξ:
    ξ∣休息∝正常(ξ∣一ξ,bξ)×π我=1n个πt吨=1T型[1−GEV公司(−μ我(t吨)；ξ)]年我(t吨)[GEV公司(−μ我(t吨)；ξ)]1−年我(t吨)
  4. 非随机缺失模型：回归系数（a）0,b0)和潜在变量年我0∗(t吨):
    [一0b0]∣休息∝π我=1n个πt吨=1T型N个(年我0(t吨)∣[1μ我(t吨)][一0b0],1)~MVN公司(嗯,M（M）)哪里M（M）=(∑我=1n个∑t吨=1T型[1,μ我(t吨)][1,μ我(t吨)]⊺)−1,米=∑我=1n个∑t吨=1T型年我0(t吨)[1,μ我(t吨)]⊺
    年我0∗(t吨)∣休息∝转弯.N个([1μ我(t吨)][一0b0],ϵ我0(t吨)){年我0∗(t吨)>0}我(年我0(t吨)=1)+转弯.N个([1,μ我(t吨)][一0b0],ϵ我0(t吨)){年我0∗(t吨)<0}我(年我0(t吨)=0)
  5. 共享的潜在因素μ我:
    μ我∣休息∝MVN公司(μ我∣x个我⊺β+ωα；∑(σ2；ρ))×πt吨=1T型[(1−GEV公司(−μ我(t吨)；ξ))年我(t吨)][(GEV公司(−μ我(t吨)))1−年我(t吨)]×πt吨=1T型q个我0(t吨)年我0(t吨)(1−q个我0(t吨))1−年我0(t吨)

3.3. MCMC工具

我们使用块Gibbs采样(Liu等人，1994年)算法，该算法从上述完整的条件分布中迭代地绘制参数/潜在变量块。参数的完整条件σ,β、和α为封闭形式；因此，我们可以根据其他参数直接绘制它们的后验样本。对于标量参数ρ和ξ（没有封闭形式的完整条件），我们可以在Gibbs采样中使用Metropolis-Hastings（M-H）算法来获得它们的后验值。计算挑战在于更新T型–潜在变量的维向量μ我对于我= 1, …,n个，它们是适度高维的、空间相关的，并且没有封闭形式的完整条件句。由于这些特性，使用within-Gibbs M-H算法更新μ我效率极低。相反，我们建议使用within-Gibbs HMC算法(Neal，2011年)实现MCMC链的更快收敛。

我们首先简要概述HMC采样，这是一种MCMC方法，它采用哈密顿动力学提出马尔可夫链中的未来状态，使马尔可夫链条能够更有效地探索高维复杂目标分布，从而加快收敛速度。哈密顿动力学是根据物体的位置描述物体运动的物理系统u个和动量v（v）时间τ，及其相关势能U型(u个)和动能K（K）(v（v）)分别为。系统的总能量是常数，称为哈密顿量H（H）(u个,v（v）)，它只是势能和动能的总和H（H）(u个,v（v）) =U型(u个) +K（K）(v（v）). 然后，可以使用微分方程组描述动能如何随时间转换为势能，反之亦然：

∂u个∂τ=∂H（H）∂v（v）=∂K（K）(v（v）)∂v（v）,∂v（v）∂τ=−∂H（H）∂u个=−∂U型(u个)∂u个.

利用这些哈密尔顿方程，我们可以预测物体在未来某个时间点的位置和动量τ′=τ+Δτ使用蛙跳法，该方法分割间隔Δτ进入之内L（左）较小的长度间隔δ并按顺序更新位置和动量变量L（左）间隔。具体来说，对于每个小间隔，蛙跳法分三步进行：

  1. v（v）(τ+δ∕2)=v（v）(τ)−(δ∕2)×∂U型(u个(τ))∂u个
  2. u个(τ+δ)=u个(τ)+δ×∂K（K）(v（v）(τ+δ∕2))∂v（v）
  3. v（v）(τ+δ)=v（v）(τ+δ∕2)−(δ∕2)×∂U型(u个(τ+δ))∂u个

蛙跳法可以L（左）近似预测的步骤u个和v（v）在τ′=τ+Δτ在MCMC采样中，哈密顿函数H（H）(u个,v（v）)使我们能够有效地探索特定的目标分布第页(u个)，并获得建议的状态u个这更有可能被接受。特别是，在我们的后验计算中，我们让位置变量为μ我动量是辅助变量ν然后，可以使用哈密顿动力学来提出μ我使用目标分布第页(μ我∣休息)，完全条件μ我MCMC取样中势能和动能函数的常见选择是

U型(μ我)=−日志(第页(μ我∣休息))=−日志[MVN公司(μ我∣x个我⊺β+ωα,∑(ρ,σ2))]−∑t吨日志[(1−GEV公司(−μ我(t吨)；ξ))年我(t吨)(GEV公司(−μ我(t吨)；ξ))1−年我(t吨)]−∑t吨日志[Φ(一0+b0μ我(t吨))年我0(t吨)(1−Φ(一0+b0μ我(t吨)))1−年我0(t吨)],K（K）(ν)=ν⊺ν2,

及其衍生物

∂U型(μ我)∂μ我=∑−1μ我−∑−1(x个我⊺β+ωα)−[年我(t吨)地理价值(−μ我(t吨)；ξ)1−GEV公司(−μ我(t吨)；ξ)+(1−年我(t吨))地理价值(−μ我(t吨)；ξ)GEV公司(−μ我(t吨)；ξ)]T型×1−[年我0(t吨)ϕ(一0+b0μ我(t吨))b0Φ(一0+b0μ我(t吨))+(1−年我0(t吨))ϕ(一0+b0μ我(t吨))b01−Φ(一0+b0μ我(t吨))]T型×1,∂K（K）(ν)∂ν=ν.

其中，gev（·）和gev（·）分别表示gev密度的pdf和cdf。更新时μ我，我们将上述导数插入三个蛙跳步骤中，运行该方法L（左）时间，并获得建议值μ我取决于其他参数。使用Metropolis验收标准，该建议值将作为马尔可夫链的下一步被确定为接受或拒绝。对于提议的μ我由于接受率高，链条运动快，因此选择合适的蛙跳步长至关重要δ，给定一个足够大的L（左）步长越大，接受率越低，而步长越小，链在目标密度勘探区的移动速度越慢，从而增加了计算时间(Neal，2011年). 运行多个测试样本后，我们选择了一个步长为δ=0.0005，并调整每个受试者的步长我在MCMC抽样的自适应阶段，结果接受率在75%到95%之间。

4 应用：GAAD数据

激励性GAAD数据（参见第1节)收集数据主要探讨PD与Hba1c或“糖化血红蛋白”测定的2型糖尿病水平之间的关系。27%的应答缺失，17%的未缺失应答属于第1类，即中度至重度PD。此外，还有几个受试者层面的协变量，即年龄（以岁为单位）、性别（1=女性，0=男性）、体重指数或BMI（以千克为单位）/米2)、吸烟状态（1=吸烟者，0=从不吸烟者）、血糖状态或Hba1c（1=高，0=受控）等。约26%的受试者是吸烟者。受试者的平均年龄约为52岁，范围为26-87岁。女性受试者在数据中占主导地位（约73%），这在该人群中并不罕见(Johnson-Spruill等人，2009年). 约74%的受试者肥胖（BMI>=30），64%的受测者Hba1c=1，这是高血糖水平的指标。此外，分析中考虑了一些牙齿水平的协变量，包括六个齿数指标T2-T7，其中T1、门牙作为参考，下颚指标（1=上颚或上颚，0=下颌或下颚）。在当前的分析中，我们从完整的数据库中随机选择了100名受试者的子集（具有完整的协变量信息）。

我们提出了四个相互竞争的模型来拟合数据，它们是

  • 模型1：具有非随机缺失的GEV链接，
  • 模型2：GEV链接缺失-随机，
  • 模型3：具有非随机缺失的probit链接，
  • 模型4：缺失随机的probit链接，

所有型号都具有中规定的CAR结构第2.2小节虽然模型1和模型2使用GEV链接解释了二进制响应中观察到的偏斜，但模型3和模型4使用对称概率链接。我们使用R进行了MCMC采样(R核心团队，2018年). 对于上的特定优先选择β和α，我们选择方法一α和一β作为0对于差异组件B类α)和B类β)，我们选择对角线元素（方差）为1000，表示弱信息信念(Wang和Dey，2011年)的β和α。对于σ2（CAR规范中的方差参数），我们最初选择一个逆Gamma(一σ,bσ)，其中一σ=bσ= 0.01. 然而，这导致了σ2以及其他参数。因此，我们使用了网格搜索方法σ2变化范围为0.1至1，增量为0.05（特别是对于模型1），并观察到以下情况获得的WAIC最低σ2= 0.25. 因此，我们修复了σ20.25，并进行MCMC计算的其余部分。空间关联参数的先验ρ紧随着严重左倾的贝塔系数（1,0.3），平均值集中在0.77，表示适度的空间关联。上一个ξ遵循Normal（0,1），而一0和b0是正常的(μ0= 0,σ02=10000). 对于所有4个模型，我们运行了2条初始值分散的链。对于每条链，10000的老化被认为是令人满意的，后验推断基于额外的10000个样本。使用跟踪图、自相关函数图和Gelman-Rubin确定收敛性R（右）^诊断。

我们使用WAIC比较了四种模型的性能(渡边，2010)准则–用于估计样本外期望的完全贝叶斯方法，以及偏差信息准则（DIC）。WAIC值基于elppd公司^W公司AIC公司=液化石油气^−第页^W公司AIC公司其中，第一项是估计的对数逐点后验预测（lppd）密度，第二项是lppd的估计方差，作为通过校正有效参数数可能出现的过拟合的惩罚项。这两项估计如下

液化石油气^=∑我=1n个∑t吨=1t吨日志(1N个∑k个=1N个第页(年我(t吨)∣μ我(t吨)k个,ξk个)),第页^W公司AIC公司=∑我=1n个∑t吨=1t吨V（V）应收账k个=1N个(日志第页(年我(t吨)∣μ我(t吨)k个,ξk个),

哪里k个是MCMC迭代编号。使其与AIC和DIC相当(Spiegelhalter等人，2002年)在偏差量表下，我们考虑WAIC公司=−2elppd公司^W公司AIC公司与AIC和DIC相比，WAIC在后验分布上取平均值，而不是以点估计为条件，因此从贝叶斯预测的角度来看，在评估实际用于新数据的预测时更为相关。

发件人表1，我们观察到模型1（具有非随机丢失的GEV链路模型）具有最小的WAIC和DIC值。我们还观察到，具有非随机缺失的模型比它们的缺失-随机对应模型好几倍。正如预期的那样，模型4（probit链接，缺少随机假设）提供了最差的拟合。这些结果表明了选择合适的链接函数对二进制数据建模的重要性。

表1：

估计4款竞争车型的WAIC和DIC。WAIC和DIC值越小，则表示模型拟合效果越好。

GEV链接Probit链接
缺少信息是的不是的不
WAIC公司3274.244053.7335564078.26
驾驶员信息中心3287.303832.273482.174076.31

通过绘制牙齿二进制响应的后验预测概率，我们进一步比较了模型拟合t吨对于随机受试者（受试者#52），以及相应的95%可信区间，对于牙齿数量（1-7），分别针对颌骨（上颌骨和下颌骨）。图2说明了具有不同链接功能和缺失假设的四个竞争模型之间的差异。总的来说，与上颌（上颌骨）相比，下颌骨（下颌骨）中的牙齿的预测估计值明显更好，可信区间更窄。我们还观察到，尽管有趣的是，缺失-随机模型的估计概率为更好的而不是非随机缺失模型。对下颌牙齿的拟合进行比较后发现，具有类似缺失假设的模型之间没有太大差异，尽管正如预期的那样，非随机缺失模型的拟合比缺失-随机模型的拟合好。

保存图片、插图等的外部文件。对象名为nihms-1534857-f0002.jpg

GAAD数据：（随机）受试者52的预测概率及其95%可信区间，在拟合4个模型后获得（如第4节). “圆圈”表示二进制响应的观测值，“实线”表示估计的后验概率，“虚线”是相应的95%逐点可信区间。上颚代表上颚的牙齿，而下颌代表下颚。

表2报告了参数的后验平均估计值及其95%可信区间，如果其可信区间不包括0，则认为参数具有显著性。从表中，我们观察到牙周健康随着年龄的增长而恶化。观察到男性PD水平高于女性的一些趋势，但效果并不显著。此外，BMI在非随机缺失模型中表现出较高PD的正向趋势，而在错误随机模型中表现为负向趋势，尽管这些影响也不显著。吸烟对较高的PD也有积极的影响，在缺少随机假设的情况下，其影响更大且显著。也有一些证据（不显著）表明，未经控制的HbA1c对PD有积极影响。仅从我们提出的模型（模型1）中可以看出，上颌牙齿的PD水平高于非上颌牙齿。此外，根据该模型，有明显证据表明，犬牙（T3）的牙周健康状况优于基线中切牙（T1），然而，PD的概率从T4（前磨牙）显著增加到T7（磨牙）。从研究对象中患病后牙（臼齿和前磨牙）占优势这一事实可以明显看出。对于T5-T7齿，非随机缺失模型的这一概率显著高于两个链路的相应缺失随机模型。后验估计b0对于非随机缺失模型，对应于潜在空间过程的斜率系数对于这两个连接函数都是正值且显著的，这证实了一种直觉，即牙周健康较差的受试者通常会有更多的缺失牙齿。GEV形状参数的估计ξis-ve（意味着有界的上尾GEV密度），在两种缺失类型下都显著，意味着疾病状态响应的二元分裂中存在偏斜。后验估计ρ所有4个模型都暗示了高度的空间关联。

表2：

在非随机缺失或随机缺失假设下，从GEV链接和概率链接模型导出的协变量系数的后验平均值和95%可信区间。

GEV链接Probit链接
非随机缺失是的不是的不
年龄0.6130.4360.7090.447
(0.356, 0.909)(0.009, 0.906)(0.423, 1.011)(0.025, 0.893)
女性−0.513−0.857−0.589−0.949
(−1.110, 0.055)(−1.884, 0.092)(−1.232, 0.067)(−1.828, 0.111)
体重指数0.118−0.2250.143−0.224
(−0.136, 0.385)(−0.715, 0.232)(−0.151, 0.436)(−0.692, 0.232)
吸烟者0.4781.4330.5871.444
(−0.066, 1.067)(0.502, 2.501)(−0.044, 1.216)(0.498, 2.406)
糖化血红蛋白0.2100.7820.2400.794
(−0.289, 0.723)(−0.087, 1.754)(−0.322, 0.807)(−0.081, 1.698)
马克西拉0.4760.3370.5410.331
(0.016, 0.968)(−0.459, 1.181)(−0.001, 1.084)(−0.474, 1.133)
T2段−0.146−0.195−0.181−0.199
(−0.389, 0.088)(−0.611, 0.206)(−0.494, 0.136)(−0.645, 0.241)
T3航站楼−0.395−0.424−0.506−0.457
(−0.704, −0.099)(−0.885, 0.038)(−0.821, −0.176)(−0.896, −0.018)
T4类0.316−0.2900.386−0.278
(0.058, 0.601)(−0.756, 0.174)(0.067, 0.697)(−0.755, 0.177)
T5类0.9440.6501.0790.668
(0.657, 1.316)(0.206, 1.202)(0.770, 1.399)(0.199, 1.137)
第6页2.1881.5012.5201.508
(1.786, 2.794)(0.928, 2.242)(2.160, 2.892)(1.004, 1.983)
第7天2.0301.6182.2991.646
(1.645, 2.576)(1.008, 2.363)(1.933, 2.676)(1.145, 2.159)
一0：缺少−0.215-−0.297-
(−0.345, −0.067)-(−0.379, −0.210)-
b0：缺少0.579-0.495-
(0.432, 0.724)-(0.416, 0.591)-
ρ 0.9960.9980.9970.998
(0.993, 0.998)(0.997, 0.999)(0.996, 0.998)(0.997, 0.999)
ξ −0.704−0.367--
(−1.118, −0.195)(−0.904, 0.021)--

使用GEV链接的一个主要障碍是其对回归系数的解释；它没有提供与逻辑回归框架中的比值比类似且有吸引力的含义。因此图3，我们提出了在各种协变量组合下观察中度到重度PD状态的后验预测概率。上面板显示了年龄为52.03岁、体重指数为37.41的随机肥胖受试者在不同性别、吸烟状况和糖化血红蛋白水平组合下上颌切牙的可能性。很明显，随机吸烟者患有非受控（高）HbA1c（FSH）的概率高于低HbA1c（FSL）的概率。高糖化血红蛋白组（MSH）男性吸烟者的相应概率也高于低糖化血红蛋白（MSL）组，这证实了先前研究的结果，这些研究探讨了吸烟状态与PD之间的关系(基纳和切斯特努特，2000年；Ah等人，1994年). 面板（b）显示了在上述各种组合下，年龄为70岁、体重指数为58的随机受试者的上颌尖牙的类似后向预测概率。这里，对于患有非受控糖化血红蛋白（FSH）的女性吸烟者来说，这种可能性再次高于低糖化血红蛋白的可能性。此外，对于吸烟和糖化血红蛋白状态相同的两组受试者，与随机的男性受试者相比，随机的女性受试者不太可能有中度到重度的PD状态。

保存图片、插图等的外部文件。对象名为nihms-1534857-f0003.jpg

后验预测概率Pr图(年我(t吨)新的=1）对于协变量的各种组合。F=女性，M=男性，S=吸烟，N=不吸烟，H=不受控（高）HbA1c，L=受控HbA1c。面板（a）显示了年龄=52.03且BMI=37.41的随机受试者的概率图，而面板（b）则显示了年龄=70且BMI=58的随机受检者的概率曲线。

5 模拟研究

在本节中，我们进行了一项小型模拟研究，以评估我们提出的模型及其比较器的有限样本性能。特别地，我们从相对偏差（RB）和均方误差（MSE）的角度研究了各种建模假设对参数推断的影响。我们认为T型=42（实际上相当于牙齿象限的站点级位置数量，包括7个牙齿），并且最初生成完整数据[年我(t吨)]n个×T型根据以下模型，使用主题级和牙齿级协变量，如下所示：

第页(年我(t吨)=1)=1−GEV公司(−μ我(t吨),ξ),μ我~MVN公司(x个我⊺β+ωα；∑(σ2；ρ)).

生成完整数据后，我们添加了如下缺失。让年我0(t吨)是缺失的指示器，因此第页(年我0(t吨) = 1) =Φ(一0+b0μ我(t吨)). 每个模拟数据集都有n个=50名受试者，主题级协变量x个我独立于logNormal（0,1）密度生成，具有相应的回归系数β= −1. 牙齿水平协变量ω由Bin（0.5）生成，具有相应的回归系数α= −2. 其他参数固定为ρ= 0.975,ξ= −0.5,一0=−1.5，以及b0= 0.5. 在这个设置下，我们生成了M（M）=100个数据集，具有倾斜的二进制结果和非随机缺失，每个数据集有50×42个二进制响应矩阵，平均缺失数据为25%，20%的观测数据为1。对于每个模拟数据集，我们拟合了相同的4个竞争模型（如第4节)依次，对参数进行相同的先验假设。

仿真结果如所示表3对于每个模型，我们总结了与RB和MSE的每个参数相对应的后验样本。与其他模型相比，模型1（具有非随机缺失的GEV链接）的回归系数的相对偏差和MSE最低β和α不出所料，忽略二进制响应中的偏度和信息缺失的模型（模型4）是四个竞争模型中估计效果最差的模型β和α与模型1相比，模型2（带有缺失随机假设的GEV链接）在估算方面表现出可比性β，但相对偏差和MSE较大α而模型3（具有非随机缺失的概率链）在估计中具有较大的RB和MSEβ和α有趣的是，模型2比模型3更好（就参数MSE和RB而言）。因此，选择合适的链接来适应二进制响应数据中的偏度特征似乎更为重要，而通过建模非随机缺失特征获得的信息则是一个优势。

表3：

模拟研究：列标签M1-M4表示拟合四个竞争模型的参数的相对偏差（RB）和均方误差（MSE），随链接函数和缺失假设的选择而变化。

参数M1级平方米立方米M4级
β 皇家银行0.003−0.0040.0250.030
MSE公司0.00290.00300.00330.0037
α 皇家银行0.0830.1170.0800.143
MSE公司0.03240.03530.07650.1010
σ 2 皇家银行0.0260.0390.0480.051
MSE公司0.00440.00430.00720.0092
ρ 皇家银行−0.007−0.007−0.002−0.002
MSE公司0.00010.00010.00020.0002
ξ 皇家银行−0.948−0.961--
MSE公司0.25800.2687--
一 0 皇家银行0.020-0.238-
MSE公司2.7280-3.3137-
b 0 皇家银行0.649-0.695-
MSE公司0.6911-0.7537-

虽然这四种模型在估计空间关联参数方面都有相似的性能ρ未指定的链接和/或缺失假设会影响空间方差的估计σ2特别是，模型2（缺少随机假设）的RB明显高于模型1，模型3和模型4的RB估计值越来越高。值得指出的是，在非随机缺失假设下，模型3（probit link）的估计RB和MSE要高得多一0和b0模型1（GEV链接）。这一结果与错误指定的链接功能可能导致模型性能更大下降的观察结果一致。

6 结论

在本文中，我们探讨了GEV链接函数对非对称二进制数据建模的潜力。在贝叶斯框架下，我们的建议同时考虑了激励数据中空间聚类和非随机缺失的其他有趣特征。我们的数据分析和仿真研究确定了选择适当的链接函数对有效估计模型参数的重要性。与对称probit或logit链接相比，GEV链接足够灵活，可以对正偏、对称或负偏二进制响应进行建模。

在我们的贝叶斯推断设置中，内吉布斯HMC算法有助于MCMC更新，比标准的随机行走Metropolis-Hastings算法更优雅。与使用HMC方法一次性更新所有参数相比，该混合方法利用了以下优点：（a）模型中大多数参数的封闭形式的完整条件，以及（b）更新潜在参数的并行性μ我在每个MCMC迭代中，从它们的完整条件开始。有条件的独立性μ我针对每个主题我给定其他参数有助于并行更新。因此，这种混合采样算法比使用全HMC一次性更新关节后方参数要高效得多。GitHub链接中提供了相关的R代码和分析数据集：https://github.com/bandyopd/GEV-Spatial网站

我们目前的模型比较仅限于模型构建和预测角度，而人们也可能对研究有影响力的观测结果的影响感兴趣(Goh和Dey，2014年)在后面。此外，为了贝叶斯计算的简单性，我们只将probit函数作为模型比较的基准对称链接。替代参数链接，如广义Pareto(霍斯金和沃利斯，1987年)，或结合非参数高斯过程先验和（参数）GEV链接函数的半参数命题(Li等人，2016b)也可以在建模具有非随机缺失的空间关联时考虑。此外，对偏斜顺序反应进行建模Wang和Dey（2011）也是可能的。所有这些都是未来可能的研究方向，将在其他地方进行。

致谢

作者感谢南卡罗来纳医科大学口腔健康研究中心提供了令人鼓舞的数据集，以及这项工作背后的背景。他们还感谢两位匿名审稿人，他们的建设性批评导致了手稿的改进。

Bandyopadhyay承认美国国立卫生研究院R01DE024984和P30CA016059（VCU Massey癌症中心支持拨款）的部分支持。

脚注

出版商免责声明：本作者接受的手稿是一份未经编辑的同行评审手稿的PDF文件，该手稿已被接受出版，但尚未被复制或更正。期刊上公布的官方记录版本是最新的，因此可能与本版本有所不同。

参与者信息

赵晓月，Amgen Inc.，加利福尼亚州千橡树，邮编：91320。

张林，明尼苏达州明尼阿波利斯市双子城明尼苏打大学生物统计学系，邮编：55455。

Dipankar Bandyopadhyay，弗吉尼亚州里士满弗吉尼亚联邦大学生物统计学系，邮编：23298。

工具书类

  • 美国儿科学会（2005）立场文件：牙周病流行病学.牙周病学杂志 76(8):1406–1419 [公共医学][谷歌学者]
  • Ah B、Michele K、Johnson GK、Kaldahl WB、Patil KD、Kalkwart KL（1994）吸烟对牙周治疗反应的影响.临床牙周病学杂志 21(2):91–97 [公共医学][谷歌学者]
  • Armitage GC（1999）牙周疾病和状况分类系统的开发.牙周病年鉴 4(1):1–6 [公共医学][谷歌学者]
  • Banerjee S、Carlin BP、Gelfand AE（2014）空间数据的层次化建模与分析，第2版。查普曼和霍尔/CRC[谷歌学者]
  • 贝萨格J（1974）空间相互作用与晶格系统的统计分析.英国皇家统计学会学报B辑（方法学）第192–236页[谷歌学者]
  • Carlin BP，Louis TA（2008）数据分析的贝叶斯方法，第3版。查普曼和霍尔/CRC[谷歌学者]
  • Chen MH、Dey DK、Shao QM（1999）二分量子响应数据的一种新的偏链模型.美国统计协会杂志 94(448):1172–1186[谷歌学者]
  • Coles S、Bawa J、Trenner L、Dorazio P（2001）极值统计建模简介，卷208.弹簧[谷歌学者]
  • Czado C，Santner TJ（1992）链接错误指定对二元回归推理的影响.统计规划与推理杂志 33(2):213–231[谷歌学者]
  • Fernandes JK、Wiegand RE、Salinas CF、Grossi SG、Sanders JJ、Lopes-Virella MF、Slate EH（2009）居住在南卡罗来纳州的非裔2型糖尿病患者的牙周疾病状况.牙周病学杂志 80(7):1062–1068[PMC免费文章][公共医学][谷歌学者]
  • Goh G，Dey DK（2014）基于函数Bregman散度的贝叶斯模型诊断.多元分析杂志 124:371–383[谷歌学者]
  • Hosking JR、Wallis JR（1987）广义Pareto分布的参数和分位数估计.技术计量学 29(三):339–349[谷歌学者]
  • Johnson-Spruill I、Hammond P、Davis B、McGee Z、Louden D（2009）南卡罗来纳州2型糖尿病患者Gullah家庭的健康：来自SuGar项目的糖尿病自我管理分析.糖尿病教育者 35(1):117–123[PMC免费文章][公共医学][谷歌学者]
  • Kim S、Chen MH、Dey DK（2008）柔性广义t吨-二进制响应数据的链接模型.生物特征 95(1):93–106[谷歌学者]
  • Kinane D，Chestnutt I（2000年）吸烟与牙周病.口腔生物学与医学评论 11(三):356–365 [公共医学][谷歌学者]
  • Li D、Wang X、Dey DK（2016a）基于广义极值分布和高斯过程先验的空间相关生存数据柔性治愈率模型.生物医学杂志 58(5):1178–1197 [公共医学][谷歌学者]
  • Li D、Wang X、Lin L、Dey DK（2016年b）高斯过程先验非参数二元回归中的柔性连接函数.生物计量学 72(三):707–719[PMC免费文章][公共医学][谷歌学者]
  • Liu JS、Wong WH、Kong A（1994）吉布斯采样器的协方差结构及其在估计量和增强方案比较中的应用.生物特征 81(1):27–40[谷歌学者]
  • McCulloch CE，Searle SR（2004）广义、线性和混合模型约翰·威利父子公司[谷歌学者]
  • McFadden D（1978）《住宅选址模型》。收件人：Karlqvist A、Lundqvist L、Snickers F、Weibull JW（编辑）空间相互作用理论与规划模型，荷兰北部，阿姆斯特丹，第75-96页[谷歌学者]
  • Neal RM（2011）MCMC使用哈密顿动力学。收件人：Brooks S、Gelman A、Jones GL、Meng XL（编辑）马尔可夫链蒙特卡罗手册，CRC Press/Taylor&Francis，佛罗里达州博卡拉顿，第5章，第113-162页[谷歌学者]
  • R核心团队（2018）R：统计计算的语言和环境。R统计计算基金会，奥地利维也纳，URLhttps://www.R-project.org/ [谷歌学者]
  • Reich BJ，Bandyopadhyay D（2010年）具有信息缺失的空间数据潜在因子模型.应用统计学年鉴 4(1):439–459[PMC免费文章][公共医学][谷歌学者]
  • Roy V，Dey DK（2014年）广义极值分布下分类模型和生存模型中后验分布的适用性.中国统计局第699–722页[谷歌学者]
  • Spiegelhalter DJ、Best NG、Carlin BP、Van Der Linde A（2002）模型复杂性和拟合的贝叶斯度量.英国皇家统计学会杂志：B辑（统计方法） 64(4):583–639[谷歌学者]
  • Tsonaka R、Verbeke G、Lesaffre E（2009）处理非单调不可忽略丢失的半参数共享参数模型.生物计量学 65(1):81–87 [公共医学][谷歌学者]
  • Wang X，Dey DK（2010年）二进制响应数据的广义极值回归：在B2B电子支付系统采用中的应用.应用统计学年鉴 4(4):2000–2023[谷歌学者]
  • Wang X，Dey DK（2011年）有序响应数据的广义极值回归.环境与生态统计 18(4):619–634[谷歌学者]
  • 渡边S（2010）奇异学习理论中贝叶斯交叉验证的渐近等价性和广泛适用的信息准则.机器学习研究杂志 11（12月）：3571–3594[谷歌学者]
  • Wiebe CB，Putnins EE（2000）美国牙周病学会牙周病分类系统——更新.加拿大牙科协会杂志 66(11):594–599 [公共医学][谷歌学者]