Bayesian kernel machine regression for count data: modelling the association between social vulnerability and COVID-19 deaths in South Carolina¶

作者: Fedelis Mutiso, Hong Li, John L Pearce, Sara E Benjamin-Neelon, Noel T Mueller et al.
来源: Journal of the Royal Statistical Society Series C
主题: 流行病学
相关性: 4/10
机构绿灯: University of California, Davis（US News 前 50，免分进入精读）
链接: https://doi.org/10.1093/jrsssc/qlad094

一、领域脉络与小综述¶

这个方向是什么¶

这个子方向是 “暴露-响应关系的半参数建模”，核心问题是：在流行病学或环境健康研究中，如何灵活地估计一组高维、相关、可能存在非线性和交互作用的暴露变量（如污染物、社会脆弱性指标）对健康结局（如死亡率、发病率）的联合影响，同时避免“维数灾难”和“模型误设”带来的偏差。当前成熟度属于“方法已提出、正在向特定数据类型（如计数数据、纵向数据）扩展”的阶段。

发展脉络（history）¶

奠基工作：核机器回归 (KMR) 与贝叶斯版本 (BKMR)
- Bobb et al. (2015, Biostatistics)：提出了 Bayesian Kernel Machine Regression (BKMR)，这是该方向的基石。它将暴露变量通过一个核函数映射到高维特征空间，从而自动捕捉非线性和交互效应，并用贝叶斯方法进行推断。作者引用句：“Bobb et al. (2015) proposed Bayesian kernel machine regression (BKMR) to model the joint effect of a mixture of exposures on a continuous outcome.” 它留下了两个主要口子：① 只处理连续结局；② 未考虑空间或时间相关性。
主要进展：向不同结局类型和复杂数据结构的扩展
- Coker et al. (2016, Environmental Health)：将 BKMR 应用于环境混合物与出生结局的研究，验证了其在真实流行病学数据中的实用性，但仍是连续结局。
- Valeri et al. (2017, Environmental Health Perspectives)：将 BKMR 扩展到二元结局（如早产），通过 probit 链接函数处理。作者引用句：“Valeri et al. (2017) extended BKMR to binary outcomes using a probit link.” 这解决了“非连续结局”的第一个口子，但留下了“计数数据”和“空间相关性”的口子。
- Antonelli et al. (2020, Biometrics)：提出了 BKMR with variable selection，通过 spike-and-slab 先验识别混合物中最重要的变量。作者引用句：“Antonelli et al. (2020) extended BKMR to include variable selection via spike-and-slab priors.” 这解决了“变量重要性识别”的问题，但未处理计数数据或空间结构。
当前 Frontier：处理计数数据、时空异质性与高维暴露
- 本文 (Mutiso et al., 2024, JRSS-C)：将 BKMR 扩展到计数数据（负二项分布），并整合了空间随机效应和时间平滑函数，以处理 COVID-19 死亡率的时空异质性。作者引用句：“In this paper, we extend BKMR to count data by proposing a negative binomial BKMR model that can accommodate overdispersion, spatial effects, and temporal trends.” 这是对 BKMR 框架的一个自然且必要的扩展，填补了其在流行病学中常见的“计数+时空”数据场景的空白。

子线索聚类¶

这些被引文献大致落在两条子线索上：

方法学扩展线索：从连续结局 (Bobb 2015) → 二元结局 (Valeri 2017) → 计数结局 (本文)。这条线索的核心是改变链接函数和似然，并调整相应的贝叶斯计算算法。
模型功能增强线索：从基础 BKMR (Bobb 2015) → 变量选择 (Antonelli 2020) → 空间/时间效应 (本文)。这条线索的核心是在核函数之外，加入额外的模型组件（如随机效应、平滑项）来处理数据结构的复杂性。

这个方向在追问的核心问题¶

如何高效识别高维暴露混合物中的关键变量及其交互作用？ 当前主流方法是 BKMR 结合 spike-and-slab 先验 (Antonelli 2020) 或后验 inclusion probability。瓶颈在于计算复杂度随暴露维度增加而急剧上升。
如何将 BKMR 扩展到更复杂的结局类型（如零膨胀计数、生存数据）？ 本文解决了计数数据，但零膨胀和生存数据仍是开放问题。
如何在 BKMR 框架下有效处理空间和时间相关性？ 本文通过加入空间随机效应和平滑时间函数来解决，但这是否是最优的（例如，是否可以通过核函数本身捕捉时空结构）仍有待探索。

