R基于贝叶斯加法回归树BART、MCMC的DLNM分布滞后非线性模型分析母婴PM2.5暴露与出生体重数据及GAM模型对比、关键窗口识别

全文链接：https://tecdat.cn/?p=38667

摘要：在母婴暴露于空气污染对儿童健康影响的研究中，常需对孕期暴露情况与健康结果进行回归分析。分布滞后非线性模型（DLNM）是一种常用于估计暴露 - 时间 - 响应函数的统计方法，尤其在假定暴露效应是非线性的情况下。以往基于DLNM的实现多采用双变量基扩展来参数化暴露 - 时间 - 响应曲面，但像样条这样的基函数假设整个暴露 - 时间 - 响应曲面是平滑的，这在实际中可能不符合实际情况，比如暴露仅在特定时间窗口与结果相关时。本文提出了一种基于贝叶斯加法回归树来估计DLNM的框架，通过一组回归树进行操作，各树在暴露 - 时间空间中假设分段常数关系。模拟实验表明，在暴露 - 时间曲面不平滑时，本文模型优于基于样条的模型，而在真实曲面平滑时两者表现相近，且所提方法方差更低，能更精确地识别暴露与未来健康结果相关的关键窗口。最后将该方法应用于估计母亲暴露于PM2.5与出生体重之间的关联。

一、引言

在很多应用场景中，人们常常关注将某个结果对之前一段时间窗口内观察到的暴露情况进行回归分析。这在环境流行病学应用中尤为常见，比如将某天的健康结果对当天及此前几天观察到的暴露（如温度或空气污染）进行回归，或者将出生或儿童健康结果对孕期每日或每周观察到的暴露情况进行回归。在母婴暴露于空气污染这一本文所关注的背景下，通常有两个推断目标：一是估计易感性的关键窗口，也就是暴露能够改变未来表型的时间段；二是估计暴露 - 时间 - 响应函数。近期诸多研究已经确定了母婴暴露于空气污染的关键窗口以及与多种结果之间的关联，包括早产、肥胖、哮喘、喘息、神经发育等情况。而且发现这种线性或非线性关联在孕期的不同周数会有所变化。

二、数据

我们分析来自出生记录以及生命统计数据。数据涵盖了2007年至2015年间（含）的活产、单胎、足月（孕周≥37周）且无已知出生缺陷的新生儿，估计了受孕日期。并且进一步将分析限定在海拔低于6000英尺的人口普查区，这样做既减少了海拔带来的潜在混杂因素影响，也降低了山地地形对暴露数据的影响。

本文关注的主要结果是胎龄别出生体重Z评分（BWGAZ），通过出生图表来获取BWGAZ值，它衡量的是出生体重相对于相同胎龄和性别的儿童预期出生体重的标准差数量。数据还包含个体层面的协变量信息，如母亲的年龄、体重、身高、收入、教育程度、婚姻状况、产前护理习惯、孕期前后是否吸烟，以及种族和是否为西班牙裔等信息。我们将分析限定在协变量信息完整的观测数据上，最终得到300463例新生儿数据。

select(bwgaz, ChildSex, MomAge, GestAge, MomPriorBMI, Race,
                    Hispanic, MomEdu, SmkAny, Marital, Income,
                    EstDateConcept, EstMonthConcept, EstYearConcept)

我们使用来自环境保护署的PM2.5暴露数据，该数据是利用降尺度数据融合的空气质量曲面数据。

sbd\_exp <- list(PM25 = sbd\_dlmtree %>% select(starts\_with("pm25\_")),
                TEMP = sbd\_dlmtree %>% select(starts\_with("temp_")),
                SO2 = sbd\_dlmtree %>% select(starts\_with("so2_")),
                CO = sbd\_dlmtree %>% select(starts\_with("co_")),
                NO2 = sbd\_dlmtree %>% select(starts\_with("no2_")))
sbd\_exp <- sbd\_exp %>% lapply(as.matrix)

三、方法

（一）DLNM框架

在介绍我们提出的方法之前，先简要回顾一下DLNM框架及标准方法。设 y_i 为样本中个体 i（i = 1,..., n）的连续型结果变量，设 x_i = [x_{i1},..., x_{iT}]^T 表示在等间隔时间 t = 1,..., T 观察到的暴露变量向量。在我们的研究中，y_i 代表BWGAZ，而 x_{it} 表示第 i 位母亲在孕期第 t 周暴露于PM2.5的情况。我们还控制一个协变量向量，记为 z_i。高斯DLNM模型如下：

