陳金建 楊 軍 陳仁超 劉起勇 歐春泉

【提 要】 目的 介紹多元meta隨機(jī)效應(yīng)模型評(píng)價(jià)多個(gè)城市暴露-反應(yīng)非線性關(guān)系的理論，并通過實(shí)例加以詳細(xì)闡述。 方法 多元meta隨機(jī)效應(yīng)模型廣泛應(yīng)用于處理多重相關(guān)參數(shù)的復(fù)雜效應(yīng)研究，本文基于該方法評(píng)估多個(gè)地區(qū)逐日氣溫對(duì)人群非意外死亡人數(shù)的非線性影響。建模過程主要基于二階段策略，首先采用分布滯后非線性模型(distributed lag non-linear model,DLNM)構(gòu)建氣溫-死亡-滯后效應(yīng)的二維交叉基，分別估計(jì)各個(gè)城市氣溫與死亡的非線性暴露-反應(yīng)關(guān)系，并提取各城市回歸系數(shù)與協(xié)方差矩陣；然后，應(yīng)用多元meta隨機(jī)效應(yīng)模型對(duì)各城市估計(jì)的參數(shù)進(jìn)行合并，得到氣溫對(duì)人群非意外死亡的綜合效應(yīng)。結(jié)果 我國16個(gè)城市的氣溫對(duì)人群非意外死亡人數(shù)存在滯后影響，總體效應(yīng)關(guān)系呈現(xiàn)U形，提示低溫與高溫對(duì)人群健康影響較大，且低溫效應(yīng)大于高溫效應(yīng)。結(jié)論 多元隨機(jī)效應(yīng)meta模型可擬合不同研究間多個(gè)復(fù)雜的相關(guān)參數(shù)，非常適合于暴露-反應(yīng)呈現(xiàn)非線性的效應(yīng)綜合研究。

氣候變化背景下，環(huán)境危險(xiǎn)因素的健康影響愈來愈受到各級(jí)政府、公眾與科學(xué)家關(guān)注。以往大量研究表明氣溫與健康結(jié)局(如：發(fā)病率與死亡率)呈現(xiàn)V、U或J型的非線性的暴露-反應(yīng)關(guān)系[1-3]。即當(dāng)氣溫高于或低于某個(gè)臨界值時(shí)，健康結(jié)局發(fā)生率會(huì)隨著溫度的升高或降低而上升。由于不同地區(qū)間氣候、人群適應(yīng)性與社會(huì)經(jīng)濟(jì)結(jié)構(gòu)存在較大差距，難以將單個(gè)城市的研究結(jié)論推廣至其他地區(qū)。因此有必要采用meta分析的技術(shù)將不同城市或地區(qū)的單獨(dú)研究結(jié)果進(jìn)行合并，從而得到更有說服力的證據(jù)。單變量的meta分析模型僅針對(duì)溫度兩端的效應(yīng)(低溫與高溫)進(jìn)行合并，無法獲得多個(gè)地區(qū)氣溫與健康的綜合反應(yīng)關(guān)系曲線，同時(shí)對(duì)溫度滯后效應(yīng)的擬合不夠充分，未能客觀反映氣溫對(duì)死亡的非線性效應(yīng)。2012年，Gasparrini[3-4]將傳統(tǒng)多元meta分析的技術(shù)引入氣溫健康效應(yīng)研究中，重新闡述該理論實(shí)現(xiàn)過程，并制作了”mvmeta”的R軟件包，為實(shí)現(xiàn)多個(gè)地區(qū)暴露-反應(yīng)非線性效應(yīng)的合并提供便利。

因此，本研究將詳細(xì)介紹多元隨機(jī)效應(yīng)meta分析的理論過程，并以中國16個(gè)城市為實(shí)例應(yīng)用該模型進(jìn)行分析。

原理與方法

暴露-反應(yīng)關(guān)系的合并通常采用標(biāo)準(zhǔn)的二階段建模分析策略：第一階段得到單個(gè)研究城市的效應(yīng)曲線，本文運(yùn)用分布滯后非線性模型(distributed lag non-linear model,DLNM)擬合各城市非線性的氣溫與人群死亡結(jié)局的暴露-反應(yīng)關(guān)系[5]；然后第二階段使用擴(kuò)展的多元meta分析方法綜合多個(gè)城市的非線性暴露-反應(yīng)關(guān)系。

1.第一階段模型

(1)廣義相加模型基本結(jié)構(gòu)[5-7]：

(1)

