Simultaneous spatial smoothing and outlier detection using penalized regression, with application to childhood obesity surveillance from electronic health records

Young-Geun Choi; Lawrence P. Hanrahan; Derek Norton; Ying-Qi Zhao

doi:10.1111/biom.13404

生物识别。作者手稿；PMC 2022 5月27日提供。

以最终编辑形式发布为：

生物计量学。2022年3月；78(1): 324–336.

2020年12月11日在线发布。数字对象标识：10.1111/生物13404

预防性维修识别码：PMC9138186

美国国立卫生研究院：美国国家卫生研究院1808902

PMID：33215685

基于惩罚回归的同步空间平滑和离群值检测，应用于电子健康记录中的儿童肥胖监测

杨根财,¹ 劳伦斯·P·汉拉罕,² 德里克·诺顿,^三和赵英琪^4,^*

作者信息版权和许可信息 PMC免责声明

关联数据

补充资料: 补充材料。
NIHMS1808902-补充-补充材料.pdf（417000）
GUID:210956FF-7FDD-438C-AC56-F185C72468AF

数据可用性声明: 由于道德限制，支持本文研究结果的数据尚未公开。感兴趣的读者可以联系威斯康星大学IRB委员会(https://www.medicine.wisc.edu/research/uw-health-sciences-irbs)了解详细信息。

总结：

近年来，电子健康记录（EHR）已成为数据驱动的细粒度监测平台。在本文中，我们利用EHR早期预防儿童肥胖。该方法同时为肥胖流行提供了平滑的疾病图和异常值信息，有助于提高公众意识和促进有针对性的干预。更准确地说，我们考虑一个惩罚多级广义线性模型。我们将区域贡献分解为平滑和稀疏信号，这些信号通过融合和对似然函数施加稀疏惩罚的组合自动识别。此外，我们对提议的可能性进行了权衡，以解释EHR数据的缺失和潜在的不代表性。我们开发了一种新的交替最小化算法，该算法计算效率高，易于实现，并保证收敛。仿真研究证明了该方法的优越性能。最后，我们将我们的方法应用于威斯康星大学人口健康信息交换数据库。

关键词：儿童肥胖监测、疾病测绘、电子健康记录、融合惩罚、离群值检测、稀疏惩罚

1.简介

预防儿童肥胖对于控制全球肥胖流行变得越来越重要。在资金有限的情况下，需要对儿童肥胖进行细致的监测，以确定和跟踪肥胖趋势，帮助设计干预措施并指导政策解决方案(Longjohn等人，2010年). 日常收集的大量健康数据库，如电子健康记录（EHR），作为评估趋势和当地儿童肥胖风险的平台，正受到关注(Friedman等人，2013年).

地理空间监测的统计方法可能包括两个方面：一是监测区域流行趋势（也称为“疾病绘图”），二是确定不同地点流行率的意外变化（也称为”热点检测“）。传统上，这两项任务是分别完成的。对于任务i），肥胖文献主要使用标准广义线性混合效应模型（GLMM）来解释个人因素和社区环境。这些方法假设区域随机效应是独立的，尽管即使调整协变量后也存在空间依赖性(Panczak等人，2016年). 为了解释空间相关性，从频率学家和贝叶斯的角度提出了平滑疾病映射的方法。在泊松对数线性模型或多级logistic模型下，区域特定效应通过核函数平滑(Ghosh等人，1999年)或花键(Ugarte等人，2010年)或被条件自回归（CAR）先验建模为相关随机向量(Besag等人，1991年;Mercer等人，2015年). 这些策略产生了“集群”风险图，增强了可解释性，但没有探索异常区域的识别。对于任务ii），最流行的方法是空间扫描统计方法(库尔多夫和纳加瓦拉，1995年;Jung，2009年). 扫描统计方法搜索预先指定的一组地理区域，并进行广义似然比检验，以测试事件的比例在区域之间、区域内部和区域外部是否均匀。然而，它可能不适合识别大小不同的多个位置。回归方法产生的残差也可用于检测区域疫情，具有较大残差的观测值被视为离群值(Farrington等人，1996年;Zhao等人，2011年). 然而，当一个异常值是一个杠杆点或存在多个异常值时，已知基于残差截止值的异常值检测会失败(她和欧文，2011年).

融合惩罚用于平滑是在最小二乘法中首次提出的(Tibshirani等人，2005年)，然后用于公共卫生研究(Wang和Rodríguez，2014年). 融合惩罚的拟合结果似乎是分段常量，从而产生拟合值的自然聚类。通过融合惩罚进行平滑，可以使用不同的惩罚进行额外的正则化，例如稀疏惩罚，这在其他平滑疾病映射方法中可能并不简单。异常值检测的稀疏惩罚与平方误差损失一起使用(Kim等人，2009年;Tibshirani和Taylor，2011年;她和欧文，2011年).Kim等人（2009年）和Tibshirani和Taylor（2011）考虑了ℓ₁罚款，以及她和欧文（2011）报告称，非凸处罚优于这两种处罚ℓ₁标准多元线性回归中基于惩罚和截断的检测方法。

我们开发了一种新的方法，可以同时生成可解释的疾病图并检测异常区域。我们建立了一个多级逻辑模型来自然地纳入风险因素。一种新的混合正则化方法包括表示区域特定效应的平滑信号和稀疏信号。平滑信号通过融合惩罚进行正则化，以便相邻位置具有相似的拟合基线肥胖率。对稀疏信号实施非凸稀疏惩罚，以便非零拟合系数表示潜在的离群值。值得一提的是，由于遗漏和非代表性，从EHR估算人口健康指标可能具有挑战性。以下Flood等人（2015），我们采用两步加权程序来解释缺失数据并调整全国代表性样本的协变量分布。

我们最初的贡献是双重的。首先，当融合和ℓ₁线性模型中考虑了惩罚(Kim等人，2009年;Tibshirani和Taylor，2011年)，据我们所知，我们是第一个采用融合惩罚和非凸惩罚来识别异常值的。其次，我们提供了一种有效的优化算法，可以保证混合正则化模型的收敛性，并可以利用现有的软件包。虽然我们的算法是用贝努利似然描述的，但它可以很容易地扩展到处理其他凸损失函数。

在第2节，我们描述了激发我们研究的威斯康星大学电子健康记录公共健康信息交换（PHINEX）数据库，并在第3节。仿真研究见第4节，这证明了我们提出的方法的优越性能。我们将我们的方法应用于PHINEX儿童肥胖监测第5节。我们在第6节.

2.数据

威斯康星大学电子健康记录公共卫生信息交换（UW eHealth PHINEX）数据库包含来自威斯康星州中南部学术医疗系统的EHR数据。它包括2007年至2012年期间在家庭医学、儿科和内科诊所进行初级护理的患者记录。所有PHINEX数据均来自Epic EHR Clarity Database（EpicCare电子病历，Epic Systems Corp.，维罗纳州）。此外，该计划将地理代码与人口普查区块组联系起来，并将EHR与社区层面的健康社会决定因素联系起来。它的创建是为了通过了解疾病风险、患者和社区的局部变化来改善临床实践和人群健康(吉尔伯特等人，2012年).

在这篇论文中，我们关注了2011-2012年间93130名年龄在2-19岁的患者。体重指数（BMI）值（英寸千克/米²)根据受试者在同一次就诊时测量的身高和体重计算得出。任何BMI达到或超过第95百分位的受试者都被归类为肥胖。在所有患者中，34852人（37.4%）缺少有效的BMI。个人层面的协变量包括性别、年龄、种族/民族、医疗服务付款人（即保险）以及2010年人口普查分组中受试者居住地的信息。特定地区的协变量包括经济困难指数（EHI）和区块组的城市化程度，其中EHI(内森和亚当斯，1989年)用于衡量区块组的社会经济地位，并对所有威斯康星州人口普查区块组进行标准化。据统计，一个人口普查分组的城市化基于其11个城市化汇总分组ESRI（2012）这些群体来源于人口普查分组人口密度、城市规模、与大都市地区的接近程度以及经济/社会中心性的数据。城市度整数值范围为1（城市最多）到11（农村最多）。