y\_i = f(x\_i) + z\_i^Tγ + ε\_i

其中 f(x_i) 是分布滞后函数，γ 是回归系数向量。分布滞后函数 f(x_i) 可以有多种线性及非线性形式。DLNM允许每个时间点的暴露与结果之间存在独特的非线性关联。一般来说，分布滞后函数定义为：

f(x\_i) = ∑\_{t = 1}^T w(x_{it}, t)

其中 w(x, t) 是将孕期第 t 周的暴露与结果相关联的暴露 - 响应函数。现有的DLNM框架利用交叉基，其中 w 在暴露浓度和时间维度上表示为双变量基扩展。惩罚样条实现允许对暴露 - 时间 - 响应的结构做出一系列假设，比如不同的岭惩罚针对特定时间的收缩，不同的差分惩罚控制曲线的平滑度。基扩展方法（如样条）可对模型进行正则化，以在预测变量存在多重共线性的情况下提高估计效应的稳定性，但这些方法也对DLNM强加了平滑性假设。

（二）树形DLNM方法

我们引入一个基于贝叶斯加法回归树（BART）框架的树和模型来估计暴露 - 时间 - 响应函数 f(x_i)。一般做法是构建二分树，在暴露浓度和时间维度上对随时间变化的暴露变量 x_i 进行划分。展示了基于暴露和时间值定义的二元规则，这些规则将暴露 - 时间空间划分为五个终端节点，记为 η_1,..., η_5。图1b展示了暴露 - 时间空间被划分为五个区域，每个区域对应一个终端节点。一棵树及其相应参数定义了一个分段常数暴露 - 响应函数：

w\_a(x\_{it}, t) = μ_{ab} ，当 x_{it} ∈ η_{ab}

树 T 的分布滞后函数形式与式（2）类似，定义如下：

f\_T(x\_i) = ∑_{t = 1}^T w\_T(x\_{it}, t)

在我们的树形DLNM（树DLNM）框架中，考虑由 A 棵回归树组成的集成。对于树 T_a（a ∈ {1,..., A}），将其 B_a 个终端节点记为 η_{a1},..., η_{aB_a}。每个终端节点 η_{ab} 在时间和暴露浓度上都有相应的限制，由树的分裂规则定义，并且有相应的参数 μ_{ab}。总体而言，每棵树的终端节点定义了暴露 - 时间 - 空间的一个划分，能够灵活地估计暴露 - 时间曲面。与式（3）一样，我们定义树 T_a 中每个暴露 - 时间组合的效应为 w_a(x_{it}, t) = μ_{ab}（如果 x_{it} ∈ η_{ab}）。集成中的每棵回归树都提供了对分布滞后非线性函数 f 的部分估计。形式上，树DLNM的暴露 - 时间 - 响应函数为：

f(x\_i) = ∑\_{a = 1}^A g(x\_i, T\_a)

其中 g(x_i, T_a) 表示树 a 给出的部分估计。树DLNM摒弃了基函数强加的平滑性假设，然而，当不同树的时间和暴露断点交错时，树的集成可以近似平滑函数。模型正则化是树先验的结果，它倾向于只有少量分裂的树，较小的树能确保模型在存在时间相关性时保持稳定，因为每个终端节点会对多个时间点进行平均。

（三）马尔可夫链蒙特卡洛（MCMC）采样器

我们使用MCMC来估计树DLNM。用于BART的MCMC方法不适用于当前模型，原因有两个：一是奇普曼等人的算法依赖于这样一个事实，即任何特定的预测变量向量 x_i 都包含在每棵树的单个终端节点中，而树DLNM将与每个观测值相关的暴露划分到各个终端节点中；二是我们对算法进行了修改，以允许对混杂变量 z 进行参数控制。由于这些差异，我们为树DLNM模型提出了一种替代的MCMC方法，具体来说，我们使用标准分析技术对 γ 进行积分，然后应用贝叶斯反向拟合同时基于每棵树 T_a 定义的划分来估计部分暴露 - 时间 - 响应的效应，MCMC采样器的详细信息可在补充材料中找到。

 