其中g(shù)為log函數(shù)；μt≡E(Y)，Y可為多種指數(shù)分布族(exponential family)，如正態(tài)分布、二項(xiàng)分布、Gamma分布、Poisson分布等；在評(píng)估環(huán)境危險(xiǎn)因素健康效應(yīng)的時(shí)間序列研究中，因變量Yt(t=1,2,…,n)通常是某一地區(qū)人群中某種陽性事件的逐日累計(jì)人數(shù)(如每日死亡人數(shù)、每日患病人數(shù)等)，與當(dāng)?shù)厝丝诨鶖?shù)相比，其為小概率事件，且時(shí)常出現(xiàn)過度離散現(xiàn)象(方差大于均數(shù))，故通常采用廣義Poisson分布；自變量xj通常是同期的逐日空氣污染物濃度、溫度、相對(duì)濕度等環(huán)境因子，uk表示其他混雜因素，βj、γk為相應(yīng)的系數(shù)；fj表示自變量xj的各種基函數(shù)，用于擬合非線性的自變量，常用的參數(shù)類型的基函數(shù)有正交函數(shù)、多項(xiàng)式、線性閾值函數(shù)、自然三次樣條和B樣條等，非參數(shù)類型的基函數(shù)有局部回歸函數(shù)和懲罰樣條函數(shù)。

(2)分布滯后非線性模型

DLNM模型以廣義線性模型和廣義相加模型等傳統(tǒng)模型的思想為基礎(chǔ)，利用交叉基(cross-basis)過程，可同時(shí)靈活擬合暴露與健康非線性及其滯后效應(yīng)，因此非常適合于溫度健康效應(yīng)研究。具體地，定義長度為N的向量X=[x1,…,xt,…,xN]T,L是最大滯后天數(shù)，l=[0,…,l,…,L]T，設(shè)定φ1(.),φ2(.),…,φvx(.)是擬合暴露-反應(yīng)維度的基函數(shù)；ψ1(.),ψ2(.),…,ψvl(.)是擬合暴露滯后效應(yīng)維度的基函數(shù)；vx,vl是暴露/滯后基函數(shù)的維度。可得如下公式[5]：

(2)

β=Mη

V(β)=MV(η)MT

其中：M=1(L+1)C?Ivx是簡化矩陣；β=(β1,β2,…,βvx)是簡化后的參數(shù)

2.第二階段模型

(1)多元meta模型

(3)

其中Si+Ψ=Σi。公式(3)包含暴露在研究內(nèi)與研究間獨(dú)立的效應(yīng)變異，前者可由抽樣誤差Nk(θi,Si)表示，服從k維正態(tài)分布；后者來自Nk(θ,Ψ)總體,Ψ為研究間待估計(jì)的協(xié)方差矩陣。當(dāng)存在有p個(gè)可用來預(yù)測結(jié)局的協(xié)變量xi=[x1i,x2i,…,xpi]T時(shí)，模型(3)可進(jìn)一步擴(kuò)展為多元meta回歸模型：

(4)

其中：

Xi(k×kp)=I(k)?xiT=

另外，若假定各個(gè)研究與總體的真實(shí)效應(yīng)之間差異的隨機(jī)效應(yīng)為ui=Nk(0,Ψ)，則模型(4)可用另一種形式表示：

(5)

(2)系數(shù)的估計(jì)

模型(3)系數(shù)的估計(jì)主要基于最大似然和限制性最大似然法。首先，模型(3)的邊際對(duì)數(shù)似然函數(shù)為：

(6)

其中：Σi(ξ)簡記為Σi，n為觀測值總數(shù)，若無缺失時(shí)其等于km。

(7)

但若Ψ未知時(shí)，另一種基于n-q個(gè)線性誤差對(duì)比的估計(jì)法更為合理，q為β中固定效應(yīng)系數(shù)個(gè)數(shù)，采用限制性最大似然法得到：

(8)

(3)假設(shè)檢驗(yàn)與模型預(yù)測

隨機(jī)效應(yīng)meta分析關(guān)注點(diǎn)在于不同研究間是否存在顯著的異質(zhì)性，其原假設(shè)H0為:Ψ=0，備擇假設(shè)H1:Ψ≠0。將Cochran Q異質(zhì)性檢驗(yàn)擴(kuò)展到多元情況得到：

(9)

進(jìn)一步地，可得到以下兩個(gè)評(píng)價(jià)異質(zhì)性的具體指標(biāo)：

(10)

其中H2為異質(zhì)性大?。籌2為協(xié)變量不能解釋部分占比。

實(shí)例分析

從中國疾病預(yù)防控制中心獲得我國16個(gè)省會(huì)城市2007-2013年逐日非意外死亡人數(shù)數(shù)據(jù)，同時(shí)期逐日平均溫度、相對(duì)濕度與大氣氣壓資料來自于中國氣象數(shù)據(jù)共享網(wǎng)。16個(gè)省會(huì)城市包括哈爾濱、長春、沈陽、北京、天津、銀川、濟(jì)南、鄭州、上海、南京、合肥、成都、武漢、重慶、長沙、廣州，各城市氣象與死亡數(shù)據(jù)統(tǒng)計(jì)描述見我們前期發(fā)表的文章[8-9]。首先采用分布滯后非線性模型擬合各城市溫度-死亡之間的非線性暴露-反應(yīng)關(guān)系，其中溫度和滯后維度均采用自由度為5的三次樣條函數(shù)，同時(shí)采用三次樣條函數(shù)控制逐日濕度、逐日大氣氣壓、時(shí)間變量的混雜，這些變量自由度分別為3、3、8/年；以上自由度采用Q-AIC最小化準(zhǔn)則進(jìn)行篩選得到；溫度最大滯后天數(shù)設(shè)定為21天[4]。得到各城市單獨(dú)效應(yīng)后，采用多元meta隨機(jī)效應(yīng)模型進(jìn)行效應(yīng)合并。