3.方法

3.1、。模型设置

我们使用双下标，ij公司(j个= 1, …,n个_我,我=1…，K（K）)以指示j个-第个主题我-第个区域。让S公司_我是的质心位置我-第个区域。让X（X）_我表示区域层面的协变量，如城市化和EHI。让Y（Y）_ij公司是肥胖的指标(ij公司)-第个主题，带有Y（Y）_ij公司=1表示肥胖。最后，让我们Z轴_ij公司是协变量的向量(ij公司)-性别、年龄、种族/民族和保险支付人等主题。

让 $对_{我 j个} = ℙ ({Y（Y）}_{我 j个} = 1 ∣ {Z轴}_{我 j个}, {X（X）}_{我})$ 。我们将模型形式化对_ij公司作为

罗吉特 (对_{我 j个}) = {Z轴}_{我 j个}^{T型} α_{1} + {X（X）}_{我}^{T型} α_{2} + β_{我} + γ_{我},

(1)

服从 \sum_{我_{1} < 我_{2}} ρ_{我_{1}, 我_{2}} | β_{我_{1}} - β_{我_{2}} | ⩽ {c（c）}_{1};

(2)

\sum_{我 = 1}^{K（K）} 我 (γ_{我} \neq 0) ⩽ {c（c）}_{2},

(3)

哪里c（c）₁,c（c）₂⩾0，逻辑(t吨)=对数{t吨/(1−t吨)}、和我（·）是指示器功能。这个β_我s代表了肥胖流行的地区贡献，而不是由个人或其他地区水平的特征来解释。由于儿童肥胖的可能性可能会受到社区环境的影响，我们预计邻近地区的个人对肥胖流行率的区域贡献相似(Panczak等人，2016年)，因此是平滑度约束(2)被强加给β_我.融合重量 $ρ_{我_{1}, 我_{2}} (ρ_{我_{1}, 我_{2}} ⩾ 0)$ 表示每对我₁和我₂。更高的值 $ρ_{我_{1}, 我_{2}}$ 将导致更相似的一对 $β_{我_{1}}$ 和 $β_{我_{2}}$ 。通过适当的调整参数选择β_我如果相似的位置被分组在一起，则可能占用的空间有限。我们可以解释β_我作为区域的分割或聚类。γ_我用于捕获潜在的异常区域，其中我-th区是一个异常值，如果γ_我≠ 0. 给定稀疏约束(三)，我们期待γ_我对于大多数区域，将为零（非异常值），但少数区域可能为非零（异常值）。我们的公式可以被视为王和罗德里格斯（2014），其中除了融合约束外，我们还在每个区域上添加了稀疏约束。这一增加使得在区域一级调整依赖性后能够捕获异常疫情。

如果满足以下条件，则模型是可识别的γ_我的距离足够远，但可能并非如此。虽然分离信号的想法在她和欧文（2011）,Kim等人（2009）,Tibshirani（2014），以及Chernozhukov等人（2017），确定模型可识别性的边界值尚未正式研究，需要进一步研究。然而，请注意 ${α_{我 j个}, β_{我} + γ_{我}}_{1 ⩽ j个 ⩽ {n个}_{我}, 1 ⩽ 我 ⩽ K（K）}$ 可识别。

3.2. 使用完整数据进行估算

表示 $N个 = \sum_{我 = 1}^{K（K）} {n个}_{我}$ , $α = {(α_{1}^{T型}, α_{2}^{T型})}^{T型}$ ,β= (β₁, …,β_K（K）)^T型和γ= (γ₁, …,γ_K（K）)^T型如果所有患者都有完整的记录，则可以通过惩罚逻辑可能性来估计参数，其中 $(\hat{α}, \hat{β}, \hat{γ}) = {argmin（最小值）}_{α, β, γ} ϕ (α, β, γ)$ ，以及目标函数ϕ是

ϕ (α, β, γ) = - 洛格里克 (α, β, γ) + {P（P）}_{λ_{1}} (β) + 问_{λ_{2}} (γ) .

(4)

归一化负对数似然函数为

- 洛格里克 (α, β, γ) = \frac{1}{N个} \sum_{我 = 1}^{K（K）} \sum_{j个 = 1}^{{n个}_{我}} [日志 {1 + 经验 ({Z轴}_{我 j个}^{T型} α_{1} + {X（X）}_{我}^{T型} α_{2} + β_{我} + γ_{我})} - {Y（Y）}_{我 j个} ({Z轴}_{我 j个}^{T型} α_{1} + {X（X）}_{我}^{T型} α_{2} + β_{我} + γ_{我})] .

(5)

第二任期 ${P（P）}_{λ_{1}} (β)$ 是源自拉格朗日量的聚变惩罚(2)，其中 ${P（P）}_{λ_{1}} (β) = λ_{1} \sum_{我_{1} < 我_{2}} ρ_{我_{1}, 我_{2}} | β_{我_{1}} - β_{我_{2}} |$ 。我们使用 $ρ_{我_{1}, 我_{2}} = 1 / d日 ({S公司}_{我_{1}}, {S公司}_{我_{2}})$ ，其中 $d日 ({S公司}_{我_{1}}, {S公司}_{我_{2}})$ 表示之间的距离 ${S公司}_{我_{1}}$ 和 ${S公司}_{我_{2}}$ 。此处，地理距离用于定义d日（·，·），但可以使用其他相似性度量。在不失一般性的情况下，我们假设 ${最大值}_{我_{1}, 我_{2}} ρ_{我_{1}, 我_{2}} = 1$ ，否则我们可以对其进行规范化。由于涉及融合惩罚的优化的计算成本在非零的数量上呈二次增加 $ρ_{我_{1}, 我_{2}} ’ 秒$ ，可能需要保留一些 $ρ_{我_{1}, 我_{2}} 秒$ 并将其他值截断为零，以便于计算。

第三学期， $问_{λ_{2}} (γ) = \sum_{我 = 1}^{K（K）} {n个}_{我} {q个}_{λ_{2}} (γ_{我}) / N个$ ，是一个稀疏的点球，是拉格朗日函数的松弛(三)，其中q个_λ（·）是一个单变量惩罚函数。特别是，我们考虑了她和欧文（2011）,q个_λ(t吨) = (λ|t吨|负极t吨²/2)我(t吨<λ) +λ²/2我(t吨⩾λ). 硬性惩罚导致非凸公式(4)，保证收敛到局部极小值。我们称一下我-第次罚款 $问_{λ_{2}} (γ)$ 通过n个_我这样，不同地区的受试者受到同等惩罚。

3.3. 优化算法

我们开发了一种交替最小化算法。它交替更新α,β，以及γ每次将其中一个最小化，同时保持其他的不变。表示当前迭代α^(t吨),β^(t吨)，以及γ^(t吨)此外，我们表示 $问_{我 j个} = ({Z轴}_{我 j个}^{T型}, {X（X）}_{我}^{T型})$ .然后 $问_{我 j个}^{T型} α = {Z轴}_{我 j个}^{T型} α_{1} + {X（X）}_{我}^{T型} α_{2}$ .

正在更新α.修复β=β^(t吨)和γ=γ^(t吨).目标函数等价于