⚠️ 作者的 framing¶

作者把缺口 frame 成什么：作者将缺口 frame 为“BKMR 尚未被扩展到计数数据，且未同时处理时空异质性”。因此，本文被定位为“显然的下一步”——一个将 BKMR 应用于 COVID-19 死亡率这一典型计数+时空数据场景的完整解决方案。
哪些竞争路线被他淡化或回避了：作者淡化了广义加性模型 (GAM) 和多元自适应回归样条 (MARS) 等非参数方法。这些方法也能处理非线性和交互，但作者在引言中仅用一句话带过，强调 BKMR 在捕捉“复杂交互”方面的优势，而未深入比较。此外，作者回避了深度学习方法（如神经网络），这些方法在捕捉高维交互方面可能更强大，但可解释性差。
什么明显该被引 / 该存在、却没出现在 intro 里？：作者没有引用任何关于 “因果推断” 的文献。虽然 BKMR 通常被用于关联研究，但 COVID-19 死亡率与社会脆弱性之间存在明显的混杂因素（如医疗资源、政策干预）。将 BKMR 与因果推断框架（如 g-computation、IPW）结合，以估计“脆弱性”的因果效应，是一个明显的、未被提及的下一步。（值得研究者去查的问题：BKMR 在因果推断中的应用现状如何？）

张力¶

未见明显对立引用。所有被引工作都在逐步扩展 BKMR 的适用性，方向一致。

二、最核心、最简单的例子 / 数学问题¶

第一步：把符号、模型、可观测数据交代清楚¶

符号：
- \( Y_i \)：第 \( i \) 个县（county）在某个时间点的 COVID-19 死亡人数。这是可观测的响应变量，是一个计数（非负整数）。
- \( \mathbf{z}_i \)：第 \( i \) 个县的 \( M \) 维暴露向量，由 Social Vulnerability Index (SVI) 的 15 个普查变量组成。这是可观测的暴露变量。
- \( \mathbf{x}_i \)：第 \( i \) 个县的 \( p \) 维协变量向量（如人口密度、中位年龄）。这是可观测的协变量。
- \( s_i \)：第 \( i \) 个县的空间位置（如经纬度或区域索引）。这是可观测的空间信息。
- \( t_i \)：第 \( i \) 个观测的时间点（如月份）。这是可观测的时间信息。
- \( \mu_i \)：第 \( i \) 个县的期望死亡人数，即 \( E[Y_i | \mathbf{z}_i, \mathbf{x}_i, s_i, t_i] \)。这是要建模的 estimand。
- \( h(\mathbf{z}_i) \)：一个未知的、非线性的“暴露-响应”函数，由核机器回归捕捉。这是核心的未知函数。
- \( \beta \)：协变量 \( \mathbf{x}_i \) 的回归系数向量。这是要估计的参数。
- \( \theta(s_i) \)：第 \( i \) 个县的空间随机效应。这是潜在变量，用于捕捉空间相关性。
- \( f(t_i) \)：时间平滑函数（如样条基函数的线性组合）。这是要估计的平滑函数。
- \( r \)：负二项分布的分散参数（overdispersion parameter）。这是要估计的 nuisance 参数。
模型：
- 数据生成机制：\( Y_i \) 服从负二项分布，即 \( Y_i \sim \text{NB}(\mu_i, r) \)，其中 \( \mu_i \) 是均值，\( r \) 是分散参数（方差 = \( \mu_i + \mu_i^2 / r \)）。
- 均值模型（链接函数为 log）：
  \[\log(\mu_i) = \mathbf{x}_i^\top \beta + h(\mathbf{z}_i) + \theta(s_i) + f(t_i)\]
