表1芝加哥1987—1997年逐日死亡及环境数据格式
Table1.The data structure of daily mortality and environmental factors in Chicago between 1987 and 1997
2020, 37(1):9-14.doi:10.13213/j.cnki.jeom.2020.19485
关键词： 随时间变化的分布滞后非线性模型; 死亡; 气温; 应用; 效应
极端气候是目前人类面临的最严重威胁之一[1]。大量研究显示气温与健康效应的暴露-反应关系为非线性，例如U型、V型或者J型，同时该效应又存在时间上的滞后性[2-4]。近年来，分布滞后非线性模型（distributed lag non-linear model，DLNM）的提出极大地促进了气象条件与健康效应的研究[5]。此外，有研究利用数十年的健康和气象数据，发现两者间的暴露-反应关系并不恒定，而是随时间发生变化[6-7]。DLNM应用于小数据集时，如单个城市的单一年份，所拟合的暴露-反应关系受极端值影响较大。为追求稳定性，多数研究将数十年数据集拆分成2个或多个子数据集来研究气温健康效应的长期趋势。为了克服这一短板，Gasparrini等[8-9]提出了随时间变化的DLNM，为研究极端气候健康效应的长期趋势提供了新的方法学思路。本文将介绍此模型的基本结构，并以R软件“dlnm”程序包中内置的芝加哥市1987—1997年每日死亡数据为例，比较随时间变化的DLNM和普通DLNM的气温-死亡累积效应结果，进一步阐述随时间变化的DLNM的应用。
DLNM由Armstrong于2006年引入气温健康效应的研究中，并简要介绍了交叉基函数的基本概念[10]。Gasparrini等于2010年进一步详细阐述了DLNM的原理以及在R软件中通过“dlnm”程序包的实现方式[5]。以气温-滞后-死亡的关系为例，该模型的基本结构如下：
$\begin{array}{l}{\rm{ \mathsf{ μ} = E}}\left( {\rm{Y}} \right)\\g\left( {{\mu _t}} \right) = \alpha + s\left( {{X_t};{\rm{ \mathsf{ η} }}} \right) + \sum\nolimits_{i = 1}^J {{h_j}\left( {{u_{jt}};{r_j}} \right)} \end{array}$ |
(1) |
其中，μt为第t天的预期死亡人数；g（.）是链接函数，考虑到离散分布，一般采用类泊松；α为模型的截距，假设气温的效应最大滞后天数为L天，自变量xt=（xt，xt-1，…，xt-L）是第t天至前L天的每日气温数据所组成的行矢量，η为相应参数；s（.）为拟合xt效应的交叉基函数；ujt是一系列随时间变化的混杂变量，如日均湿度、风速和空气污染数据等，rj为相应参数，hj（.）为相应的函数。
s（.）通过DLNM来拟合气温-滞后-死亡的关系，即利用交叉基函数整合气温-死亡的非线性关系和气温的非线性滞后效应这两个维度。假设气温-死亡的非线性关系用基函数Ø1（.），…，Øνx（.）来描述，非线性滞后效应用基函数φ1（.），…，φνl（.）来描述，则DLNM可表达为：
$s\left( {{x_t};{\rm{ \mathsf{ η} }}} \right)=\sum\nolimits_{j = 1}^{{v_x}} {\sum\nolimits_{k = 1}^{{v_l}} {{r_{tj}}{c_{\rm{k}}}{\eta _{jk}}} } $ |
(2) |
其中，rtj=[Øj（xt），…，Øj（xt-L）]是通过基函数Øj（.）转换气温维度xt得到的矢量；ck=[φk（0），…，φk（L）]是通过基函数φk（.）转换滞后维度（0，…，L）得到的矢量；η=（η11，…，ηνxνl）则是交叉基函数的相应参数。基函数Øj（.）和φk（.）多采用样条函数，例如自然三次样条，νx和νt分别为相应自由度。
相较于普通DLNM，随时间变化的DLNM加入了时间变量与交叉基函数的交互项，以研究随时间变化的气温-滞后-死亡非线性关系。该设计的重要假设是气温-滞后-死亡效应是随时间线性变化的。其基本结构如下：
$\begin{array}{l}{\rm{ \mathsf{ μ} = E}}\left( {\rm{Y}} \right)\\g\left( {{\mu _t}} \right) = \alpha + s\left( {{x_t}, {\mathit{T}_\mathit{t}};{\rm{ \mathsf{ η} }}} \right) + \sum\nolimits_{i = 1}^J {{h_j}\left( {{u_{jt}};{r_j}} \right)} \\s\left( {{x_t}\mathit{, }{\mathit{T}_\mathit{t}};{\rm{ \mathsf{ η} }}} \right) = \sum\nolimits_{j = 1}^{{v_x}} {\sum\nolimits_{k = 1}^{{v_l}} {{r_{tj}}{c_{\rm{k}}}{\eta _{jk0}}} } + {T_t}\sum\nolimits_{j = 1}^{{v_x}} {\sum\nolimits_{k = 1}^{{v_l}} {{r_{tj}}{c_\mathit{k}}{{\rm{ \mathsf{ η} }}_{jk1}}} } \end{array}$ |
(3) |
交叉基函数的参数η=（η110，…，ηνxνl0）+（η111，…，ηνxνl1）。其中，（η110，…ηνxνl0）同公式（2）中η=（η11，…，ηνxνl），是交叉基函数的参数；（η111，…，ηνxνl1）则为交叉基函数s（xt；η）与时间变量Tt交互项的参数。
实例分析采用R软件“dlnm”程序包中内置的芝加哥市1987年1月1日—1997年12月31日每日居民全死因死亡及环境数据（daily mortality weather and pollution data for Chicago，chicagoNMMAPS）。该数据集包含了每日的死亡、气象、空气污染和时间数据，如日期（date）、时间变量（time）、星期几（dow）、死亡总数（death）、平均气温（temp）、平均相对湿度（rh）、可吸入颗粒物（PM10）。数据格式见表 1。
表1芝加哥1987—1997年逐日死亡及环境数据格式
Table1.The data structure of daily mortality and environmental factors in Chicago between 1987 and 1997
使用R软件3.5.1的dlnm程序包分析随时间变化的气温-滞后-死亡关系。基本模型选择广义相加模型，具体代码及相应参数的解释如下：
#载入dlnm程序包#
library（dlnm）
#载入芝加哥市1987—1997年数据（共4 018 d），计算rh和PM10滞后0~1 d的移动平均值rh01和pm01#
data < - subset（chicagoNMMAPS，year<=1997）
data＄rh01 < -filter（data$rhum，c（1，1）/2，side=1）
data＄pm01 < -filter（data$pm10，c（1，1）/2，side=1）
#建立交叉基函数bs.t，根据经验最大滞后天数设置为30 d，气温和滞后维度的自由度均为4#
kn < - quantile（data＄temp，c（0.1，0.5，0.9））
bs.t < - crossbasis（data＄temp，lag=30，argvar=list（fun="ns"，knots=kn），arglag=list（knots=logknots（30，2）））
#建立time与交叉基函数的两个交互项，并把time中心锚定在1987年7月1日（即第182天）和1997年7月1日（即第3 835天），由于采用线性假设，bs.t的系数可以代表全年的气温效应[8, 11] #
int1 < -（（data$time-182）/4018）*bs.t
int2 < -（（data$time-3835）/4018）*bs.t
#建立三个广义相加模型，通过时间变量time（自由度为每年8个）控制逐日死亡人数的季节趋势和长期趋势，同时控制rh、PM10以及dow的混杂效应，参照点cen设置为死亡相对风险最低的气温值（minimum mortality temperature，MMT）；提取相应的气温-死亡关系在滞后0~30 d的累积效应；模型fit为普通DLNM，交叉基函数bs.t的参数代表1987— 1997年的平均气温-滞后-死亡效应；模型fit1为随时间变化的DLNM，额外添加交互项int1，则bs.t的参数代表1987年7月1日的气温-滞后-死亡效应，由于采用线性交互项，bs.t的参数亦可视为1987年全年的平均效应；模型fit2为随时间变化的DLNM，额外添加交互项int2，则bs.t的参数代表1997年的气温-滞后-死亡效应#
fit < - gam（death~ns（time，length（unique（data＄year））*8）+bs.t+ns（rh01，3）+ pm01+as.factor（dow），family= quasipoisson，data=data）
cross_all < -crosspred（bs.t，model=fit，cen= MMT）
fit1 < - gam（death~ns（time，length（unique（data＄year））*8）+bs.t+ns（rh01，3）+ pm01+as.factor（dow）+ int1，family= quasipoisson，data=data）
cross_all1 < -crosspred（bs.t，model=fit1，cen= MMT1）
fit2 < - gam（death~ns（time，length（unique（data＄year））*8）+bs.t+ns（rh01，3）+pm01+as.factor（dow）+ int2，family= quasipoisson，data=data）
cross_all2 < -crosspred（bs.t，model=fit2，cen= MMT2）
为比较随时间变化DLNM的结果，本研究采用普通DLNM拟合1987—1989年、1995—1997年、1987年和1997年的累积气温-死亡关系。R代码同模型fit，为方便比较随时间变化的DLNM的结果，参照点cen的参数分别同MMT1和MMT2。温度效应采用相对危险度（relative risk， RR）及其95%可信区间（confidence interval，CI）描述。使用荟萃回归检验不同年份的冷、热效应各自的统计学差异，检验水准α=0.05[12-13]。
图 1是将1987—1997年数据作为整体分析的结果。图 1A是采用普通DLNM拟合芝加哥市1987—1997年日均气温与死亡累积0~30 d的暴露-反应关系，呈近似V型，19.2℃时死亡风险最低，为MMT。图 1B是采用随时间变化DLNM的拟合结果，显示该市1987年和1997年累积0~30 d的日均气温与死亡的暴露-反应关系，MMTs分别为18.2℃和19.8℃。与图 1A相比，图 1B的置信区间同样较窄，提示拟合结果比较稳定。
图 1日均气温与总死因死亡的累积暴露-反应关系（lag0~30 d）
以MMT为参照点，定义日均气温第1百分位（-15.8℃）和第99百分位（28.9℃）累积0~30 d的RR值分别为冷温和热温效应。1987—1997年的平均冷温效应和热温效应的RR值分别为1.53（95% CI：1.34~1.75）和1.32（95% CI：1.17~1.48）。冷温效应RR值由1987年的1.59（95% CI：1.25~2.01）降低至1997年的1.50（95% CI：1.13~1.98），但差异无统计学意义（P=0.756）。热温效应RR值由1987年的1.04（95% CI：0.85~1.28）升高至1997年的1.75（95% CI：1.39~2.21），差异有统计学意义（P=0.001）。
图 2均为采用普通DLNM对不同年份小数据集拟合的结果。图 2A是为提高模型稳定度，拟合3年数据（即1987—1989年和1995—1997年）的结果。结果显示，同一冷温（-15.8℃）效应RR值在1987—1989年间为1.33（95% CI：1.03~1.71），在1995—1997年间为1.54（95% CI：1.14~2.08），差异没有统计学意义（P=0.471）。同一热温（28.9 ℃）效应RR值在1987— 1989年间为1.24（95% CI：1.00~1.52），在1995—1997年间为2.32（95% CI：1.78~3.04），差异有统计学意义（P < 0.001）。图 2B为1987年和1997年单一年份的拟合结果，结果显示可信区间较大（特别是冷温效应），提示模型稳定度较低。
图 2采用普通DLNM拟合的日均气温与总死因死亡的累积暴露-反应关系（lag0~30 d）
随时间变化DLNM将1987—1997年作为整体进行分析，不破坏时间序列数据；而普通DLNM把时间序列数据分为不同时间段进行分析，本文选取了1987—1989年和1995—1997年两个时间段纳入分析，结果显示采用普通DLNM（图 2A）拟合3年数据的结果类似于随时间变化的DLNM（图 1B）的拟合结果，即两者都显示随着年份的延后，低温效应RR值无明显变化，而高温效应RR值升高。随时间变化DLNM的数据拟合结果（图 1B）比普通DLNM采用1987年和1997年单一年份的拟合结果（图 2B）更为稳定。该结果表明随时间变化DLNM能很好地拟合极端温度的效应随时间的变化趋势。
本研究发现芝加哥市1987—1997年在0~30d的累积气温-死亡关系为非线性，近似V型，即冷温和热温均会引起超额死亡风险。亚洲（例如中国、日本和韩国）、欧洲和澳洲（澳大利亚）等国家和地区的研究有类似的发现[14-15]。此前的研究显示，人体的自主温度调节系统不能完全抵消冷温和热温的影响，引起多种病理生理改变进而导致死亡[16-17]。例如冷温可导致血液黏稠度升高、血脂和心率异常等[18]。热温可以引起一系列高温疾病，如热中风、脱水、呼吸道炎性反应和急性肾衰竭等[19-22]。
本研究利用随时间变化的DLNM，发现在1987— 1997年间气温-死亡非线性关系并非恒定，而是随着时间的延后热温效应明显升高。该发现提示存在对热温的不适应性。此前不同国家的发现各有异同，例如一项日本的研究显示，随着时间的延后冷温和热温导致的死亡风险均下降[23]。另一项多国家研究显示，在1993—2006年的14年间，热温效应在美国和西班牙等国家明显下降，在澳大利亚升高，在韩国和英国无明显变化[8]。由于缺乏必要的个体信息，很难直接解释这种随时间变化的气温效应的机制。但是有研究推测，这可能与社会经济和地理因素有关，例如人口老龄化、地区差异、空调使用率以及经济水平等[23]。
本研究同时比较了随时间变化的DLNM、普通DLNM单一年份和普通DLNM三年数据的分析结果。普通DLNM单一年份的结果明显更不稳定，特别是冷温效应。随时间变化的DLNM的结果与普通DLNM三年数据的结果相似，均为冷温效应无明显变化而热温效应明显升高，反映了前者在气温效应随时间变化研究中的可靠性。相比较普通DLNM三年的结果，随时间变化的DLNM另一个优势在于能很好地估计某特定气温的效应随时间变化速率，特别是单一城市的小数据集。
随时间变化的DLNM是在普通DLNM基础上发展而来，因此也面临DLNM相同的问题。例如目前缺乏统一的标准确定每年几个自由度控制长期趋势和季节趋势，暴露和滞后两个维度上自由度或者节点的数目，最大滞后时间和模型拟合效果等。此外，应用随时间变化的DLNM的重要假设是气温效应随时间为线性变化。实际情况有可能出现气温的效应为非线性变化。
本研究显示，随时间变化的DLNM能很好地拟合极端气温的效应随时间的线性变化趋势。该模型还可以推广至其他时间序列研究中，探索多种暴露因素与结局的时间变化特征。
表1 芝加哥1987—1997年逐日死亡及环境数据格式
Table 1The data structure of daily mortality and environmental factors in Chicago between 1987 and 1997
图 1 日均气温与总死因死亡的累积暴露-反应关系（lag0~30 d）
Figure 1The association between daily mean temperature and all-cause mortality (over lag 0-30 days)
[注] A为采用普通DLNM拟合的1987—1997年的平均效应；B为采用随时间变化DLNM拟合的1987年（蓝色）和1997年（红色）的效应。 [Note] A shows the average association between 1987 and 1997 modelled using ordinary DLNM; B shows the associations in 1987 (blue colour) and 1997 (red colour) modelled using time-varying DLNM.
图 2 采用普通DLNM拟合的日均气温与总死因死亡的累积暴露-反应关系（lag0~30 d）
Figure 2The association between daily mean temperature and all-cause mortality (over lag 0-30 d) modelled using ordinary DLNM
[注] A为1987—1989年（蓝色）和1995—1997年（红色）的效应；B为1987年（蓝色）和1997年（红色）单一年份的效应。 [Note] A shows the associations in 1987-1989 (blue colour) and 1995- 1997 (red colour); B shows the single-year associations in 1987 (blue colour) and 1997 (red colour).| [1] |
IPCC. Summary for policymakers[M]//FIELD C B, BARROS V R, DOKKEN D J, et al. Climate change 2014: impacts, adaptation, and vulnerability. Part A: global and sectoral aspects. Cambridge, UK: Cambridge University Press, 2014: 1-32. |
| [2] |
ZHAO Q, ZHANG Y, ZHANG W, et al. Ambient temperature and emergency department visits:time-series analysis in 12 Chinese cities[J]. Environ Pollut, 2017, 224:310-316. DOI: 10.1016/j.envpol.2017.02.010 |
| [3] |
GUO Y, BARNETT AG, PAN X, et al. The impact of temperature on mortality in Tianjin, China:a case-crossover design with a distributed lag nonlinear model[J]. Environ Health Perspect, 2011, 119(12):1719-1725. DOI: 10.1289/ehp.1103598 |
| [4] |
GUO Y, GASPARRINI A, ARMSTRONG B, et al. Global variation in the effects of ambient temperature on mortality:a systematic evaluation[J]. Epidemiology, 2014, 25(6):781-789. DOI: 10.1097/EDE.0000000000000165 |
| [5] |
GASPARRINI A, ARMSTRONG B, KENWARD M G. Distributed lag non-linear models[J]. Stat Med, 2010, 29(21):2224-2234. DOI: 10.1002/sim.3940 |
| [6] |
ONOZUKA D, HAGIHARA A. Variation in vulnerability to extreme-temperature-related mortality in Japan:a 40-year time-series analysis[J]. Environ Res, 2015, 140:177-184. DOI: 10.1016/j.envres.2015.03.031 |
| [7] |
NORDIO F, ZANOBETTI A, COLICINO E, et al. Changing patterns of the temperature-mortality association by time and location in the US, and implications for climate change[J]. Environ Int, 2015, 81:80-86. DOI: 10.1016/j.envint.2015.04.009 |
| [8] |
GASPARRINI A, GUO Y, HASHIZUME M, et al. Temporal variation in heat-mortality associations:a multicountry study[J]. Environ Health Perspect, 2015, 123(11):1200-1207. DOI: 10.1289/ehp.1409070 |
| [9] |
GASPARRINI A, GUO Y, HASHIZUME M, et al. Changes in susceptibility to heat during the summer:a multicountry analysis[J]. Am J Epidemiol, 2016, 183(11):1027-1036. DOI: 10.1093/aje/kwv260 |
| [10] |
ARMSTRONG B. Models for the relationship between ambient temperature and daily mortality[J]. Epidemiology, 2006, 17(6):624-631. DOI: 10.1097/01.ede.0000239732.50999.8f |
| [11] |
ZHANG Y, YU Y, PENG M, et al. Temporal and seasonal variations of mortality burden associated with hourly temperature variability:a nationwide investigation in England and Wales[J]. Environ Int, 2018, 115(6):325-333. |
| [12] |
GUO Y. Hourly associations between heat and ambulance calls[J]. Environ Pollut, 2017, 220:1424-1428. DOI: 10.1016/j.envpol.2016.10.091 |
| [13] |
ZHAO Q, LI S, COELHO MS, et al. Geographic, demographic, and temporal variations in the association between heat exposure and hospitalization in Brazil:a nationwide study between 2000 and 2015[J]. Environ Health Perspect, 2019, 127(1):017001. DOI: 10.1289/EHP3889 |
| [14] |
MA W, WANG L, LIN H, et al. The temperature-mortality relationship in China:an analysis from 66 Chinese communities[J]. Environ Res, 2015, 137:72-77. DOI: 10.1016/j.envres.2014.11.016 |
| [15] |
GASPARRINI A, GUO Y, HASHIZUME M, et al. Mortality risk attributable to high and low ambient temperature:a multicountry observational study[J]. Lancet, 2015, 386(9991):369-375. DOI: 10.1016/S0140-6736(14)62114-0 |
| [16] |
PARSONS K. Human thermal environments:the effects of hot, moderate, and cold environments on human health, comfort, and performance[M]. 3rd ed. New York:CRC Press, 2014. |
| [17] |
GASPARRINI A, ARMSTRONG B, KOVATS S, et al. The effect of high temperatures on cause-specific mortality in England and wales[J]. Occup Environ Med, 2012, 69(1):56-61. DOI: 10.1136/oem.2010.059782 |
| [18] |
LIU C, YAVAR Z, SUN Q. Cardiovascular response to thermoregulatory challenges[J]. Am J Physiol Heart Circ Physiol, 2015, 309(11):H1793-H1812. DOI: 10.1152/ajpheart.00199.2015 |
| [19] |
WILLIAMS S, NITSCHKE M, WEINSTEIN P, et al. The impact of summer temperatures and heatwaves on mortality and morbidity in Perth, Australia 1994-2008[J]. Environ Int, 2012, 40:33-38. DOI: 10.1016/j.envint.2011.11.011 |
| [20] |
HANSEN A L, BI P, RYAN P, et al. The effect of heat waves on hospital admissions for renal disease in a temperate city of Australia[J]. Int J Epidemiol, 2008, 37(6):1359-1365. DOI: 10.1093/ije/dyn165 |
| [21] |
KALDUR T, UNT E, ÖÖPIK V, et al. The acute effects of passive heat exposure on arterial stiffness, oxidative stress, and inflammation[J]. Medicina, 2016, 52(4):211-216. DOI: 10.1016/j.medici.2016.06.001 |
| [22] |
MICHELOZZI P, ACCETTA G, DE SARIO M, et al. High temperature and hospitalizations for cardiovascular and respiratory causes in 12 European cities[J]. Am J Respir Crit Care Med, 2009, 179(5):383-389. DOI: 10.1164/rccm.200802-217OC |
| [23] |
CHUNG Y, YANG D, GASPARRINI A, et al. Changing susceptibility to non-optimum temperatures in Japan, 1972-2012:the role of climate, demographic, and socioeconomic factors[J]. Environ Health Perspect, 2018, 126(5):057002. DOI: 10.1289/EHP2546 |
[作者简介]
[收稿日期] 2019-07-19
引用格式
赵琦
,
李珊珊
,
郭玉明
.
随时间变化的分布滞后非线性模型应用介绍：以气温与死亡关系为例[J].环境与职业医学,
2020, 37(1): 9-14.
doi:10.13213/j.cnki.jeom.2020.19485.
ZHAO Qi , LI Shan-shan , GUO Yu-ming . Time-varying distributed lag non-linear model: Using temperature-mortality association as an example.Journal of Environmental & Occupational Medicine, 2020, 37(1): 9-14. doi:10.13213/j.cnki.jeom.2020.19485.