ϕ (α, β^{(t吨)}, γ^{(t吨)}) = \frac{1}{N个} \sum_{我 = 1}^{K（K）} \sum_{j个 = 1}^{{n个}_{我}} [日志 {1 + 经验 (问_{我 j个}^{T型} α + μ_{我 j个}^{(t吨)})} - {Y（Y）}_{我 j个} (问_{我 j个}^{T型} α + μ_{我 j个}^{(t吨)})]

具有 $μ_{我 j个}^{(t吨)} = β_{我}^{(t吨)} + γ_{我}^{(t吨)}$ ，对应于上的经典逻辑回归N个个人。可以运行标准包（例如glm公司在里面R（右）)以获得α^(t吨+1).

正在更新β.修复α=α^(t吨+1)和γ=γ^(t吨)，然后

ϕ (α^{(t吨 + 1)}, β, γ^{(t吨)}) = \underset{= : 我 (β)}{\underset{︸}{\frac{1}{N个} \sum_{我 = 1}^{K（K）} \sum_{j个 = 1}^{{n个}_{我}} [日志 {1 + 经验 (β_{我} + θ_{我 j个}^{(t吨)})} - {Y（Y）}_{我 j个} (β_{我} + θ_{我 j个}^{(t吨)})]}} + λ_{1} \sum_{我_{1} < 我_{2}} ρ_{我_{1}, 我_{2}} | β_{我_{1}} - β_{我_{2}} |,

哪里 $θ_{我 j个}^{(t吨)} = 问_{我 j个}^{T型} α^{(t吨 + 1)} + γ_{我}^{(t吨)}$ 对于每个我和j个。为了简单起见，请定义ψ(β) =ϕ(α^(t吨+1),β,γ^(t吨))，在中是凸的β。要更新β^(t吨)，我们建议最小化代理目标函数，其中我(β)被其周围的局部二次近似所取代β^(t吨).

写出的二阶泰勒展开式我(β)在β^(t吨)作为

\tilde{我} (β; β^{(t吨)}) = 我 (β^{(t吨)}) + \nabla_{β} 我 {(β^{(t吨)})}^{T型} (β - β^{(t吨)}) + \frac{1}{2} {(β - β^{(t吨)})}^{T型} \nabla_{β β}^{2} 我 (β^{(t吨)}) (β - β^{(t吨)}),

其中+_β和 $\nabla_{β β}^{2}$ 是关于β.将代理目标函数定义为 $\tilde{ψ} (β; β^{(t吨)}) = \tilde{我} (β; β^{(t吨)}) + {P（P）}_{λ_{1}} (β)$ 。我们计算 $\tilde{β} = 参数 {最小值}_{β} \tilde{ψ} (β; β^{(t吨)})$ ，其中

\tilde{β} = \underset{β}{argmin（最小值）} [\frac{1}{2} \sum_{我 = 1}^{K（K）} {A类}_{我}^{(t吨)} {(β_{我} - {B类}_{我}^{(t吨)})}^{2} + λ_{1} \sum_{我_{1} < 我_{2}} ρ_{我_{1}, 我_{2}} | β_{我_{1}} - β_{我_{2}} |],

具有

{A类}_{我}^{(t吨)} = \sum_{j个 = 1}^{{n个}_{我}} \frac{经验 (β_{我}^{(t吨)} + θ_{我 j个}^{(t吨)})}{{1 + 经验 (β_{我}^{(t吨)} + θ_{我 j个}^{(t吨)})}^{2}}; {B类}_{我}^{(t吨)} = β_{我}^{(t吨)} - \frac{1}{{A类}_{我}^{(t吨)}} \sum_{j个 = 1}^{{n个}_{我}} [\frac{经验 (β_{我}^{(t吨)} + θ_{我 j个}^{(t吨)})}{1 + 经验 (β_{我}^{(t吨)} + θ_{我 j个}^{(t吨)})} - {Y（Y）}_{我 j个}] .

用于计算 $\tilde{β}$ ，我们应用了由Yu等人（2015），它产生了一个稳定的解决方案，并且易于实现。

确保 $ψ (β^{(t吨)}) ⩾ ψ (\tilde{β})$ ，我们采用Lee等人（2016）的一步修改 $\tilde{β}$ ：如果 $ψ (β^{(t吨)}) ⩾ ψ (\tilde{β})$ ，让 $β^{(t吨 + 1)} = \tilde{β}$ ; 否则， $β^{(t吨 + 1)} = \tilde{小时} \tilde{β} + (1 - \tilde{小时}) β^{(t吨)}$ ，其中 $\tilde{小时} = {argmin（最小值）}_{小时 \in [0, 1]} ψ (小时 \tilde{β} + (1 - 小时) β^{(t吨)})$ 。我们将在中显示提议1那个 $\tilde{小时}$ 始终存在并且ψ(β^(t吨)) ⩾ψ(β^(t吨+1))暂停迭代。

正在更新γ.鉴于此α=α^(t吨+1)和β=β^(t吨+1),

ϕ (α^{(t吨 + 1)}, β^{(t吨 + 1)}, γ) = \frac{1}{N个} \sum_{我 = 1}^{K（K）} \sum_{j个 = 1}^{{n个}_{我}} [日志 {1 + 经验 (γ_{我} + ν_{我 j个}^{(t吨)})} - {Y（Y）}_{我 j个} (γ_{我} + ν_{我 j个}^{(t吨)})] + \frac{1}{N个} \sum_{我 = 1}^{K（K）} {n个}_{我} {q个}_{λ_{2}} (γ_{我}),

哪里 $ν_{我 j个}^{(t吨)} = 问_{我 j个}^{T型} α^{(t吨 + 1)} + β_{我}^{(t吨 + 1)}$ .稍微滥用符号，我们定义了一个单变量目标函数ϕ_我(γ)和损失函数我_我(γ) (我= 1, …,K（K）)作为

ϕ_{我} (γ) = \underset{我_{我} (γ)}{\underset{︸}{\sum_{j个 = 1}^{{n个}_{我}} [日志 {1 + 经验 (γ + ν_{我 j个}^{(t吨)})} - {Y（Y）}_{我 j个} (γ + ν_{我 j个}^{(t吨)})]}} + {n个}_{我} {q个}_{λ_{2}} (γ) .

很明显 $ϕ (α^{(t吨 + 1)}, β^{(t吨 + 1)}, γ) = {N个}^{- 1} \sum_{我 = 1}^{K（K）} ϕ_{我} (γ_{我})$ 因此，只需优化K（K）单变量函数ϕ_我(·),我= 1, …,K（K）.尽管每个ϕ_我(γ)是非凸的，我们可以找到ϕ_我如下所示。让 $\tilde{t吨} = {argmin（最小值）}_{t吨 \in ℝ} 我_{我} (t吨)$ .自 ${q个}_{λ_{2}} (\cdot)$ 在[−外为常数λ₂,λ₂]，最小值为ϕ_我（·）位于[−λ₂,λ₂]或等于 $\tilde{t吨}$ 因此，我们提出了一种网格搜索方法。让{t吨₁, …,t吨_T型} ⊆ [−λ₂,λ₂]、和 ${\hat{γ}}_{我}^{(t吨 + 1)} = {argmin（最小值）}_{γ \in {\tilde{t吨}, {t吨}_{1}, \dots, {t吨}_{T型}}} ϕ_{我} (t吨)$ .

完整的算法见Web附录C。所提出的算法保证了以下属性。

