荷兰地区层面温度相关超额死亡率的分析

《Environmental Epidemiology》：District-level analysis of temperature-related excess mortality in the Netherlands

　　背景：虽然大多数流行病学研究都是在较大的地理尺度上评估温度与死亡率之间的关联，但最近的研究指出，由于暴露程度、脆弱性和适应能力的局部差异，温度相关的死亡风险在空间上是异质性的。本研究的目的是评估环境空气温度对荷兰各区全因死亡率的影响。

方法：采用三阶段设计来评估2010-2019年期间3116个区以及两个年龄组（0-74岁、75岁以上）的日平均气温与死亡率之间的关联。首先，使用具有分布滞后非线性模型的案例时间序列设计来估计355个荷兰市镇的温度-死亡率关联。其次，利用市镇级别的 demographic、经济社会、气候和地形预测因子，通过多变量元回归汇总各市镇的估计结果。最后，使用元回归模型将风险估计值降尺度到区级，并估算标准化的热相关和冷相关超额死亡率。

结果：熱相关的超额死亡率在城市地区较高，而与寒冷相关的超额死亡率在荷兰北部和南部地区更为明显。在研究期间，热与每年平均346例超额死亡（经验95%置信区间：44-558）相关，而冷与每年11,988例超额死亡（7,554-15,609）相关。这相当于每10万人口年标准化超额死亡率分别为2.2（0.3-3.5）和73.3（46.8-94.7）。75岁以上组的死亡风险比年轻组更为显著。

结论：我们的结果表明，与温度相关的标准化死亡率存在地理异质性，城市地区的热相关死亡率较高，而寒冷相关的死亡率则表现出明显的区域聚集现象。这些发现有助于地方政府根据特定地区的脆弱性实施气候适应措施，以预防与热和寒冷相关的超额死亡。

本研究的意义：通过描述小区域内的温度相关死亡负担，本研究可能为荷兰的公共卫生政策提供参考，并支持未来对气候变化下温度相关健康风险的评估。此外，尽管我们的分析主要集中在荷兰，但其优势在于使用了高空间分辨率的暴露和健康数据，以及应用了新的建模技术，从而能够估算小区域的死亡风险。

引言：大量证据表明，短期暴露于次优环境温度与死亡风险增加有关。全球范围内，次优温度每年导致9.4%的死亡案例，其中0.9%可归因于热，8.5%归因于寒冷。由于气候变化，预计在未来几十年内极端天气事件（如热浪）的频率和强度都会增加。目前的预测表明，在温暖气候下，与热相关的死亡率将会增加。同时，虽然暴露于寒冷天气的情况可能会减少，但由于人口老龄化，脆弱人群预计会增加。因此，适应和准备应对极端温度对于减少全球范围内的热相关和冷相关健康负担仍然至关重要。在过去的几十年中，流行病学研究主要关注在较大地理范围内估计温度-死亡率风险，这些研究通过汇总暴露数据和健康结果来获得足够的统计能力。然而，越来越多的文献表明，温度的健康影响在空间上是异质性的，因为温度暴露、人口脆弱性和气候适应能力会因地区特定的生物生理、社会文化、基础设施和技术因素而有所不同。因此，小区域分析更适合识别对温度的局部脆弱性并计算空间详细的健康影响。随着收集的高精度地理分辨率数据的增加，发展出了新的小区域流行病学分析方法。这些方法已在英格兰和威尔士的社区层面、挪威的市镇层面以及法国的伪县层面成功应用。在荷兰进行的研究主要依赖于全国汇总的时间序列数据来评估温度-死亡率关联。因此，温度暴露和人口脆弱性的空间差异被忽略了，而这可能是导致全国范围内温度-死亡率关联异质性的原因之一。本研究的目的是使用Gasparrini等人提出的三阶段设计，在荷兰的小区域层面上评估环境空气温度对全因死亡率的影响，以识别小规模的风险模式。首先，我们估计了所有荷兰市镇的日平均气温与全因死亡率之间的关联。接着，我们使用市镇级别的demographic、经济社会、气候和地形预测因子，通过多变量元回归汇总了各市镇的具体风险估计值。最后，我们应用元回归模型将风险估计值降尺度到区级，并估算了标准化的热相关和冷相关超额死亡率。通过这种具有两级分析（从市镇到区级）的新研究设计，我们能够在不对每个区分别建立模型的情况下，以尽可能高的分辨率分析数据并估算健康影响。