模拟研究

模拟设置

我们开展模拟研究，旨在对比树DLNM和树DLNMse（树DLNM的一种变体）与已有的DLNM方法（如使用惩罚和非惩罚样条的那些方法）的实际表现。模拟的关键在于对比代表平滑和非平滑暴露 - 时间 - 响应函数的模拟场景下各模型的表现。
我们依据相关公式模拟数据，样本量等同于我们的暴露数据（n = 300463）。为准确体现空气污染暴露数据中存在的自相关性，我们使用了数据分析中的PM2.5暴露数据。用每个观测值的连续37周来模拟足月妊娠情况，以此模拟DLNM曲面。我们考虑了四种模拟的暴露 - 时间 - 响应函数，每种对应不同的真实模型（树DLNM、树DLNMse、基于样条的平滑DLNM以及线性DLM）。这四种DLNM场景分别如下：

• A场景：在第11 - 15周暴露呈现分段常数效应。

• B场景：在第11 - 15周暴露呈现线性效应。

• C场景：在第11 - 15周暴露呈现平滑、非线性效应（逻辑斯蒂形状）。

• D场景：在第11 - 15周暴露呈现平滑、非线性效应（逻辑斯蒂形状），且时间上的平滑效应在第13周达到峰值，向两边各延伸大约五周。
  我们使用对数转换后的暴露数据生成结果，所有场景都以对数暴露值1为中心。暴露 - 时间曲面的几个横截面展示在图2和图3中，每个场景下DLNM曲面的代数细节以及图形表示可在补充材料C.1部分找到。
  我们生成了一组协变量（五个标准正态分布变量、五个概率为0.5的二项分布变量）以及来自标准正态分布的相应系数。通过使用臭氧数据引入季节性趋势，具体而言，每5周（第5、10、…、35周）添加一个随机臭氧效应，每次的臭氧测量值以均值为零进行中心化，缩放至标准差为1，然后乘以从均值为0、方差为0.04的正态分布中抽取的值。这样，每个模拟数据集都有不同的季节性趋势，且与暴露（PM2.5）和结果都相关。我们设置误差方差，使得Var[f(xi)]/σ2 = 1/1000，以此代表现实的信噪比，并且在每个场景中运行500次模拟重复。

模拟估计量与对比

树DLNM和树DLNMse采用了描述的先验设置。在所有对数暴露值的0.01百分位数到99.9百分位数之间，选取30个等间距的值作为暴露维度上潜在的分割点。经过5000次迭代的预热期后，每个模型运行15000次迭代，每隔10次抽取一次结果。
我们将树DLNM和树DLNMse与几种基于样条的惩罚和非惩罚DLNM模型进行对比，这些模型介绍如下：

• GAM：基础模型，在暴露和时间维度上由秩为10的惩罚三次B样条平滑器定义，带有二阶惩罚，采用REML（限制最大似然估计）进行估计。

• DLM：在暴露浓度上采用线性假设的GAM模型。

• GLM - AIC：在暴露和时间维度上选择最优数量（自由度从1到10）的非惩罚二次B样条，通过最小化AIC（赤池信息准则）来确定。

• GAMcr：通过用三次回归样条替换GAM中的三次B样条基，并对二阶导数添加惩罚来定义。

• GAM - exp：GAM模型，用可变岭惩罚替换二阶惩罚。
  为评估模型性能，我们将每个模型的DLNM以对数暴露值1为中心，并在一系列点组成的网格上评估估计的DLNM。在每个模型中，我们纳入所有10个模拟协变量以及年份和月份的指示变量，以此控制额外的季节性趋势。对暴露浓度值进行对数转换，以减少暴露数据的偏态，并使样条基模型中的节点等间距分布。对于树DLNM而言，对响应进行对数转换与否不会影响其结果，但对树DLNMse的平滑性有影响。

模拟结果

模型性能的汇总指标展示在表中。我们通过整个暴露 - 时间曲面的均方根误差（RMSE）以及分解到模拟关键窗口内外的RMSE来对比各模型。还展示了95%置信区间的经验覆盖率以及平均置信区间宽度。此外，对比了在模拟关键窗口内网格点上识别非零效应的概率（TP）、在模拟关键窗口外网格点上错误地识别出非零效应的概率（FP），以及正确识别非零效应的精度（TP/(TP + FP)）。我们将真实暴露 - 时间曲面中的非零效应定义为在 - 0.005到0.005区间之外的任何效应，以此考虑B、C和D场景的情况（B场景在第11 - 15周各处都有非零效应，D场景各处都有非零效应）。