P（P）提议1：假设每个我，存在j个₁,j个₂这样的话 ${Y（Y）}_{我 {j个}_{1}} = 0$ 和 ${Y（Y）}_{我 {j个}_{2}} = 1$ .任何选择α^(t吨),β^(t吨)，以及γ^(t吨)，更新的迭代α^(t吨+1),β^(t吨+1)，以及γ^(t吨+1)中的算法1Web附录C满足单调递减性质：ϕ(α^(t吨),β^(t吨),γ^(t吨)) ⩾ϕ(α^(t吨+1),β^(t吨),γ^(t吨)) ⩾ϕ(α^(t吨+1),β^(t吨+1),γ^(t吨)) ⩾ϕ(α^(t吨+1）,β^(t吨+1),γ^(t吨+1)).

证据推迟到Web附录A假设表明，幼稚流行率 $\sum_{j个 = 1}^{{n个}_{我}} {Y（Y）}_{我 j个} / {n个}_{我}$ 每个都依赖于（0，1）我，这对于保证每个步骤中存在最优值至关重要。由提议1，的任何极限点{(α^(t吨),β^(t吨),γ^(t吨))}是一个静止点，如果ϕ是连续的。由于目标函数ϕ由于该算法是非凸的，因此只能保证收敛到局部最优解，并且需要仔细选择初始点。我们可以使用热启动策略，其中上一个优化参数下的解决方案用作下一个优化选项的初始点。该策略在我们的数值研究中实施时表现良好。此外，将所述算法扩展到其他（多级）广义线性模型也很简单。我们仍然可以解决α-使用非现成软件包（例如。glm公司在里面R（右）)、和β-和γ-使用相同策略的步骤。

3.4. 调谐参数的选择

我们实现了一个模型选择过程来调整λ₁和λ₂。我们使用了在她和欧文（2011）, ${银行识别码}^{*} (λ_{1}, λ_{2}) = - 2 N个 \cdot 洛格里克 (\hat{α}, \hat{β}, \hat{γ}) + DF公司 \cdot (1 + 日志 N个)$ .给， $洛格里克 (\hat{α}, \hat{β}, \hat{γ})$ 定义于(5)，自由度（DF）是通过组合套索回归和融合套索回归中计算的DF来计算的(Tibshirani等人，2005年;邹等人，2007;Tibshirani和Taylor，2011年)，其中

DF公司 = (的尺寸 \hat{α}) + (# 不同值的 \hat{β}) + (# 的非零值 \hat{γ}) .

(6)

我们搜索了(λ₁,λ₂)在最小化BIC的候选集中*(λ₁,λ₂).

3.5. 考虑缺失和选择偏差的权重

如前几节所述，我们的数据集涉及大量肥胖指标的缺失值(Y（Y）_ij公司). 此外，数据可能无法直接与国家样本进行比较。我们考虑两步加权程序来调整缺失的BMI值和选择偏差。

第一步是解释BMI的缺失。我们假设BMI的随机缺失（MAR），其中缺失BMI的概率独立于其对协变量的响应条件(Little和Rubin，2014年). 让R（右）_ij公司=1，如果Y（Y）_ij公司被观察到并且R（右）_ij公司否则=0。体重被定义为观察BMI的逆概率， $ℙ ({R（右）}_{我 j个} = 1 ∣ {Z轴}_{我 j个}, {X（X）}_{我})$ ，可以通过logistic回归进行估计。第二步是根据年龄、性别和种族/民族的人口分布进行调整。我们使用2012年全国人口普查数据进行了分层后校正。每个受试者的最终权重是逆概率权重和分层后权重的乘积。相应地修改了目标函数和后续程序。此外，我们使用了一个具有一阶正态近似的bootstrap方法（参见示例。Puth等人，2015年和Efron和Tibshirani，1994年)构造置信区间，以解释缺失和选择偏差的不确定性，其中对每个重采样数据集重新计算权重和模型估计。有关详细信息，请参阅Web附录B.

4.模拟研究

我们将该方法与经典的GLMM、具有条件自回归随机效应的GLMM（GLMM-CAR）以及Jung（2009）（扫描统计信息）。GLMM假设 $罗吉特 (对_{我 j个}) = {Z轴}_{我 j个}^{T型} α_{1} + {X（X）}_{我} α_{2} + {第页}_{我} + δ$ 其中(第页₁, …,第页_K（K）)^T型~MVN(0,我_K（K）)是独立随机效应，以及δ是全局截距。GLMM-CAR的模型为 $罗吉特 (对_{我 j个}) = {Z轴}_{我 j个}^{T型} α_{1} + {X（X）}_{我} α_{2} + {b条}_{我} + {第页}_{我} + δ$ ，其中(b条₁, …,b条_K（K）)^T型~MVN(0,Σ)表示空间平滑随机效果，以及Σ享受由提出的CAR模型中的形式Besag等人（1991年）。我们使用了该函数S.CARMillivel（）属于R（右）包裹CARBayes公司使用默认选项实现GLMM-CAR。一旦通过GLMM和GLMM-CAR拟合模型，我们就使用基于截断的方法来识别异常区域。让 ${\hat{第页}}_{我}$ 是预测的随机效应我-第个区域。这个我-第th个区域被声明为异常值，如果 $| {\hat{第页}}_{我} | > 2.5 \hat{σ}$ ，其中 $\hat{σ}$ 是估计的标准偏差第页_我.切断 $2.5 \hat{σ}$ 是文学中的一个流行选择(她和欧文，2011年). 协变量调整扫描统计方法(Jung，2009年)假设 $罗吉特 (对_{我 j个}) = {Z轴}_{我 j个}^{T型} α_{1} + {X（X）}_{我} α_{2} + 我 (我 \in S公司) θ + δ$ .给，S公司表示一组区域。对于每个 $S公司 \in S公司$ 该方法反复拟合模型并计算似然比检验统计量进行检验H（H）₀:θ= 0. 然后该方法选择 ${S公司}_{0} \in S公司$ 如果相应的轻轨统计数据最大，则作为热点。我们还包括三个“oracle”版本的方法，其中ϕ相对于以下之一最小化α,β，或γ而其他两个设置为真值：相对于α（甲骨文α)；关于β（甲骨文β)；关于γ（甲骨文γ).

我们考虑过K（K）(K（K）=20，40）个地区，每个地区的受试者人数n个(n个= 50, 100). 我们生成了问_ij公司= (Z轴_ij公司,X（X）_我)^T型，其中Z轴_ij公司和X（X）_我从概率为0.5的伯努利分布中得出。我们设置了α= (α₁,α₂) = (−0.2, 0.2). 为了简单起见，我们模拟了K（K）一维直线上的位置S公司_我~统一（5，95），我= 1, …,K（K）.β_我如果为5⩽，则设置为logit（0.4）S公司_我<35，如果35⩽，则logit（0.5）S公司_我<65，如果65⩽，则为logit（0.6）S公司_我⩽ 95. 我们随机选择K（K）_{O（运行）}地区，其中γ_我对于⌊=2K（K）₀/2个地区（t吨⌋是不大于的最大整数t吨)和γ_我=−2表示其余。因此K（K）_{O（运行）}区域，具有γ_我≠0是离群值。我们改变了异常值的数量，以便K（K）_{O（运行）}/K（K）= 0%, 5%, 10%, 15%. 对于每个场景，我们重复生成1000个数据集。我们对每个数据集应用了不同的方法，并评估了性能指标。性能度量平均超过1000次复制。使用建议的修改BIC，在预定义的候选集中选择调谐参数λ₁∈ [2⁻², 2¹²]和λ₂∈ [2⁻⁵, 2²].

为了比较离群值检测的性能，我们使用了马修斯相关系数（MCC），定义为

电动机控制中心 = \frac{TP（转移定价） \cdot TN公司 - FP公司 \cdot FN公司}{\sqrt{(TP（转移定价） + FP公司) (TP（转移定价） + FN公司) (TN公司 + FP公司) (TN公司 + FN公司)}} .