研究区域和时期：荷兰拥有受北海影响的温带海洋性气候，地形相对平坦。2019年，该国人口约为1730万，面积为大约41,543平方公里。该国共有355个市镇，平均每个市镇由九个区组成。2019年，有3,116个有人居住的区。平均而言，每个区的面积为11.18平方公里（中位数=6.49平方公里），人口为5,544人。选择2010-2019年作为研究时期，是为了避免2019年冠状病毒病（COVID-19）大流行的潜在混杂效应，同时之前的流行病学研究已发现2010年荷兰首次启动国家热浪计划前后温度相关死亡风险的差异。

死亡率和暴露数据：我们从荷兰统计局收集了2010年1月1日至2019年12月31日的全因死亡病例数据。基于这些个体死亡病例，为荷兰355个市镇内的3,116个区创建了特定的全因死亡率日时间序列。日温度数据来自Bussalleu等人开发和验证的高分辨率时空模型。使用两阶段混合建模方法，估算了2003-2020年期间欧洲的每日平均、最低和最高环境温度。简而言之，地面以上2米处的日环境温度以1×1公里的精细空间分辨率进行预测，结合了卫星数据、基于站点的温度测量、大气再分析数据和土地利用信息。对于当前分析，我们通过将网格化温度数据与2019年的行政区边界叠加，得到了2010-2019年期间的特定区日平均温度值。对于每个区，我们计算了与其边界相交的所有网格单元的日平均温度的面积加权平均值，权重与网格单元与区的交叉面积成比例。此外，为了进行敏感性分析，我们还从ERA5再分析数据集中获取了日温度数据。ERA5数据集提供了基于28×28公里尺度的历史气候和天气数据。地表温度参数对应于陆地、海洋或内陆水域上方2米处的空气温度。使用相同的面积加权平均方法，根据与每个区相交的网格单元，得出了特定区的日平均温度值。

元预测因子：我们编制了一个包含可能与温度脆弱性相关的区级和市镇级特征的数据集。我们识别了可能解释极端温度脆弱性地理差异的因素，并利用区和市镇级别的数据将其具体化。所得数据集包括21个关于人口统计、建筑、土地利用和天气特征的变量。人口统计和建筑特征的数据在市镇和区级别获得。土地利用和天气特征的数据使用了2019年的行政边界进行评估。有关元预测因子的更多信息及其定义和来源，请参见表S1（https://links.lww.com/EE/A488）。21个元预测因子的空间分布显示出高度碎片化的模式（图S1）。变量之间的相关性在市镇和区级别上相似（图S2）。为了降低元预测因子数据集的维度并确保这些变量的预测能力，我们进行了主成分分析（PCA）。这种方法产生了一组新的无关变量，称为主成分（PCs）。通过PCA减少指标的数量，我们能够识别出由多个小区域特征组合形成的脆弱性模式，避免了处理长列表变量时可能出现的任意选择或计算挑战。

