基于分布滯后非線性模型的歸因風(fēng)險評估方法及應(yīng)用*

寧波市疾病預(yù)防控制中心（315010） 谷少華 賀天鋒陸蓓蓓 徐倩倩 梅秋紅 張思恒

基于分布滯后非線性模型的歸因風(fēng)險評估方法及應(yīng)用*

寧波市疾病預(yù)防控制中心（315010） 谷少華 賀天鋒△陸蓓蓓 徐倩倩 梅秋紅 張思恒

目的介紹基于分布滯后非線性模型的歸因風(fēng)險評估方法，并運用該方法評估寧波市氣溫暴露造成人群死亡的歸因風(fēng)險。方法分布滯后非線性模型通過交叉基函數(shù)實現(xiàn)同時描述因變量在自變量維度與滯后維度的分布，使其能夠同時評估出暴露因素的滯后效應(yīng)和非線性效應(yīng)。收集寧波市2009－2014年人群死亡和氣象資料，利用時間序列分析結(jié)合分布滯后非線性模型，評估氣溫造成人群死亡的歸因死亡人數(shù)和人群歸因分值。結(jié)果寧波市2009－2014年日均氣溫與總死亡的累積暴露－反應(yīng)關(guān)系曲線近似呈L型，26℃為最適宜溫度。歸因于氣溫暴露造成的死亡人數(shù)為29037例（95%CI：19181～38074），占總死亡的13．39%（95%CI：9．19%～17．49%）。低溫的歸因風(fēng)險大于高溫，歸因死亡人數(shù)分別為27088例和1977例，歸因分值分別為12．49%和0．91%。結(jié)論無論高溫或低溫均與人群死亡增加相關(guān)，低溫的歸因風(fēng)險更大。

歸因風(fēng)險 分布滯后非線性模型 氣溫 死亡

隨著極端天氣和空氣污染事件不斷增多，研究者越來越關(guān)注氣溫或大氣顆粒物等環(huán)境暴露因子對人群健康的影響，其中開展環(huán)境因素的歸因風(fēng)險評估是流行病學(xué)研究中的重要環(huán)節(jié)。根據(jù)暴露因素與健康結(jié)局的關(guān)聯(lián)程度，結(jié)合暴露人群數(shù)量和暴露水平，能夠評估出暴露因素造成的人群歸因風(fēng)險，這個指標(biāo)比以往研究中采用的相對危險度（relative risk，RR）或比值比（odds ratio，OR）更能反映出整體疾病負(fù)擔(dān)［1］。但是，由于氣溫或大氣顆粒物與健康結(jié)局的暴露－反應(yīng)關(guān)系往往呈非線性關(guān)系，同時效應(yīng)又存在滯后性和持續(xù)性，這些特征對準(zhǔn)確評估歸因風(fēng)險提出了挑戰(zhàn)［2－3］。Gasparrini等基于分布滯后非線性模型提出了新的歸因風(fēng)險評估方法，能夠同時擬合暴露－反應(yīng)關(guān)系的非線性效應(yīng)和滯后效應(yīng)，為開展環(huán)境因素的歸因風(fēng)險評估提出了新的思路［4－5］。本研究將對此方法的基本理論進(jìn)行介紹，并運用此方法開展寧波市氣溫對人群死亡影響的歸因風(fēng)險評估。

原理與方法

1．歸因風(fēng)險計算方法

歸因風(fēng)險的基礎(chǔ)指標(biāo)是人群歸因分值（attribute fraction，AF），表示如果人群不再暴露于某風(fēng)險因子，相應(yīng)健康結(jié)局減少的數(shù)量占該健康結(jié)局的比例［6］。如果已知健康結(jié)局的總數(shù)和人群歸因分值，則可計算出歸因人數(shù)（attribute numbers，AN），這些指標(biāo)同樣可理解為對風(fēng)險因子采取干預(yù)措施后可以達(dá)到的效果。對于最簡單的二分類暴露模式（暴露或非暴露），效應(yīng)指標(biāo)可用相對危險度（RR）表示，人群歸因分值基本計算公式如下［1，6］：

其中Pe表示人群暴露于待研究因素的比例，對于全人群暴露（如空氣污染或氣溫）的風(fēng)險評價，可認(rèn)為Pe為1。公式簡化為［1，5］：

βx為暴露因素x的回歸系數(shù)。

由于人群暴露于某風(fēng)險因子通常是一個持續(xù)的過程，評估時可以將暴露分為不同的水平，分別計算相對于基線暴露水平時人群的風(fēng)險，然后進(jìn)行風(fēng)險累加。則公式修正為［5－6］：