这里，TP代表真阳性，其中检测到的异常区域确实是异常值；TN代表真阴性，其中标记的正常区域是正常的；FP代表假阳性，其中检测到的异常区域实际上是正常的；FN代表假阴性，其中标记的正常区域实际上是一个异常值。MCC值越大越好，其中MCC=1表示分类器越完善，MCC=0表示随机猜测。我们根据提议的GLMM、GLMM-CAR、Scan Statistic和Oracle评估了MCCγ。MCC见图1.该方法的MCC与oracle的MCC相当γ。它在增加的基础上有所改善K（K）并在异常值比例增加时趋于稳定。相反，当异常值比例或K（K）增加。我们进一步给出了该模型的真阳性率（敏感性）和真阴性率（特异性）附录D.1的补充材料.图S1和S2系列结果表明，当异常值比例增加时，GLMM和GLMM-CAR的灵敏度均较低。这一发现与现有文献一致，例如，她和欧文（2011）这表明基于cutoff的离群值检测在多个离群值情况下可能无法很好地运行，因为最佳的cutoff选择取决于通常未知的真实残差分布。扫描统计的MCC约为零，即使在n个和K（K）增加，表明扫描统计无法检测到多个异常区域。总之，我们的方法在识别离群值方面表现出了良好的性能，特别是当离群值的比例增加时。

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

在单独的窗口中打开

图1。

MCC，在1000次复制中改变异常值的数量。

我们进一步将该方法与GLMM、GLMM-CAR和Oracle进行了比较α根据个体水平协变量效应的偏差 $({\hat{α}}_{1} - α_{1})$ 以及社区层面的协变量效应 $({\hat{α}}_{2} - α_{2})$ ; 的经验覆盖概率 ${\hat{α}}_{1}$ 和 ${\hat{α}}_{2}$ ; 以及区域级患病率的均方根误差（RMSE） $\sqrt{\sum_{我 = 1}^{K（K）} {({\hat{对}}_{我} - 对_{我})}^{2} / K（K）}$ .给， $对_{我} = E类 ({Y（Y）}_{我 j个} ∣ {X（X）}_{我})$ 和 ${\hat{对}}_{我}$ 是估计的个体水平流行率估计值的经验平均值， ${\hat{对}}_{我 j个}$ ，已接管j个= 1, …,n个_我结果见表1扫描统计被排除在比较之外，因为它没有提供α和对_我.估算偏差α₁，在所提出的方法中，个体水平的协变量效应接近于零，特别是当两者同时存在时n个和K（K）增加。不同方法的置信区间具有可比性。的偏见 ${\hat{α}}_{2}$ ，区域层面的协变量效应被降低为K（K）在建议的方法中增加。所提方法和Oracle的性能α在偏见方面是相似的。覆盖率略低α₂当K（K）=20可能是由于相对较小的K（K）用于估计大量区域级参数(β和γ). 当K（K）增加到40。RMSE ${\hat{对}}_{我}$ 虽然正如预期的那样，与Oracle相比，所提议的方法比两个GLMM中的要小α.

表1

偏差（±标准误差）和经验覆盖概率 ${\hat{α}}_{1}$ 和 ${\hat{α}}_{2}$ ，RMSE的 ${{\hat{对}}_{我}}_{我 = 1}^{K（K）}$ ，超过1000次复制，改变了真实异常区域的比例(K（K）_{O（运行）}/K（K）).

K（K）_{O（运行）}/K（K）	方法	n个=每个地区50					n个=每个地区100
		${\hat{α}}_{1}$		${\hat{α}}_{2}$		RMSE的 ${{\hat{对}}_{我}}$	${\hat{α}}_{1}$		${\hat{α}}_{2}$		RMSE的 ${{\hat{对}}_{我}}$
		偏见	人物配对关系	偏见	人物配对关系	RMSE的 ${{\hat{对}}_{我}}$	偏见	人物配对关系	偏见	人物配对关系	RMSE的 ${{\hat{对}}_{我}}$
K（K）=20个地区
	建议的内容	−.004 ± .008	.942	.002 ± .011	.909	.053	.003 ± .006	.950	−.011 ± .009	.871	.040
0%	单位：百万加仑	−.004 ± .008	.950	−.008 ± .013	.809	.057	.003 ± .006	.951	−.013 ± .011	.684	.044
	GLMM-CAR公司	−.004 ± .008	.957	−.013 ± .011	.943	.052	.002 ± .006	.956	−.014 ± .009	.956	.038
	甲骨文公司α	−.002± .007	.948	.003 ± .006	.958	.016	.003±.005	.953	−.003 ± .005	.953	.011
	提出	−.006 ± .008	.949	−.004 ± .012	.931	.054	.002 ± .006	.957	−.010 ± .010	.898	.041
5%	GLMM公司	−.006 ± .008	.951	−.013 ± .016	.663	.062	.001 ± .006	.953	−.026 ± .017	.465	.046
	GLMM-CAR公司	−.003 ± .008	.955	−.013 ± .016	.912	.062	.001 ± .006	.951	−.017 ± .015	.940	.045
	甲骨文公司α	−.003±.007	.944	.006 ± .006	.959	.016	.003 ± .005	.953	−.003 ± .005	.945	.011
	提出	−.003±.008	.936	−.002 ± .013	.943	.055	.004 ± .006	.950	−.011 ± .010	.918	.041
10%	GLMM公司	−.003 ± .008	.936	−.012±.020	.608	.063	.003 ± .006	.951	−.025 ± .020	.466	.046
	GLMM-CAR公司	−.005 ± .008	.951	.002 ± .019	.920	.063	.002 ± .006	.956	−.006 ± .019	.934	.046
	甲骨文公司α	−.001 ± .007	.943	.005 ± .007	.961	.016	.003 ± .005	.953	−.004 ± .005	.951	.011
	提出	−.001 ± .009	.937	−.003 ± .014	.947	.056	.003 ± .006	.961	−.010 ± .010	.907	第041页
15%	GLMM公司	−.002 ± .009	.930	−.012 ± .023	.545	.063	.002 ± .006	.956	−.025 ± .023	.400	.046
	GLMM-CAR公司	−.005 ± .008	.948	.010 ± .023	.911	.063	.002 ± .006	.955	.002 ± .023	.924	.046
	甲骨文公司α	.001 ± .007	.950	.005 ± .007	.954	.016	.002±.005	.956	−.004 ± .005	.940	.011
K（K）=40个地区
	提出	−.004 ± .006	.948	.010 ± .007	.919	.043	−.000 ± .004	.946	−.002±.005	.923	.033
0%	GLMM公司	−.004 ± .006	.948	.004 ± .009	.804	.055	−.001 ± .004	.950	−.003 ± .008	.683	.043
	GLMM-CAR公司	−.003 ± .006	.943	.002 ± .007	.961	.047	−.001 ± .004	.949	−.001 ± .005	.983	.035
	甲骨文公司α	−.005 ± .005	.953	.005 ± .005	.950	.012	.000 ± .003	.950	−.001 ± .003	.939	.008
	提出	−.003 ± .006	.938	.010 ± .008	.948	.045	−.000 ± .004	.950	.000 ± .005	.942	.033
5%	GLMM公司	−.004 ± .006	.934	.005 ± .011	.639	.061	−.001 ± .004	.950	−.001±.011	.507	.046
	GLMM-CAR公司	−.004 ± .006	.948	.002 ± .010	.947	.059	−.000 ± .004	.945	−.004 ± .010	.944	.044
	甲骨文公司α	−.004 ± .005	.948	.004 ± .005	.938	.012	.000 ± .003	.951	−.000 ± .003	.932	.008
	建议的内容	−.006 ± .006	.941	.011 ± .008	.963	.046	−.001 ± .004	.943	.003 ± .006	.950	.034