统计分析：我们采用三阶段设计来评估短期温度暴露与区级全因死亡率之间的关联。在第一阶段，我们应用案例时间序列设计来估计市镇级别的温度-死亡率关联。案例时间序列设计允许应用为时间序列分析开发的灵活建模技术，而无需将数据汇总到单个时间序列中。因此，对于355个市镇中的每一个，我们使用条件伪泊松回归模型建模了特定区的温度和死亡率序列。我们对两个年龄组（0-74岁和75岁以上）分别进行了建模。温度-死亡率关联通过分布滞后非线性模型（DLNMs）进行建模。DLNMs是一种广泛用于建模即时和延迟非线性暴露-响应关联的框架。由于日温度与死亡率之间的关系表现出滞后和非线性模式，DLNMs在文献中被广泛用于估算温度-死亡率关联。在DLNMs中，风险关联通过跨基项进行建模，该基项定义了暴露-响应和滞后-响应关系。我们使用三次B样条函数对暴露-响应进行建模，该函数在日市镇温度分布的第10、75和90百分位数处设置了三个结点。滞后-响应使用自然三次样条函数进行建模，其内部结点在滞后对数尺度上均匀分布。滞后基的设定没有包括截距。我们定义了0-21天的滞后期，以捕捉温度对死亡的延迟和间接影响。这种暴露-滞后响应关系的参数化方法来自之前的温度-死亡率研究。模型中还包括了一个考虑基线风险和时间趋势变化的项。这项调整缩小了温度分布的第51百分位数和第99百分位数之间的范围。然后，我们计算了各地区的特定死亡率相对风险（RR）和95%置信区间（CI），这些数据是基于第99百分位数和第1百分位数的温度与MMT（最小平均温度）进行比较的，以获得与热和冷相关的死亡率RR。我们还使用之前描述的方法，估算了年龄特定的年度超额死亡率，包括总死亡人数以及由于温度高于（与热相关的超额死亡率）和低于（与冷相关的超额死亡率）MMT而导致的死亡人数比例。这些估计值被用来根据2019年荷兰统计局的人口数据，计算每10万人中每年因热和冷导致超额死亡率的具体比率。随后，使用年龄特定的标准权重对这些超额死亡率进行了标准化。相关的不确定性通过蒙特卡洛模拟（n = 500次）量化为95%的经验置信区间。我们测试了四个地区级变量对热和冷相关死亡率的影响：平均收入、树木覆盖率、燃气使用量和地址密度。对于每个变量，地区被分为两组——低于和高于中位数。然后，我们检查了两组之间的总标准化死亡率是否是否存在显著差异。

我们进行了几次敏感性分析。首先，我们修改了第一阶段交叉基础项中的暴露-反应函数的参数化：在一个模型中，结节位置被移动到了第20、70和85百分位数；而在另一个模型中，二次B样条被替换为自然样条。其次，我们将滞后周期缩短到3天，以评估温度的更直接影响。最后，我们使用ERA5再分析数据集的温度数据重复了分析。

所有统计分析都是在R（版本4.5.0）中使用的，使用了dlnm、gnm和mixmeta包。

最终分析包括了2010年至2019年期间所有原因导致的1,444,448例死亡。2019年平均每个地区每10万人中有869例死亡（图1）。平均日平均温度显示出适度的地理梯度，从东北部的9.9°C到南部和西部的11.5°C（图1）。最低日平均温度和最高日平均温度也存在地理梯度，沿海地区的最低温度较高，而内陆地区的最高温度较低（图S3，https://links.lww.com/EE/A488）。使用ERA5再分析数据估算的各地区的平均、最低和最高日平均温度在补充材料中报告（图S4，https://links.lww.com/EE/A488）。从Bussalleu等人26和ERA5数据集估算的各地区日温度高度相关。各地区的平均皮尔逊相关系数（标准差）为0.992（0.01）。

图2显示了按年龄组汇总的总体累积暴露-反应关系的估计值，在元预测值的平均值处进行预测。0-74岁人群的汇总曲线显示出温度对死亡率的中等到无影响。特别是在极端温度下，曲线没有显示死亡风险的增加。75岁以上年龄组的曲线从第5百分位数开始呈现出典型的V形。然而，从最低温度到第5百分位数，死亡率随着温度的降低而减少。两个年龄组的MMT都相对较高，分别对应于第99百分位数和第94百分位数（表S2，https://links.lww.com/EE/A488）。75岁以下人群的与热相关的RR [95% CI] 分别为1.00 [1.00–1.00] 和 1.08 [0.99–1.18]（表S2，https://links.lww.com/EE/A488）。75岁以上人群的RR分别为1.05 [1.01–1.11] 和 1.18 [1.11–1.25]（表S2，https://links.lww.com/EE/A488）。