RRi為和基線水平相比，各暴露水平下的相對危險度；βxi為暴露水平為i時的效應(yīng)參數(shù)，當(dāng)暴露－反應(yīng)關(guān)系為線性時，可表示為回歸系數(shù)（β）和暴露水平（xi）的乘積。

2．分布滯后非線性模型（distributed lag non-linear model，DLNM）

由于人群健康不僅受到當(dāng)天環(huán)境因素的影響，還可能與幾天前的暴露水平有關(guān)，為了評估這種滯后效應(yīng)，Zanobetti［3］和Armstrong［2］等人開始將分布滯后模型運用于空氣污染或氣溫等環(huán)境因素短期效應(yīng)研究中。在以往研究基礎(chǔ)上，Gasparrini等［4，7］利用交叉基（cross-basis）函數(shù)闡述了分布滯后非線性模型的理論，本文將對此模型進(jìn)行簡要介紹。

模型的基本結(jié)構(gòu)如下：

g為鏈接函數(shù)族；Y為結(jié)局變量；xi為自變量；fj為自變量xi的各種基函數(shù)，如線性閾值函數(shù)或樣條函數(shù)等；μk為其他混雜因素；βj和γk為方程中相應(yīng)的參數(shù)。

如果自變量和因變量的暴露－反應(yīng)關(guān)系用函數(shù)f（x）表示，滯后－反應(yīng)關(guān)系用w（l）表示，將兩個函數(shù)合并即可得到雙維度的暴露－滯后－反應(yīng)關(guān)系函數(shù)f·w（x，l），計算過程中f·w（x，l）可以簡化表達(dá)為βxt，l。因此，自變量不同滯后時間的累積風(fēng)險s（xt；η）計算公式如下：

η為方程中相應(yīng)的參數(shù)，L為最大滯后時間。

分布滯后非線性模型通過交叉基函數(shù)給暴露－反應(yīng)關(guān)系添加滯后維度，實現(xiàn)同時描述因變量在自變量維度與滯后維度的分布，使其能夠同時評估出暴露因素的滯后效應(yīng)和非線性效應(yīng)。此外，通過累積效應(yīng)的計算能夠發(fā)現(xiàn)效應(yīng)最低時的暴露水平，稱最適宜暴露水平，可以作為歸因風(fēng)險評估的基線水平［5］。

3．基于DLNM的歸因風(fēng)險計算方法

Gasparrini等［5］認(rèn)為可以通過兩種方式計算滯后效應(yīng)，第一種是“從后往前看”，認(rèn)為第t天的風(fēng)險是前一段時間（t－l0，…，t－L）暴露效應(yīng)的累加，可稱為“后向視角（backward perspective）”；第二種是“從前往后看”，認(rèn)為第t天的暴露造成了未來一段時間（t+l0，…，t+L）的風(fēng)險，稱為“前向視角（forward perspective）”。結(jié)合分布滯后非線性模型的原理，歸因分值和歸因人數(shù)的計算公式修改如下［5］：

后向視角：

前向視角：

L為暴露因素的最長滯后時間；nt為第t日的人群某健康結(jié)局發(fā)生總數(shù)。

“后向視角”是研究中常用的對滯后效應(yīng)的解釋［2，4］，但是其計算過程較為復(fù)雜。“前向視角”的原理和計算均較為簡單，但是由于在計算歸因人數(shù)時，nt采用了滯后期間總健康結(jié)局人數(shù)的平均值，因此可能會低估實際的風(fēng)險大?。?］。實例分析中將主要報告“后向視角”的計算結(jié)果，同時比較兩種計算方法的差異。

應(yīng) 用

寧波市2009年1月1日至2014年12月31日人群逐日死亡數(shù)據(jù)來源于寧波市疾病預(yù)防控制中心，同期的氣象資料來源于寧波市氣象局。利用Excel 2013軟件進(jìn)行數(shù)據(jù)整理，變量包括逐日的人群總死亡數(shù)（death）、日均氣溫（temp）、日均相對濕度（rh）、日均氣壓（press）、日期（date）、時間變量（time）、星期幾（dow）等，資料保存名稱為“NBdeath2009－2014．csv”，數(shù)據(jù)整理格式見表1。

表1 寧波市2009－2014年數(shù)據(jù)整理情況（前10條記錄）