10%	GLMM公司	−.006 ± .006	.941	.004 ± .014	.589	.062	−.002 ± .004	.943	−.003 ± .014	.423	.046
	GLMM-CAR公司	−.004 ± .006	.949	−.002 ± .013	.944	.062	−.002 ± .004	.943	−.007 ± .013	.943	.045
	甲骨文公司α	−.006 ± .005	.951	.006 ± .005	.935	.012	−.001 ± .003	.950	.000±.003	.946	.008
	提出	−.007 ± .006	.945	.010 ± .008	.966	.047	.000 ± .004	.947	.001 ± .006	.951	.035
15%	GLMM公司	−.006±.006	.948	.006 ± .016	.512	.063	−.001 ± .004	.949	.005 ± .017	.383	.046
	GLMM-CAR公司	−.004 ± .006	.950	−.002 ± .016	.950	.063	−.001 ± .004	.945	−.005 ± .016	.945	.045
	甲骨文公司α	−.006 ± .005	.950	.005±.005	.948	.012	.000 ± .003	.952	−.000 ± .003	.955	.008

在单独的窗口中打开

然后，我们比较了β, $\sqrt{\sum_{我 = 1}^{K（K）} {({\hat{β}}_{我} - β_{我})}^{2} / K（K）}$ 根据建议的方法、GLMM-CAR和Oracleβ。其余方法未包括在比较中，因为它们没有提供β.图2显示的RMSE $\hat{β}$ 减少的时间n个或K（K）增加，并且随着离群值比例的增加而略有增加。该方法与Oracle类似β表现优于GLMM-CAR。这表明所提出的方法提供了对基线肥胖率的良好估计。

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

在单独的窗口中打开

图2。

RMSE的 $\hat{β}$ ，在1000次复制中改变异常值的数量。

在Web附录D.2我们报告了不同方法的效果Y（Y）可能不见了。正如预期的那样，在存在缺失的情况下，估计系数的偏差较大，但所提出的方法总体上优于竞争对手。

5.应用于PHINEX数据库

我们将人口普查区块组视为地理单元，并根据行为风险因素监测系统的指导方针，排除了某些样本量较小的区块组(CDC，2016年). 个体协变量Z轴_ij公司包括性别、截至2012年的年龄、种族/民族、保险状况和地区层面的协变量X（X）_我包括城市化和EHI。年龄分为3组：2-4岁、5-9岁和10-14岁。种族和民族被合并成一个单一的协变量，并被分为4组：西班牙裔、非西班牙籍白人、非西班牙裔黑人和非西班牙裔其他人。包括有商业健康服务付款人或医疗补助的患者，并排除了一些没有保险的受试者。城市化从1到11分为3类：城市（1-4）、郊区（5-8）或农村（9-11）。我们对EHI进行了标准化，以确保所提算法的数值稳定性。儿童肥胖症的罕见频率见Web附录E.

职位S公司_我由质心的经纬度矢量定义我-第th块组。我们建造了 $ρ_{我_{1}, 我_{2}}$ 作为测地距离的倒数， $ρ_{我_{1}, 我_{2}} = 1 / {d日}^{地理} ({S公司}_{我_{1}}, {S公司}_{我_{2}})$ ，其中 ${d日}^{地理} ({S公司}_{我_{1}}, {S公司}_{我_{2}})$ 表示之间的较大圆距离 ${S公司}_{我_{1}}$ 和 ${S公司}_{我_{2}}$ 。对于我₁-第个区域，我们保留了L（左）最大的 $ρ_{我_{1}, 我_{2}} 秒$ 并将其他值截断为零，我们在这里处理L（左）作为调谐参数。上的网格搜索λ₁∈ [2⁻¹, 2¹⁷]，λ₂∈ [2⁻⁵, 2²]和L（左）∈{3,5,7}是为了找到使BIC*最大化的最佳调谐参数组合。每个参数的置信区间是使用超过1000次重复的bootstrap构建的。

估计 $\hat{α} 秒$ 总结如下表2总体而言，所提出的方法比GLMM和GLMM-CAR具有更宽的置信区间。考虑到我们的方法有更多的参数需要估计，这一结果是可以预期的，这可能导致更高的变异性。除郊区效应外，我们模型的估计系数与GLMM和GLMM-CAR的估计系数相当。女性的肥胖率低于男性，而年龄较小的儿童肥胖率较低。非拉美裔白人和其他非拉美籍人群的肥胖率均低于非拉美族黑人和拉美裔患者。与商业保险的受试者相比，医疗补助受试者的肥胖患病率更高。EHI与估计的肥胖率呈正相关。

表2

协变量效应的拟合系数和置信区间（括号内）。

	模型
	提出	GLMM公司	GLMM-CAR公司
个体水平协变量
性别（基数：女性）
~ 男性	.235 (.123, .347)	.226 (.147, .305)	.227 (.182, .274)
2012年年龄（基数：学龄前）
~ 学龄儿童	.568 (.449, .687)	.562 (.448, .675)	.558 (.489, .629)
~青少年	.875 (.773, .977)	.869 (.757, .981)	.864 (.798, .930)
种族/民族（基础：白人、非西班牙裔）
~ 黑人、非西班牙裔	.437 (.269, .605)	.434 (.315, .553)	.440 (.360, .519)
~ 其他，非西班牙裔	.035 (−.198, .269)	.042 (−.107, .191)	.054（−.061，.166）
~ 西班牙裔	.680 (.538, .822)	.667 (.556, .779)	.672 (.609, .734)
保险状态（基数：商业）
~ 医疗补助	.522 (.417, .627)	.509 (.408, .611)	.509 (.452, .567)
社区水平协变量
城市化（基础：城市）
~ 郊区的	−.134 (−.217, −.051)	.037 (−.058, .132)	−.030 (−.071, .106)
~ 农村	.125 (−.011, .262)	.237 (.095, .379)	.090 (−.084, .279)
经济困难指数（标准化）	.120 (.083, .156)	.143 (.105, .181)	.096（.044、.149）

在单独的窗口中打开

拟合基线肥胖率， ${罗吉特}^{- 1} ({\hat{β}}_{我})$ ，的左上角显示图3这似乎与大麦迪逊地区的经验知识相吻合。发病率最低的地区包括麦迪逊西部、米德尔顿和维罗纳地区。众所周知，这些地区是最近开发和扩大的，包括一般较年轻的人，与周边地区相比，他们在社会经济上更具优势。中等流行区包括麦迪逊大中部和东部地区，是该地区较为成熟的历史性地区，众所周知，该地区的中产阶级公民较为定型。发病率最高的地区显然是距离麦迪逊市中心地理位置最远的地区，也都在包括麦迪逊在内的戴恩县以外。

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

在单独的窗口中打开

图3。

儿童肥胖监测中的估计基线患病率和确定的异常值。每个多边形代表一个普查区块组。左上：建议方法的结果。异常值标记为趋势上方的黑色（黄色）， ${\hat{γ}}_{我} > 0$ （低于趋势， ${\hat{γ}}_{我} < 0$ ). 右上角：GLMM的结果。异常值标记为趋势上方的黑色（黄色）， ${\hat{第页}}_{我} > 2.5 \hat{σ}$ （低于趋势， ${\hat{第页}}_{我} < 2.5 \hat{σ}$ ). 左下：GLMM-CAR的结果。异常值的标记方式与GLMM相同。右下：扫描统计发现的聚类具有最高的似然比，低于趋势。

该方法识别了几个异常值。肥胖率高于趋势的异常地点( ${\hat{γ}}_{我} > 0$ )和低于趋势( ${\hat{γ}}_{我} < 0$ )分别显示为黑色和黄色，单位为图3我们确定6%的区块组为高于趋势的异常值，8%为低于趋势的异常值。结果显示于表3包括：