在第二阶段模型中，我们包括了五个主要特征（PCs）作为元回归变量，这些变量解释了78%的方差（图S5，https://links.lww.com/EE/A488）。五个PC的地理分布以及PC与21个变量之间的相关性在图S6和S7中展示（https://links.lww.com/EE/A488）。尽管包括了元预测变量，元回归模型仍然显示出一些残余异质性（Cochran’s Q检验：P = 0.01；I2 = 5.0%），表明部分风险变异无法通过模型解释。似然比检验也表明，年龄和PC1在元回归模型中是统计上显著的（表S3，https://links.lww.com/EE/A488）。

图3显示了各地区和市镇因热和冷导致的年度标准化超额死亡率的地理分布。与热和冷相关的超额死亡率显示出不同的模式，但并非完全相反。地区级地图显示，人口密集和城市化的地区热相关超额死亡率较高。冷相关死亡率的大小大于热相关死亡率；它们显示出中等的空间梯度，北部、西南部和东南部地区的脆弱性较高，而内陆地区的脆弱性较低（图3）。与地区地图相比，市镇地图显示了更为连贯和易于解释的超额死亡率空间模式（图3）。

表S5显示了整个国家每年因热和冷导致的超额死亡人数和每10万人年的标准化超额死亡率。平均而言，在整个研究期间，高于MMT的温度每年导致346例超额死亡（经验95% CI = 44,558），而低于MMT的温度每年导致11,988例超额死亡（7,554–15,609），相当于因热导致的标准化超额死亡率为每10万人年2.2例（经验95% CI = 0.3, 3.5），因冷导致的标准化超额死亡率为每10万人年73.3例（46.8–94.7）。

图4显示了按地区平均收入、树木覆盖率、燃气使用量和地址密度分层的每10万人年的标准化超额死亡率。不同收入水平的地区之间的标准化死亡率差异很小，高收入地区的热相关死亡率较高，而低收入地区的冷相关死亡率较高。尽管地区和市镇级别的地图表明城市化地区的热相关标准化死亡率较高，但在低地址密度和高地址密度的地区之间没有发现标准化死亡率的显著差异。

改变暴露-反应参数化中结点的位置以及使用空间分辨率更粗糙的ERA5温度数据并未改变主要发现（图S10，https://links.lww.com/EE/A488）。通过自然样条建模的暴露-反应曲线在极端寒冷温度下没有显示死亡率的降低。应用3天的滞后周期也影响了暴露-反应曲线的形状，年轻人和老年人的热相关风险增加，而寒冷温度没有相关风险（图S10，https://links.lww.com/EE/A488）。

这项研究提供了2010年至2019年间荷兰3,116个地区温度与死亡率关联的全国性评估。我们的地图显示，城市地区的热相关超额死亡率较高，而北部和南部地区的冷相关超额死亡率较高。老年人（75岁以上）的热和冷相关死亡率风险增加，而年轻人（0-74岁）的关联为零到中等。在轻微寒冷条件下，死亡风险适中，而在极端寒冷条件下，死亡风险较低。我们发现高收入地区的热相关标准化超额死亡率有所增加，而低收入地区的冷相关标准化超额死亡率也有所增加。与地区级的树木覆盖率、燃气使用量和地址密度无关，标准化死亡率没有显著差异。在研究期间，高温（即超过MMT的环境温度）每年导致平均346例死亡，而低温导致11,988例死亡。这相当于每10万人年的热相关标准化超额死亡率为2.2 [0.3–3.5]，冷相关率为73.3 [46.83–94.74]。

地区级的温度相关标准化死亡率地图显示了热相关死亡率的明显空间异质性。这种空间变化可能不仅由特定地点对温度的社会生理脆弱性差异引起，还可能由方法学因素造成。特别是，估算的各地区MMT决定了哪些日子被归类为热天或冷天，这取决于MMT在每个地区温度分布中的位置。MMT较高的地区热天较少，因此，根据定义，这些地区的热相关超额死亡率较低。在市镇层面汇总效应后发现，热相关死亡率的空间模式更加一致，人口密集的市镇死亡率较高。英格兰和威尔士以及法国的城市化地区也发现了高的热相关死亡率；相比之下，挪威的农村地区热相关死亡率较高。我们的地图显示，荷兰北部、西南部和东南部的冷相关死亡率较高。同样，在英格兰和威尔士，偏远地区（北部英格兰和西南部地区）观察到了较高的冷相关死亡率；而在挪威和法国，冷相关死亡率没有明显的地理聚集现象。