利用R軟件（3．1．1版本）中的“dlnm”程序包評估氣溫對總死亡的歸因風(fēng)險，首先通過分布滯后非線性模型計算出氣溫與死亡的暴露－滯后－反應(yīng)關(guān)系，再根據(jù)累積效應(yīng)最小判斷出最適宜氣溫作為基線評估水平，結(jié)合寧波市的總死亡人數(shù)和溫度分布范圍，最終評估出氣溫造成的人群歸因分值和歸因總死亡人數(shù)［5，8］?；灸Ｐ瓦x擇廣義線性模型，通過時間變量（time＝1，2，3，…，2191）控制死亡人數(shù)的長期變化趨勢和季節(jié)趨勢，并控制了相對濕度、氣壓、星期幾效應(yīng)等混雜因素的影響。根據(jù)以往文獻(xiàn)，選擇日均氣溫作為氣溫的代表指標(biāo)，基函數(shù)選用自然三次樣條函數(shù)（nature cubic spline），最長滯后時間為14天，時間變量的自由度7/年，交叉基函數(shù)中的自由度利用赤池信息準(zhǔn)則（Akaike information criterion，AIC）確定［4－5，8］。主要編程代碼如下：

研究發(fā)現(xiàn)寧波市2009－2014年日均氣溫和總死亡的關(guān)聯(lián)有統(tǒng)計學(xué)意義（P＜0．05），累積暴露－反應(yīng)關(guān)系曲線近似呈L型，26℃為最適宜溫度，約在全年氣溫的第78百分位。高溫的效應(yīng)出現(xiàn)早，持續(xù)時間短，滯后時間約為0～2天；低溫的效應(yīng)出現(xiàn)晚，持續(xù)時間長，滯后時間約為2～10天。以最適宜溫度作為參考，氣溫第1百分位和第99百分位累積0～14天的相對危險度分別為1．51（95%CI：1．32～1．72）和1．15（95%CI：1．06～1．25），低溫的效應(yīng)大于高溫，見圖1。

圖1 日均氣溫與總死亡的暴露－反應(yīng)關(guān)系

采用“后向視角”評估方法，日均氣溫26℃作為基線暴露水平，累積滯后0～14天時，寧波市2009－2014年歸因于氣溫暴露造成的死亡人數(shù)為29037例（95%CI：19181～38074），占總死亡的13．39%（95%CI：9．19%～17．49%）；高溫和低溫的歸因總死亡人數(shù)分別為1977例和27088例，歸因分值分別為0．91%和12．49%。如圖2所示，當(dāng)只分析高溫與總死亡的關(guān)系時，“前向視角”計算的歸因死亡人數(shù)隨氣溫波動而波動，且均大于或等于零；而“后向視角”的計算結(jié)果則比氣溫的變化較為滯后，數(shù)值波動更大，并出現(xiàn)小于零的情況。但是，兩種評估方法的最終歸因風(fēng)險結(jié)果非常接近，“前向視角”的結(jié)果偏小，詳情見表2。

討 論

本研究發(fā)現(xiàn)無論高溫或低溫均與人群死亡增加相關(guān)，但是歸因于低溫效應(yīng)的比例（12．49%）明顯大于高溫的作用（0．91%）。兩者效應(yīng)的差異可能是由于影響機制不同，高溫能夠引起機體心率升高、血液粘度增加、水鹽代謝失調(diào)等改變，效應(yīng)出現(xiàn)快且持續(xù)時間短；而低溫則主要引起血管收縮、血壓改變、炎癥反應(yīng)等，效應(yīng)持續(xù)時間較長，造成的影響可能也更大［8－9］。同時部分文獻(xiàn)也發(fā)現(xiàn)，隨著人們采取更多的適應(yīng)措施，一些地區(qū)高溫的效應(yīng)可能正在不斷降低，而低溫的效應(yīng)則相對穩(wěn)定［10－11］。氣候條件、經(jīng)濟水平、人群特征、生活習(xí)慣等原因均有可能導(dǎo)致不同地區(qū)氣溫歸因風(fēng)險差異，寧波市氣溫的人群歸因分值（13．39%）高于中國的其他城市（11．00%）［8］，這提示寧波市人群對氣溫可能更為敏感，當(dāng)?shù)卣枰獙Υ送度敫嗟年P(guān)注。

圖2 寧波市2009年6－9月高溫的歸因總死亡人數(shù)和氣溫分布趨勢

盡管氣溫的累積效應(yīng)均為正值，但是“后向視角”評估出高溫造成的歸因風(fēng)險出現(xiàn)了負(fù)值，這可能是由“收獲效應(yīng)（harvest effect）”導(dǎo)致的。一次極端高溫短時間內(nèi)會造成大量人群死亡，其中多為老年人、慢性病患者等脆弱人群，這些人群的死亡使得一段時間內(nèi)整個人群中脆弱人群的比例降低，人群對氣溫變化的敏感程度也隨之降低，稱為“收獲效應(yīng)”，這種現(xiàn)象也在以前的多次研究中被提到［5，8］?！扒跋蛞暯恰眲t無法觀測到這種現(xiàn)象，由于其可能低估風(fēng)險大小，因此不適宜用于整體歸因風(fēng)險評估。但是，“前向視角”的原理和計算方法簡單，容易評估出不同暴露范圍（如高溫或低溫）對健康影響的差異，因此比較適用于分段風(fēng)險評估［5］。