- 核心假设：\( h(\cdot) \) 由一个正定核函数 \( K(\cdot, \cdot) \) 生成，并位于一个再生核希尔伯特空间 (RKHS) 中。这意味着 \( h(\cdot) \) 可以表示为 \( h(\mathbf{z}) = \sum_{j=1}^n \alpha_j K(\mathbf{z}, \mathbf{z}_j) \)。在贝叶斯框架下，这等价于对 \( h(\cdot) \) 施加一个高斯过程先验，其协方差函数为 \( K \)。
可观测数据：
- 研究者实际能观测到的是：\( \{ (Y_i, \mathbf{z}_i, \mathbf{x}_i, s_i, t_i) \}_{i=1}^n \)，即每个县的死亡人数、SVI 变量、协变量、位置和时间。
- 想要但观测不到的是：① 真实的暴露-响应函数 \( h(\cdot) \)；② 空间随机效应 \( \theta(s_i) \)；③ 时间平滑函数 \( f(t) \) 的精确形式。这些都需要通过模型假设和贝叶斯推断来“识别”或“估计”。

第二步：讲最小内核¶

本文的最小内核是 “如何将 BKMR 从连续结局扩展到计数结局”。我们剥去空间和时间效应，只考虑一个最简单的特例：没有空间效应、没有时间趋势、没有协变量，只有暴露 \( \mathbf{z}_i \) 和计数结局 \( Y_i \)。

最简特例：假设 \( Y_i \sim \text{Poisson}(\mu_i) \)（即 \( r \to \infty \) 的负二项特例），且模型为：
\[\log(\mu_i) = h(\mathbf{z}_i)\]
其中 \( h(\cdot) \) 服从一个均值为 0、协方差为 \( \tau^2 K(\cdot, \cdot) \) 的高斯过程先验。
核心思路：在连续结局的 BKMR 中，\( Y_i = h(\mathbf{z}_i) + \epsilon_i \)，由于误差 \( \epsilon_i \) 是高斯分布，\( h(\cdot) \) 的后验分布可以通过共轭性（高斯过程先验 + 高斯似然 = 高斯过程后验）解析地得到。但在计数数据中，似然是 Poisson，与高斯先验不共轭，后验无法解析求解。
本文的关键想法：使用 数据增广 (Data Augmentation) 技巧。具体来说，利用 Pólya-Gamma 分布 的性质：对于 Poisson 回归，可以引入一个潜变量 \( \omega_i \)，使得条件似然 \( p(Y_i | \omega_i, h(\mathbf{z}_i)) \) 变成高斯形式。这样，在给定 \( \omega_i \) 的条件下，\( h(\cdot) \) 的后验更新就变成了一个标准的高斯过程回归问题，可以解析计算。然后，再通过 Gibbs 采样从 \( \omega_i \) 的条件后验中采样。
为什么成立：Pólya-Gamma 数据增广技巧（Polson, Scott, & Windle, 2013）提供了一个“共轭性”的桥梁。它将非共轭的 Poisson 回归问题，转化为一个条件共轭的、可迭代采样的高斯过程回归问题。因此，整个贝叶斯推断可以通过一个高效的 Gibbs 采样器完成，而无需使用 Metropolis-Hastings 等低效算法。
一句话总结：本文的核心数学贡献是将 Pólya-Gamma 数据增广技巧应用于 BKMR 框架，从而将 BKMR 从连续结局推广到计数结局，并保持了计算上的可处理性。

三、这篇论文做了什么¶

三句话¶

研究了什么问题：研究了社会脆弱性（由 SVI 的 15 个变量衡量）与 COVID-19 死亡率之间的关联，并开发了一个能够处理计数数据、非线性/交互效应、空间相关性和时间趋势的统计模型。
核心工具 / 方法：提出了 负二项 Bayesian Kernel Machine Regression (NB-BKMR) 模型，该模型结合了核机器回归（捕捉非线性/交互）、负二项似然（处理过离散计数）、空间随机效应（ICAR 先验）和时间平滑函数（B-spline），并使用 Pólya-Gamma 数据增广的 Gibbs 采样器进行推断。
主要结论：模拟研究表明 NB-BKMR 能有效恢复暴露-响应函数、识别关键变量，并优于忽略空间或时间效应的简化模型。在南卡罗来纳州 2021 年 COVID-19 死亡数据应用中，模型识别出“社会经济地位”和“住房类型与交通”是影响死亡率的关键脆弱性维度，并发现脆弱性效应随时间变化。