我们的结果与之前的研究一致，表明年龄对温度相关死亡率有显著影响。虽然在0-74岁年龄组中发现的效应为中等至无，但75岁以上年龄组在极端温度下显示出高热相关和冷相关死亡率风险，暴露-反应曲线在第五百分位数之后呈现出典型的V形。直到第五百分位数，暴露-反应函数表明死亡风险下降。这种异常可能反映了室外温度与实际人群暴露之间的差异，因为在极端寒冷温度下，人们更倾向于待在室内。我们发现高收入地区的热相关标准化超额死亡率有所增加，而低收入地区的冷相关标准化超额死亡率也有所增加。与其他研究结果类似，其他研究也发现社会经济贫困程度较高和教育水平较低的社区中冷相关死亡率较高。与我们的结果一致，另外两项小范围研究也未发现树木覆盖率、燃气使用量和地址密度对标准化死亡率有显著影响。例如，与农村地区相比，城市地区的热致死亡率有所上升，而在教育水平较低的市政区域内这种趋势更为明显。12,17,18 在我们的分析中，尽管没有观察到与背景因素相关的超额死亡率存在显著差异，但一个可能的解释是，在研究区域内，对温度的脆弱性是由多种脆弱性因素共同作用的结果，而不仅仅是由单一的地区级特征决定的。de Schrijver等人的研究支持了这一观点，12 他们发现，在农村地区，人口老龄化和社会经济因素在调节与热相关的风险中起着更重要的作用，而在城市地区，热致死亡的主要驱动因素是环境因素，如温度和空气污染。另一个解释是，使用地区级变量可能会掩盖个体层面的重要差异。因此，我们的发现并不排除个体层面存在效应修饰的可能性。此外，由于我们仅比较了两个组别（高于和低于中位数的组别）之间的死亡率差异，可能会忽略那些差异更大的组别。

与全球研究结果一致，我们的研究发现寒冷对死亡率的影响明显大于热。1,34 这一结果尤为重要，因为流行病学预测表明，在不同的气候变暖情景下，由于人口老龄化，荷兰的寒冷相关死亡率可能会保持不变甚至增加。43 在解释当前分析中估计的寒冷相关死亡率负担时，必须考虑一些因素。首先，寒冷温度与死亡率之间的关联似乎对滞后期的选择很敏感。具体来说，当应用3天的滞后期时，寒冷相关的相对风险（RR）几乎消失，这表明寒冷相关死亡率可能主要由寒冷暴露的间接和延迟效应驱动，而不是其直接影响。此外，热相关和寒冷相关死亡率的归因在很大程度上取决于参考温度（MMT）在整体温度分布中的位置。在我们的研究中，年轻组和老年组的国家级MMT分别为23.4（第99百分位）和19.7（第94百分位）。因此，在整个研究期间，大多数天数都导致了寒冷相关的超额死亡率，而热效应仅限于温度高于MMT的相对较少的几天。挪威（女性和男性的第93和第94百分位）以及英格兰和威尔士（第95百分位）也观察到相当高的MMT。17 相比之下，早期的荷兰研究使用汇总的温度数据得出的MMT较低。Huynen等人的研究21发现65岁以上年龄组的MMT为16.5°C，Folkerts等人的研究23发现65岁以上年龄组的MMT为15.3°C。与之前的研究相比，观察到的较高MMT可能是由于方法上的不同，例如使用了模型化的地区级温度数据，而不是早期研究中使用的监测站温度测量数据。

我们的发现也与最近使用类似方法评估国家层面小区域死亡风险的欧洲研究一致。17,18 这些研究分别报告了英格兰和威尔士以及挪威的热相关超额死亡率为每10万人年1.6 [1.2–1.9]例和0.8 [0.6–1.1]例，寒冷相关超额死亡率为每10万人年122.3 [112.9–131.5]例和73.4 [69.7–77.1]例。多国和多城市（MMC）研究也得出了类似的死亡率估计。例如，Masselot等人的研究44报告了2000年至2019年间欧洲854个城市地区的热相关标准化死亡率为每10万人年9 [7–10]例，寒冷相关死亡率为每10万人年83 [74–92]例。MMC研究中的热相关标准化率平均高于我们研究中的率。这部分是由于参考温度位置的不同。如上所述，使用汇总数据的研究报告的MMT较低，可能导致热相关死亡率较高。