表2 不同研究方法中的人群歸因分值（AF）和歸因總死亡人數(shù)（AN）及其95%置信區(qū)間

基于分布滯后非線性模型的歸因風(fēng)險評估方法能夠同時考慮暴露因素的滯后效應(yīng)和非線性效應(yīng)，這種方法同樣可以應(yīng)用于評估空氣污染或其他環(huán)境因素對健康的影響，還可推廣到任何探究預(yù)測變量與結(jié)局關(guān)系的時間序列研究［7］。通過疾病負(fù)擔(dān)評估能夠提高人們對暴露因素危害的認(rèn)識，并為決策者進(jìn)行健康效益分析，制定相應(yīng)的防護(hù)策略提供理論依據(jù)。

［1］Steenland K，Armstrong B．An overview of methods for calculating the burden of disease due to specific risk factors．Epidemiology，2006，17（5）：512-519．

［2］Armstrong B．Models for the relationship between ambient temperature and daily mortality．Epidemiology，2006，17（6）：624-631．

［3］Zanobetti A，Wand MP，Schwartz J，et al．Generalized additive distributed lag models：quantifying mortality displacement．Biostatistics，2000，1（3）：279-292．

［4］Gasparrini A，Armstrong B，Kenward MG．Distributed lag non-linear models．Statistics in medicine，2010，29（21）：2224-2234．

［5］GasparriniA，Leone M．Attributable risk from distributed lagmodels．BMC medical research methodology，2014，14（1）：55．

［6］Baker D，Nieuwenhuijsen M J．Environmental Epidem iology：Study Methods and Application．Oxford University，2008．

［7］楊軍，歐春泉，丁研，等．分布滯后非線性模型．中國衛(wèi)生統(tǒng)計，2012，29（5）：772-773．

［8］Gasparrini A，Guo Y，Hashizume M，et al．Mortality risk attributable to high and low ambient temperature：a multicountry observational study．Lancet，2015，386（9991）：369-375．

［9］Basu R．High ambient temperature and mortality：a review of epidem iologic studies from 2001 to 2008．Environmental health：a global access science source，2009，8（1）：40．

［10］Barnett AG．Temperature and cardiovascular deaths in the US elderly：changes over time．Epidem iology，2007，18（3）：369-372．

［11］Gasparrini A，Guo Y，Hashizume EM，et al．Temporal Variation in Heat-Mortality Associations：A Multicountry Study．Environmental health perspectives，2015，123（11）：1200-1207．

（責(zé)任編輯：鄧 妍）

M easures and App lication for Attributable Risk from Distributed Lag Non-linear M odel

Gu Shaohua，He Tianfeng，Lu Beibei，et al
（Ningbo Municipal Center for Disease Control and Prevention（315010），Ningbo）

ObjectiveTo introducemeasures of attributable risk from distributed lag non－linearmodel（DLNM），and to apply thesemethods for estimating themortality risk attributable to outdoor temperature in Ningbo city．MethodsDLNM is based on a cross-basis function that describes simultaneously the shape of the relationship along both the space of the predictor and the lag dimension of its occurrence，and could assess the potentially non-linear and lag effects．The daily data on mortality and meteorological factorswere collected from 2009 to 2014 in Ningbo city．A time series study using a DLNM was used to estimate the attributable number and fraction to the effectof temperature onmortality．ResultsThe overall cumulative exposure-response curve between temperature and mortality was L-shaped at lag 0～14 days，and the m inimum-mortality temperature was 26℃．In total，13．39%（95%CI：9．19%～17．49%）of totalmortality was attributable to outdoor temperature，while the attributable number was 29037（95%CI：19181－38074）．More attributable deaths were due to cold，w ith a fraction of 12．49%corresponding to 27088 deaths，compared w ith 0．91%and 1977 deaths for heat．ConclusionBoth heat and cold were associated w ith an increased risk of daily mortality，and mostmortality burden were caused by cold．

Attributable risk；Distributed lag non-linearmodel；Temperature；Mortality

寧波市科技局創(chuàng)新團隊項目（編號：2012B82018）；浙江省公益技術(shù)應(yīng)用研究計劃（編號：2016C33194）；浙江省醫(yī)藥衛(wèi)生科技計劃項目（編號：2014KYA202）

△通信作者：賀天鋒，E-mail：469345174＠qq．com