关键设定与假设¶

在第二节最小记号的基础上，补全完整设定：

完整模型：
\[Y_i \sim \text{NB}(\mu_i, r)\]

\[\log(\mu_i) = \mathbf{x}_i^\top \beta + h(\mathbf{z}_i) + \theta(s_i) + f(t_i)\]
- \( h(\cdot) \)：高斯过程先验，\( h \sim \mathcal{GP}(0, \tau^2 K) \)，核函数 \( K \) 为高斯核（RBF）或 Matern 核。\( \tau^2 \) 是方差参数，控制 \( h \) 的幅度。
- \( \theta(s_i) \)：空间随机效应，采用 条件自回归 (CAR) 先验，特别是 Intrinsic CAR (ICAR) 先验。这意味着 \( \theta(s_i) \) 的条件分布依赖于其相邻县（由邻接矩阵定义）的效应值，从而引入空间平滑。
- \( f(t_i) \)：时间平滑函数，用 B-spline 基函数 的线性组合表示，\( f(t) = \sum_{k=1}^K \gamma_k B_k(t) \)，并对系数 \( \gamma_k \) 施加随机游走先验（如二阶随机游走）以实现平滑。
相比已有文献的强化：
- 相比 Bobb (2015)：将结局从连续扩展到计数（负二项），并加入了空间和时间效应。
- 相比 Valeri (2017)：将结局从二元扩展到计数，并加入了空间和时间效应。
- 相比 Antonelli (2020)：加入了空间和时间效应，但本文没有包含变量选择（spike-and-slab 先验）。这是一个弱化：本文通过后验的“变量重要性”指标（如将某个变量固定后预测变化）来识别关键变量，而非正式的贝叶斯变量选择。
关键假设：
- SUTVA-like 假设：一个县的死亡人数不受其他县暴露水平的影响（除了通过空间随机效应捕捉的相关性）。
- 无未测量混杂：模型中的协变量 \( \mathbf{x}_i \) 和空间/时间效应足以解释暴露 \( \mathbf{z}_i \) 与结局 \( Y_i \) 之间的混杂。这是一个很强的、未明确讨论的假设。
- 核函数选择：高斯核的带宽参数是固定的或通过先验估计的，其选择会影响模型表现。

主要结果¶

模拟研究：
- 设定：模拟了 46 个县（模仿南卡罗来纳州），4 个时间点，15 个 SVI 变量。生成了三种真实暴露-响应函数：线性、非线性和非线性+交互。
- 核心结论：NB-BKMR 模型（包含空间和时间效应）在估计 \( h(\cdot) \) 和预测 \( \mu_i \) 方面，始终优于忽略空间效应、忽略时间效应或两者都忽略的简化模型。具体地，均方根误差 (RMSE) 和偏差显著更低。
- 变量重要性：通过“后验包含概率”（PIP）或“变量固定后的预测变化”指标，模型能正确识别出对结局影响最大的前几个变量。
真实数据应用：南卡罗来纳州 COVID-19 死亡 (2021)
- 数据：46 个县，12 个月，15 个 SVI 变量（分为 4 个主题：社会经济地位、家庭组成与残疾、少数族裔与语言、住房类型与交通），协变量包括人口密度、中位年龄、ICU 床位等。
- 结果：
  1. 关键脆弱性维度：模型发现，“社会经济地位”（如贫困率、失业率）和“住房类型与交通”（如拥挤住房、无车辆家庭）是影响 COVID-19 死亡率的最重要维度。而“少数族裔与语言”和“家庭组成”的重要性相对较低。
  2. 脆弱性效应的时间变化：模型估计的“脆弱性效应”（即 \( h(\mathbf{z}_i) \) 的预测值）在 2021 年初（疫苗推广前）较高，随后下降，但在 Delta 变体流行期间（夏末秋初）再次上升。这提示脆弱性的影响是动态的。
  3. 空间模式：空间随机效应揭示了某些县（如农村县）具有高于平均水平的基线风险，即使控制了 SVI 和协变量。