为评估模型性能，我们将每个模型的 DLNM 以对数暴露值 1 为中心，并在一系列点组成的网格上评估估计的 DLNM。在每个模型中，我们纳入所有 10 个模拟协变量以及年份和月份的指示变量，以此控制额外的季节性趋势。对暴露浓度值进行对数转换，以减少暴露数据的偏态，并使样条基模型中的节点等间距分布。对于 树DLNM 而言，对响应进行对数转换与否不会影响其结果，但对 树DLNMse 的平滑性有影响。

模拟结果

模型性能的汇总指标展示在表 1 中，各指标从不同角度反映了不同模型在不同场景下的表现情况，以下是详细介绍：

整体均方根误差（RMSE）体现了模型估计值与真实值之间的偏差程度，从整体来看，在 A、C 和 D 场景中，树DLNM 和 树DLNMse 总体 RMSE 与其他模型相当甚至更好，如在场景 A 中，树DLNM 的整体 RMSE 为 0.086，相对一些基于样条的模型（如 GAM 的 0.294、GLM - AIC 的 1.531 等）表现出更好的估计准确性。

无效应区域 RMSE 则反映了在暴露 - 时间 - 响应曲面中效应为零区域各模型的偏差情况，在所有场景下，基于树的方法（树DLNM 和 树DLNMse）在这方面具有优势，比如场景 D 中，树DLNM 的无效应区域 RMSE 为 0.041，相对部分其他模型更低，意味着在无效应区域它们的估计更接近真实情况，这得益于树特定参数上的收缩先验降低了方差，从而在无效应区域导致更低的 RMSE。

效应区域 RMSE 展示了在存在效应的区域各模型的表现，像在场景 A 中，树DLNM 在效应区域的 RMSE 为 0.213，在该场景下比部分基于样条的模型更优，不过在场景 B 中，树DLNM 和 树DLNMse 在这方面的 RMSE 高于部分基于样条的模型，因为基于样条的模型在数据点稀少的极端暴露值处插值效果更好。

总体置信区间宽度指标能看出模型估计的精确程度，我们的模型（树DLNM 和 树DLNMse）平均置信区间宽度相对较小，尤其在时间边界或极端暴露浓度处优势明显，这在一定程度上是因为基于树的模型估计中缺乏 “波动”，例如在场景 C 中，树DLNM 的总体 CI 宽度为 0.94，相较于部分其他模型更窄，有助于缩窄置信区间以及降低 RMSE，尤其在零效应区域。

真阳性率（TP）表示在模拟关键窗口内网格点上正确识别非零效应的概率，误报率（FP）则是在模拟关键窗口外网格点上错误地识别出非零效应的概率，精度（TP/(TP + FP)）综合体现了模型正确识别非零效应的能力。树DLNM 和 树DLNMse 在所有模拟场景中的精度都是较高的，像在场景 A 中精度能达到 1.00（树DLNM）和 0.98（树DLNMse），高精度源于误报率（FP）近乎为零，但在 B 和 D 场景中存在真阳性率（TP）较低的权衡情况。

图2和图3展示了使用树DLNM、树DLNMse和GAMcr模型估计得到的暴露 - 时间 - 响应曲面的横截面。图中显示的非零估计值意味着对于具有特定时间和暴露 - 浓度值的任何观测值，响应发生了变化。在A、C和D场景中，树DLNM和树DLNMse总体RMSE与其他模型相当甚至更好。在所有场景下，基于树的方法在暴露 - 时间 - 响应曲面的零效应区域具有最低的RMSE。图2突出显示了树DLNM和树DLNMse能够清晰区分有无效应的时间段的能力，树特定参数上的收缩先验降低了方差，从而在无效应区域导致更低的RMSE。在非零效应区域，我们的模型在A场景下比基于样条的模型具有更低的RMSE，在C和D场景下则与之相当。在B场景中，树DLNM和树DLNMse在非零效应区域的RMSE更高，因为基于样条的模型在数据点稀少的极端暴露值处插值效果更好。