值得注意的是，使用分辨率较低的温度数据进行敏感性分析得出的估计结果与主要分析的结果相当（图S10，https://links.lww.com/EE/A488）。这一结果特别值得注意，因为它表明较低的分辨率温度数据也可以充分描述与温度相关的风险，而不需要高度详细的温度空间信息。Mistry等人45也报告了类似的发现，他们表明从ERA5再分析数据得出的温度相关风险估计与基于站点观测的数据相当。

本研究的一个关键优势是使用了在小区域级别汇总的高分辨率温度数据。这种方法克服了基于在大地理范围内汇总的温度数据的分析的局限性，后者无法捕捉温度-死亡率关联的细尺度变化。此外，使用成熟的方法论框架使我们能够得出多个效应估计，从而便于将我们的结果与荷兰其他地区以及其他欧洲研究进行直接比较。另一个方法论优势是在元回归模型中使用主成分（PCs）作为元预测因子，而不是单独的脆弱性变量。这种方法能够将相关的局部特征汇总为较少数量的复合预测因子，减少了数据集的维度并最小化了变量之间的共线性。

尽管有这些优势，仍需承认一些局限性。首先，暴露-响应曲线的参数化可能引入了热相关RR估计的不确定性。具体来说，因为大多数地区的MMT高于第90百分位，所以在第10、75和90百分位处设置节点可能会影响模型捕捉热相关RR的能力。因此，估计的热相关风险以及随之而来的热相关超额死亡率具有更大的不确定性。其次，分析是在市政级别进行的，随后向下缩放到地区级别，假设同一市政区域内各地区的基线死亡风险是均匀的，地区间的差异主要由元预测因子的值驱动。虽然这种假设简化了建模框架并有助于有意义的地理聚类——这对公共卫生干预措施很有用——但它可能掩盖了局部风险的差异性。第三，使用时间固定的地区级元预测因子限制了我们考察研究期间不断变化的局部特征对死亡率影响的能力。第四，我们仅估计了两个广泛年龄组的关联。因此，可能忽略了更细分类别之间的脆弱性差异。然而，对更多更细分类别的温度-死亡率关联进行评估将导致统计功效问题，尤其是对于年轻人群体。

最后，关于在元回归模型中使用PCs作为元预测因子有两个需要考虑的问题。首先，尽管我们纳入了五个PCs作为元预测因子，解释了78%的区域级特征方差，但元回归仍显示出市政之间的中等残余异质性。这表明我们模型中包含的PC集和潜在的区域级变量并未完全解释市政间基线风险的变化。其次，尽管PCA有效降低了数据复杂性，但它阻止了对特定脆弱性因素的独立效应的评估。然而，不同树木覆盖率、气体使用情况和居住密度的地区之间死亡率没有显著差异，这可能表明至少在荷兰，复合指标可能比单一变量更能捕捉对温度暴露的脆弱性的多因素性质。

地区级的温度相关标准化死亡率地图显示，城市化地区的热相关死亡率较高，而北部、西南部和东南部的寒冷相关死亡率较高。75岁以上年龄组的死亡风险比年轻组更为显著。低收入和高收入地区的热相关和寒冷相关标准化死亡率之间存在微小差异，而不同树木覆盖率、气体使用情况或居住密度的地区之间没有差异。重要的是，使用小规模数据使我们能够在细小的空间分辨率上进行风险估计，揭示了在更广地理尺度分析中可能被忽略的局部异质性。这些发现可能有助于指导荷兰的公共卫生政策，并支持未来对气候变化下温度相关健康风险的评估。

作者感谢Danny Houthuijs、Frank den Hertog、Cor Jacobs、Guus Velders、Lisbeth Hall、Fons van der Lucht和Ronald Gijsen在研究概念化过程中所做的贡献。