粗略肥胖率， ${\hat{对}}_{我}^{原油} = \frac{1}{\sum_{j个} {w个}_{我 j个} 我 ({R（右）}_{我 j个} = 1)} \sum_{j个} {w个}_{我 j个} 我 ({R（右）}_{我 j个} = 1) {Y（Y）}_{我 j个}$ ;
基线肥胖率， ${\hat{对}}_{我}^{bsl公司} = {罗吉特}^{- 1} ({\hat{β}}_{我})$ ;
根据协变量和异常值调整的肥胖率，
${\hat{对}}_{我}^{形容词} = {\hat{E类}}_{γ_{我} = 0} ({Y（Y）}_{我 j个} ∣ {X（X）}_{我}) = \frac{1}{\sum_{j个} {w个}_{我 j个} 我 ({R（右）}_{我 j个} = 1)} \sum_{j个 = 1}^{{n个}_{我}} {w个}_{我 j个} 我 ({R（右）}_{我 j个} = 1) \cdot {罗吉特}^{- 1} ({Z轴}_{我 j个}^{T型} {\hat{α}}_{1} + {X（X）}_{我}^{T型} {\hat{α}}_{2} + {\hat{β}}_{我}),$
和检测频率B类=1000个引导复制。我们注意到，异常值标识为相对的以适应潮流。例如，区块组212的估计粗肥胖率为普通水平（0.180）。然而，粗略比率远高于预期肥胖患病率的拟合值（0.085）。存在无法解释的信息，可能导致发病率升高。因此，它被宣布为高于趋势的异常值。基于bootstrap的检测频率提供了畸变不确定性的一瞥。趋势上方的异常区域的频率往往高于趋势下方的异常区域。

表3

通过该方法识别的异常数据块组的匿名ID、样本大小、粗肥胖率、拟合基线肥胖率、调整肥胖率和检测频率B类=1000个引导复制。

块组ID	未加权样本量(n个_我)	粗肥胖率 $({\hat{对}}_{我}^{原油})$	合适的基线肥胖率 $({\hat{对}}_{我}^{bsl公司})$	调整后的肥胖率 $({\hat{对}}_{我}^{形容词})$	检测频率 $(\sum_{我} 我 ({\hat{γ}}_{我} \neq 0) / B类)$
高于趋势
7	60	.432	.097 (.061, .150)	.264 (.153, .351)	.701
23	91	.234	.048 (.042, .057)	.129 (.110, .167)	.483
24	93	.206	.048（.042，.057）	.113 (.095, .128)	.515
25	104	.207	.048 (.042, .057)	.115 (.094, .139)	.496
83	96	.291	.058 (.047, .064)	.174 (.133, .179)	.689
85	91	.356	.058 (.047, .064)	.187（.139，.196）	.895
100	62	.288	.058 (.047, .064)	.149 (.117, .155)	.763
102	53	.335	.058 (.047, .064)	.217 (.147, .246)	.579
124	66	.207	.058 (.047, .064)	.117 (.090, .136)	.535
200	93	.218	.058（.047，.064）	.133 (.102, .151)	.474
212	100	.180	.048 (.042, .057)	.085 (.076, .094)	.634
244	71	.203	.053 (.044, .064)	.121 (.096, .135)	.446
245	94	.248	.053 (.044, .064)	.148 (.105, .184)	.573
252	68	.278	.056 (.046, .073)	.148 (.108, .179)	.706
254	82	.204	.058 (.048, .071)	.103 (.080, .111)	.698
264	74	.257	.048 (.042, .057)	.130 (.092, .157)	.726
低于趋势
22	110	.063	.048 (.042, .057)	.123 (.109, .151)	.363
32	221	.045	.048 (.042, .057)	.092 (.080, .101)	.079
35	146	.067	.048（.042，.057）	.126 (.109, .151)	.300
70	67	.168	.058 (.047, .064)	.283 (.203, .292)	.269
73	109	.175	.058 (.047, .064)	.273 (.202, .279)	.136
82	259	.173	.058 (.047, .064)	.329 (.209, .342)	.106
94	134	.105	.058 (.047, .064)	.198 (.159, .210)	.383
115	68	.061	.058 (.047, .064)	.157 (.120, .161)	.468
118	192	.141	.058（.047，.064）	.248 (.155, .253)	.062
125	81	.047	.058 (.047, .064)	.115 (.091, .115)	.277
127	295	.051	.058 (.047, .064)	.114 (.094, .121)	.195
128	125	.054	.058 (.047, .064)	.154 (.122, .175)	.698
136	466	.038	.048 (.042, .057)	.114 (.100, .132)	.684
153	469	.045	.048 (.042, .057)	.095 (.085, .103)	.047
159	202	.081	.058（.047，.064）	.158 (.130, .176)	.345
168	139	.117	.056 (.046, .071)	.221 (.169, .247)	.418
186	66	.097	.058 (.047, .064)	.167 (.118, .184)	.204
229	133	.139	.077 (.048, .127)	.299（.189、.395）	.599
235	61	.139	.077 (.048, .127)	.235 (.139, .312)	.365
246	210	.060	.053 (.044, .064)	.138 (.106, .158)	.283
262	90	.076	.081 (.050, .114)	.189 (.116, .237)	.459

在单独的窗口中打开

由于检测异常区域的基本机制不同，因此通过各种方法识别的异常值并不完全重叠。粗略方法确定比例最高/最低的区域，而不考虑协变量效应。GLMM方法基于不排除平滑区域效应的最高/最低残差效应来检测异常值。另一方面，该方法和GLMM-CAR根据平滑区域效应进行调整后的剩余区域效应将区域声明为异常值。与我们的方法不同，GLMM和GLMM-CAR都识别出很少的异常值：两个区域由GLMM识别，一个区域由GLMM-CAR识别。如中所述第4节，这种保守行为与预期相符。另一个值得注意的是，扫描统计数据仅确定麦迪逊中西部地区为异常低。

我们模型中的局部疫情可以在颗粒水平上进行比较研究。肥胖患病率由患者的人口统计学特征、行为和社区环境因素的相互作用决定。我们的模型只考虑了其中的一个子集。离群值可能代表环境与预期存在显著差异的社区（例如，与杂货店、公园等的平均通行情况相比，环境要好或差得多），和/或它可能代表社区成员的行为与预期有很大不同（例如，更多或更少的体力活动，更好或更差的饮食习惯，等等）。根据我们的研究结果，医疗保健专业人员可以研究异常值中的风险因素，并将这些因素与相邻的区块组进行比较。

6.结束语

基于使用常规收集的EHR数据进行儿童肥胖监测的动机，我们开发了一个多级惩罚logistic回归模型，其中融合了融合和非凸稀疏惩罚，以同时进行区域平滑和异常值检测。虽然我们只考虑了空间监视，但我们有兴趣将该方法推广到用于时空监视的纵向数据设置。

在我们的论文中，我们假设BMI是MAR，这在观测数据中是不可测试的。EHR数据MAR的可行性是一个活跃的研究领域(Snyder等人，2018年). 未来，我们可以开发敏感性分析技术，研究结果对非受控混杂的敏感性(格陵兰，2004年). 另一个未来的方向是为所建议的工作开发原则推理程序。我们可以潜在地使用推理程序中的思想来处理惩罚广义线性模型(Taylor和Tibshirani，2018年).

补充材料

单击此处查看。^{（417K，pdf）}

确认