表1 各個(gè)城市第一階段模型的系數(shù)估計(jì)

*：b1～b5是由于擬合氣溫變量時(shí)自由度為5，待估計(jì)的參數(shù)為5個(gè)

表1展示了第一階段中估計(jì)得到的氣溫的模型系數(shù)。在得到各個(gè)城市系數(shù)估計(jì)及相應(yīng)的協(xié)方差矩陣后，可利用第二階段多元meta隨機(jī)效應(yīng)進(jìn)行效應(yīng)合并。

圖1 16個(gè)省會(huì)城市氣溫對(duì)非意外死亡的效應(yīng)曲線

圖1提示16個(gè)城市溫度-死亡效應(yīng)曲線呈U型，最適(死亡風(fēng)險(xiǎn)最低)溫度為17℃，低溫與高溫均使人群死亡風(fēng)險(xiǎn)上升。相對(duì)于最適溫度，第1百分位數(shù)氣溫(低溫)與第99百分位數(shù)氣溫(高溫)的相對(duì)危險(xiǎn)度分別為2.039(95%CI:1.744～2.384)與1.390(95%CI:1.171～1.650)。進(jìn)一步，以秦嶺-淮河南北分界線為界，得到南方低溫與高溫的相對(duì)危險(xiǎn)度分別為2.309(95%CI:1.773～2.995)與1.175(95%CI:1.030～1.341)，北方分別為(95%CI:1.521～2.277)與1.566(95%CI:1.268～1.2.099)。

討 論

本文介紹的多元meta隨機(jī)效應(yīng)模型具有綜合效應(yīng)關(guān)系的功能。它是傳統(tǒng)的標(biāo)準(zhǔn)兩階段設(shè)計(jì)的擴(kuò)展，可對(duì)信息量更豐富的多元關(guān)聯(lián)參數(shù)進(jìn)行綜合估計(jì)。Sauerbrei和Royston[10]提出另一種通過對(duì)暴露的效應(yīng)進(jìn)行網(wǎng)格化，以期構(gòu)建總體非線性關(guān)系的meta分析方法，稱為meta平滑(meta-smoothing)。該方法應(yīng)用前提是假設(shè)各網(wǎng)格的效應(yīng)為獨(dú)立的，但通常情況該假設(shè)條件難以滿足，因此在暴露效應(yīng)值的標(biāo)準(zhǔn)誤和置信區(qū)間的估計(jì)可能存在偏倚。多元meta隨機(jī)效應(yīng)模型通過估計(jì)構(gòu)建協(xié)方差矩陣來允許各城市的效應(yīng)估計(jì)之間存在相關(guān)性?？舍槍?duì)模型擬合與比較、異質(zhì)性和預(yù)測分別進(jìn)行統(tǒng)計(jì)檢驗(yàn)，并且在參數(shù)估計(jì)方面有更好的統(tǒng)計(jì)性能，比如可通過估計(jì)研究間協(xié)方差矩陣結(jié)構(gòu)來提高精度。但也有一定的不足[3]，需要估計(jì)的模型參數(shù)繁多，包括k維固定效應(yīng)參數(shù)向量，p個(gè)個(gè)體的參數(shù)個(gè)數(shù)為kp，協(xié)方差參數(shù)個(gè)數(shù)為k(k+1)/2；另外，要求所有個(gè)體研究取相同的暴露-反應(yīng)函數(shù)，可能造成部分信息的缺失。

實(shí)例應(yīng)用中探究了中國16個(gè)城市的每日氣溫和非意外死因之間的關(guān)系，得出氣溫-死亡綜合的暴露-反應(yīng)關(guān)系曲線呈U形，提示低溫與高溫會(huì)在一定程度上增大人群死亡率，并且低溫比高溫的影響更大，這與以往研究相一致[1，11-12]。

綜上所述，多元meta分析可非常適合應(yīng)用于多個(gè)地區(qū)暴露與健康呈現(xiàn)復(fù)雜的非線性關(guān)系的環(huán)境流行病學(xué)研究。另外，該方法還可應(yīng)用于多站點(diǎn)臨床試驗(yàn)、多中心隊(duì)列、多個(gè)癌癥登記研究或多城市時(shí)間序列研究等；其研究目的還包括合并生存曲線[13-14]、縱向剖面[15]，ROC曲線[16]和熱浪效應(yīng)[17]等不同研究問題，具有非常良好的應(yīng)用前景。