这个例子想说明什么：① NB-BKMR 能处理真实世界复杂的计数数据；② 它能提供比传统加性模型更丰富的洞察（如非线性、交互、时空变化）；③ 其输出（变量重要性、脆弱性效应图）对公共卫生决策有实际指导意义。

证明路线与技术技巧¶

本文是应用 / 方法型论文，没有严格的渐近理论证明。其“证明”在于贝叶斯推断算法的推导和有效性。

整体路线（Gibbs 采样器）：
1. 数据增广：利用 Pólya-Gamma 分布，为每个观测 \( i \) 引入一个潜变量 \( \omega_i \)。条件于 \( \omega_i \)，负二项似然转化为一个加权高斯似然。
2. 更新 \( h(\cdot) \)：在给定 \( \omega_i \) 和其他参数后，\( h(\cdot) \) 的后验是一个高斯过程，其均值和协方差有闭式解（类似于高斯过程回归）。通过“有限维近似”（如将 \( h \) 表示为基函数展开或使用诱导点）进行采样。
3. 更新空间效应 \( \theta \)：在给定其他参数后，\( \theta \) 的条件后验是一个多元高斯分布，其精度矩阵由 ICAR 先验和似然贡献构成。可以通过稀疏矩阵运算高效采样。
4. 更新时间效应 \( f \)：在给定其他参数后，B-spline 系数 \( \gamma \) 的条件后验是一个多元高斯分布，其先验是随机游走。
5. 更新其他参数：回归系数 \( \beta \)、方差参数 \( \tau^2 \)、分散参数 \( r \) 等，都可以从它们的条件后验中直接采样（通常是共轭分布）。
关键跳跃点：
- 难点：如何将非共轭的负二项似然与高斯过程先验结合。
- 解决办法：Pólya-Gamma 数据增广。这是整个算法的基石，它将一个复杂的非共轭问题转化为一系列条件共轭的、易于采样的高斯问题。
技术技巧点名：
- Pólya-Gamma 数据增广：用于处理负二项似然，是核心技巧。
- ICAR 先验：用于空间随机效应建模，捕捉空间相关性。
- B-spline + 随机游走先验：用于时间平滑函数建模。
- Gibbs 采样器：整个推断框架，通过迭代采样所有参数的条件后验分布。

🔎 结论是否比证明窄¶

是。作者在引言和结论中声称模型能“识别关键 SVI 变量”，但模型本身并未包含正式的变量选择机制（如 spike-and-slab 先验）。所谓的“变量重要性”是通过后验预测来评估的：将某个变量固定在其某个分位数（如中位数），然后观察预测的死亡率变化。这是一种描述性的、而非推断性的变量选择方法。作者在方法部分明确写道：“To assess variable importance, we compute the posterior predictive distribution of the outcome when a single exposure is fixed at a specific quantile...”，这比“识别”的措辞要弱。因此，结论中“identify the relative importance”的说法比其实际证明的方法要宽泛。

四、开放问题¶

正式的贝叶斯变量选择：本文未包含 spike-and-slab 先验。能否将 Antonelli et al. (2020) 的变量选择方法整合到 NB-BKMR 中，以实现对计数数据暴露混合物中关键变量的正式推断？（扎根于：本文方法部分“Variable importance”一节，以及引言中对 Antonelli 2020 的引用。）
因果推断：本文模型是关联性的。如何将 NB-BKMR 嵌入到因果推断框架（如 g-computation、IPW 或 DML）中，以估计社会脆弱性对 COVID-19 死亡率的因果效应？这需要处理未测量混杂，并明确定义因果 estimand。（扎根于：引言中未引用任何因果推断文献，这是一个明显的 gap。）
计算可扩展性：Gibbs 采样器在处理大量县或长时间序列时，计算成本如何？能否引入变分推断或更高效的 MCMC 算法（如 HMC）来加速？（扎根于：本文计算部分仅描述了 Gibbs 采样器，未讨论大规模数据下的计算瓶颈。）
核函数选择：模型对高斯核带宽参数的敏感性如何？是否存在数据驱动的、更鲁棒的核函数选择方法？（扎根于：模型设定部分对核函数 \( K \) 的描述，以及模拟研究中可能未充分探讨的敏感性分析。）

Maintained by 陈星宇 · Homepage · Source on GitHub