我们感谢联合主编、副主编和所有审稿人提出的有帮助和有见地的意见，这些意见帮助改进了手稿。这项工作得到了美国国立卫生研究院授予的R21HD086754、P30CA015704和S10OD020069、韩国国立研究基金会授予的2020R1G1A1A01006229和韩国首尔女子大学研究拨款授予的1-1903-2017的支持。作者感谢Donghyeon Yu（韩国仁川仁川英华大学）共享了优化最小化算法的R源代码，感谢Albert Y.Kim（马萨诸塞州阿默斯特市阿默斯特学院）共享了R（右）包裹空间Epi.

脚注

支持信息

中引用的Web附录、表和图第3节,4，以及5可在威利在线图书馆的生物统计学网站上获取此论文。

数据可用性声明

由于道德限制，支持本文研究结果的数据尚未公开。感兴趣的读者可以联系威斯康星大学IRB委员会(https://www.medicine.wisc.edu/research/uw-health-sciences-irbs)了解详细信息。

工具书类

Besag J、York J和Molli A（1991年）。贝叶斯图像恢复及其在空间统计中的两个应用.统计数学研究所年鉴,43(1):1–20.[谷歌学者]
CDC（2016）。行为风险因素监测系统，2015年BRFSS数据的可比性（版本#1-修订日期：2016年6月）技术报告。[谷歌学者]
Chernozhukov V、Hansen C和Liao Y（2017）。对密集和稀疏信号总和恢复的熔岩攻击.统计年鉴,45(1):39–76.[谷歌学者]
Efron B和Tibshirani RJ（1994）。引导程序简介查普曼和霍尔/CRC。[谷歌学者]
ESRI（2012）。环境系统研究所挂毯分割参考指南.CA：雷德兰。[谷歌学者]
Farrington CP、Andrews NJ、Beale AD和Catchpole MA（1996年）。传染病暴发早期检测的统计算法.英国皇家统计学会杂志。A系列（社会统计学）,159(三):547–563.[谷歌学者]
Flood TL、Zhao Y-Q、Tomayko EJ、Tandia A、Carrel AL和Hanrahan LP（2015）。儿童肥胖的电子健康记录和社区健康监测.美国预防医学杂志,48(2):234–240.[PMC免费文章][公共医学][谷歌学者]
Friedman DJ、Parrish RG和Ross DA（2013年）。电子健康记录与美国公共卫生：现状与未来展望.美国公共卫生杂志,103(9):1560–1567.[PMC免费文章][公共医学][谷歌学者]
Ghosh M、Natarajan K、Waller LA和Kim D（1999年）。空间数据分析的层次Bayes{GLMs}：疾病制图应用.统计规划与推断杂志,75(2)：305–318。[谷歌学者]
格陵兰岛S（2004）。先前分布对非受控混杂和反应偏差的影响.美国统计协会杂志,98(461):47–54.[谷歌学者]
Guilbert TW、Arndt B、Temte J、Adams A、Buckingham W、Tandia A、Tomasallo C、Anderson HA和Hanrahan LP（2012年）。UW ehealth-PHINEX临床电子健康记录公共健康信息交换的理论与应用.WMJ公司,111(三):124–33. [公共医学][谷歌学者]
Jung I（2009）。协变量调整空间扫描统计的广义线性模型方法.医学统计学,28(7):1131–1143. [公共医学][谷歌学者]
Kim S-J、Koh K、Boyd S和Gorinevsky D（2009年）。L1趋势过滤.SIAM审查,51(2):339–360.[谷歌学者]
Kulldorff M和Nagarwalla N（1995年）。空间疾病集群：检测和推断.医学统计学,14(8):799–810. [公共医学][谷歌学者]
Lee S、Kwon S和Kim Y（2016年）。惩罚优化问题的改进局部二次逼近算法.计算统计与数据分析,94:275–286.[谷歌学者]
Little RJ和Rubin DB（2014年）。缺失数据的统计分析约翰·威利父子公司。[谷歌学者]
Longjohn M、Sheon AR、Card-Higginson P、Nader PR和Mason M（2010年）。从国家儿童肥胖监测中学习.卫生事务,29(三):463–472. [公共医学][谷歌学者]
Mercer LD、Wakefield J、Pantazis A、Lutambi AM、Masanja H和Clark S（2015）。复杂调查数据的时空平滑：儿童死亡率的小区域估计.应用统计学年鉴,9(4):1889–1905.[PMC免费文章][公共医学][谷歌学者]
Nathan RP和Adams CF（1989）。关于城市困境的四种观点.政治科学季刊,104(三):483.[谷歌学者]
Panczak R、Held L、Moser A、Jones PA、Rühli FJ和Staub K（2016）。寻找大人物：瑞士男性义务兵肥胖的小面积绘图和空间建模.BMC肥胖,三(1):1–12.[PMC免费文章][公共医学][谷歌学者]
Puth MT、Neuhäuser M和Ruxton GD（2015年）。关于用bootstrapping计算置信区间的各种方法.动物生态学杂志,84(4):892–897. [公共医学][谷歌学者]
She Y和Owen AB（2011年）。基于非凸惩罚回归的异常检测.美国统计协会杂志,106(494)：626–639。[谷歌学者]
Snyder JW、Bauer CR、Beaulieu-Jones BK、Pendergrass SA、Lavage DR和Moore JH（2018年）。描述和管理电子健康记录中缺失的结构化数据：数据分析.JMIR医学信息学,6(1)：e11。[PMC免费文章][公共医学][谷歌学者]
Taylor J和Tibshirani R（2018年）。l1-惩罚似然模型的选择后推理.加拿大统计杂志,46(1):41–61.[PMC免费文章][公共医学][谷歌学者]
Tibshirani R、Saunders M、Rosset S、Zhu J和Knight K（2005）。通过融合套索实现轻盈流畅.英国皇家统计学会杂志。B系列：统计方法,67:91–108.[谷歌学者]
Tibshirani RJ（2014）。基于趋势滤波的自适应分段多项式估计.统计年鉴,42(1):285–323.[谷歌学者]
Tibshirani RJ和Taylor J（2011）。广义套索的解路径.统计年鉴,39(三):1335–1371.[谷歌学者]
Ugarte MD、Goicoa T和Militino AF（2010年）。基于惩罚样条的死亡风险时空建模.环境计量学,21(3–4):270–289.[谷歌学者]
Wang H和Rodríguez A（2014）。使用对数线性模型和广义拉索惩罚确定佛罗里达州的儿童癌症簇.统计与公共政策,1(1):86–96.[PMC免费文章][公共医学][谷歌学者]
Yu D、Won J-H、Lee T、Lim J和Yoon S（2015）。基于优化最小化和并行处理的高维融合拉索回归.计算与图形统计杂志,24(1):121–153.[谷歌学者]
Zhao Y、Zeng D、Herring AH、Ising A、Waller A、Richardson D和Kosorok MR（2011年）。使用局部时空方法检测疾病暴发.生物计量学,67(4):1508–1517.[PMC免费文章][公共医学][谷歌学者]
Zou H、Hastie T和Tibshirani R（2007年）。关于套索的“自由度”.统计年鉴,35(5):2173–2192.[谷歌学者]

基于惩罚回归的同步空间平滑和离群值检测，应用于电子健康记录中的儿童肥胖监测

Young-Geun Choi先生

劳伦斯·P·汉拉罕

德里克·诺顿

赵应琪

关联数据

总结：

1.简介

2.数据

3.方法

3.1、。模型设置

3.2. 使用完整数据进行估算

3.3. 优化算法

3.4. 调谐参数的选择

3.5. 考虑缺失和选择偏差的权重

4.模拟研究

表1

5.应用于PHINEX数据库

表2

表3

6.结束语

补充材料

补充材料

确认

脚注

数据可用性声明

工具书类