plot(tDLNM, plot.type = "slice",

图对比了基于树的模型如何在暴露值边界处减弱效应，而GAMcr则线性延续趋势。除B场景外，基于树的模型的覆盖率接近名义覆盖率。在B场景中，所有模型的覆盖率都低于名义覆盖率，但树DLNM和树DLNMse表现最佳，各自的曲面覆盖率为87%。此外，我们的模型平均置信区间宽度最小，这在时间边界或极端暴露浓度处尤为显著，因为基于样条的模型在这些地方的“波动”更明显（如图2和图3所示）。基于树的模型估计中缺乏“波动”，这有助于缩窄置信区间以及降低RMSE，尤其在零效应区域。而且，树DLNM和树DLNMse在模拟重复间的差异要小得多。
B场景虽然对DLM来说看似自然，但存在几个难点。首先，树DLNM要进行恰当估计，需要在正确的关键窗口内有很多跨越暴露浓度的分割树。其次，当数据稀疏时（如本场景中的高、低浓度情况），树DLNM会减弱效应。第三，在高浓度时，存在从零到较大效应的跳跃，平滑方法无法适应这种情况，特别是DLM由于平滑性假设会将关键窗口延伸到远超真实效应的时间段。树DLNM和树DLNMse在所有模拟场景中的精度都是最高的（见表1），高精度源于误报率（FP）近乎为零，但在B和D场景中存在真阳性率（TP）较低的权衡。图中的横截面图展示了树DLNM和树DLNMse适应非平滑暴露 - 时间响应曲面的能力，补充图显示了每周至少在一个暴露值上检测到非零效应的概率，结果表明基于样条的方法在真实关键窗口外误分类周数的概率要高得多，而基于树的模型能适应暴露 - 时间 - 响应曲面的变化平滑性，很少在真实关键窗口外检测到非零效应，关键在于树DLNM和树DLNMse所检测到的关键窗口大概率是正确的。

数据分析

我们使用树DLNM和树DLNMse来估计孕期前37周母亲暴露于PM2.5与儿童BWGAZ之间的关系。通过使用每周暴露数据，我们将任何方法能够识别关键窗口的时间分辨率限制在周的层面。为进行对比，我们还应用了使用惩罚三次回归样条的DLNM（GAMcr）以及DLM方法，并且控制了母亲的基线特征以及季节和长期趋势。母亲的特征包括：孕前年龄（二次拟合）、体重、是否吸烟（孕期前后）、收入、教育程度、首次接受产前护理的时间、种族、是否为西班牙裔、海拔以及居住的郡县。我们不控制胎儿性别或孕周，因为结果变量BWGAZ已经针对这些因素进行了调整。此外，我们使用受孕年份和月份的指示变量来调整季节性效应。
对于树DLNM和树DLNMse，我们使用与模拟中相同的超参数，让模型运行5000次迭代的预热期，接着运行15000次迭代，并保留马尔可夫链蒙特卡洛（MCMC）采样器中每10次抽取的结果。我们在暴露维度上指定30个等间距的潜在分割点，范围从对数暴露值的0.1百分位数到99.9百分位数，尝试过不同数量的潜在分割点，但结果并无差异。在树DLNMse中，我们将平滑参数σx设置为对数暴露值标准差的一半。GAMcr和DLM模型使用与模拟中相同的设置，所有模型的DLNM估计都以暴露值中位数（约7 µg/m3 ）为中心，将包含暴露 - 时间 - 响应中95%置信区间不包含零的区域的任何一周定义为关键窗口。

DLNM结果

树DLNMse的后验平均暴露 - 时间 - 响应估计值展示在图4a中，低于中位数的PM2.5暴露与BWGAZ的增加相关，高于中位数的暴露浓度则表明BWGAZ略有下降，但95%可信区间并不能让我们认为这种下降与零有显著差异，这种模式在整个孕期各周都存在。

图4b展示了在第5、15、25和35周的暴露 - 时间 - 响应曲面的横截面，表明存在一个跨越整个孕期的关键窗口。基于树DLNMse，孕期内从中位数（7.0 µg/m3 ）的PM2.5暴露变化到第25百分位数（5.89 µg/m3 ），将导致BWGAZ的累积平均增加量为0.0132（95%置信区间[0.0003, 0.0354]），换算成实际出生体重大约增加5.74g（95%置信区间[0.11, 15.41]）（这是近似值，因为BWGAZ考虑了孕周和胎儿性别因素）。非线性关联表明，PM2.5暴露进一步降低到第10百分位数（5.02 µg/m3 ）时，BWGAZ平均增加量为0.055（95%置信区间[0.016, 0.090]），换算成实际出生体重大约增加24.1g（95%置信区间[7.155, 39.10]），这些结果意味着将PM2.5降低到当前国家环境空气质量标准以下，会使该人群的平均出生体重增加。

GAMcr的平均暴露 - 时间 - 响应估计值展示在图4a中，与树DLNMse的估计值很相似，正如在模拟中看到的那样，二者在尾部表现上存在差异，GAMcr随着效应延续趋势且区间较大，尽管在低暴露水平下GAMcr的点估计值较大，但较大的置信区间包含了零。相反，树DLNMse逐渐减小并估计出一个较小的效应，其区间明显更小且不包含零，树DLNMse在边界附近区间较小是因为这些边界区域被划分到包含内部区域的终端节点中，因此得到相同的估计值。我们发现PM2.5增加与BWGAZ降低之间存在关联，这与以往文献一致，孙等人2016年的一项荟萃分析发现，孕期PM2.5增加10 µg/m3与出生体重降低15.9g相关（95%置信区间[ - 26.8, - 5]），孕中期和孕晚期暴露增加也被确定与出生体重存在非零负相关，朱等人2016年在另一项荟萃分析中也报告了类似结果，斯特里克兰等人2019年发现，在整个孕期，对于出生体重分布的较高百分位数，PM2.5与出生体重之间的关联程度会增加，最后，一项研究调查PM2.5的单个化学成分发现，孕期每个阶段母亲暴露于PM2.5都存在非零的低出生体重风险增加情况。

对比灵活性稍差的模型替代方案

为进行对比，我们拟合了树DLNM、DLM以及几种线性模型来对比结果，这些模型的结果都与树DLNMse的结果相符，关于这些方法的更多细节可在补充材料D部分找到。

讨论

在这项工作中，我们提出了一种基于树的DLNM方法，用于估计随时间变化的一系列污染暴露与连续出生结果之间的关联。树DLNM消除了暴露 - 时间响应曲面中的平滑性假设，树DLNMse仅在暴露 - 浓度维度施加平滑性而不在时间维度施加，树DLNM还有可能在暴露 - 响应函数内考虑测量误差。通过放松时间维度的平滑性假设，我们的新方法能更精确地识别易感性关键窗口，有着重要的应用价值和理论意义。

参考文献

Berrocal, V. J., Gelfand, A. E., Holland, D. M., & Statisti-Cian, S. (2010). A Spatio-Temporal Downscaler for Output From Numerical Models. Journal of Agricultural, Biological, and Environmental Statistics, 15(2), 176–197.

Bose, S., Rosa, M. J., Mathilda Chiu, Y.-H., Leon Hsu, H.-H., Di, Q., Lee, A., Kloog, I., Wilson, A., Schwartz, J., Wright, R. O., Morgan, W. J., Coull, B. A., & Wright, R. J. (2018). Prenatal nitrate air pollution exposure and reduced child lung function: Timing and fetal sex effects. Environmental Research, 167, 591–597.

Carvalho, C. M., Polson, N. G., & Scott, J. G. (2010). The horseshoe estimator for sparse signals. Biometrika, 97(2), 465–480.

Chang, H. H., Reich, B. J., & Miranda, M. L. (2012). Time-to-Event Analysis of Fine Particle Air Pollution and Preterm Birth: Results From North Carolina, 2001-2005. American Journal of Epidemiology, 175(2), 91–98.

Chang, H. H., Warren, J. L., Darrow, L. A., Reich, B. J., Waller, L. A., & Chang, H. H. (2015). Assessment of critical exposure and outcome windows in time-to-event analysis with application to air pollution and preterm birth study. Biostatistics, 16(3), 509–521.

Chipman, H. A., George, E. I., & McCulloch, R. E. (1998). Bayesian CART model search. Journal of the American Statistical Association, 93(443), 935–948.