時間序列分析在空氣污染與健康領(lǐng)域的應(yīng)用及其R軟件實現(xiàn)

路 鳳 李亞偉 李成橙 陳 晨 宋士勛 郭玉明 施小明△

【提 要】 目的 探討時間序列分析在空氣污染健康影響研究領(lǐng)域的應(yīng)用與軟件實現(xiàn),為我國相關(guān)工作提供實用的操作方法參考。方法 利用美國芝加哥市1987年至2000年大氣污染與死亡數(shù)據(jù),分別采用廣義線性模型、廣義相加模型,并結(jié)合分布滯后模型,介紹時間序列分析的基本理論與R軟件實現(xiàn)步驟。結(jié)果 時間序列分析可在R軟件中方便實現(xiàn)。結(jié)論 R軟件為時間序列分析在空氣污染流行病學(xué)研究中的應(yīng)用提供了相對成熟的軟件包,在實際研究中值得推廣。

在空氣污染與健康研究領(lǐng)域,經(jīng)常需要用時間序列方法將隨時間變化的污染物暴露資料和隨時間變化的事件發(fā)生數(shù)資料聯(lián)系起來,分析人群健康結(jié)局與暴露水平之間的關(guān)系[1]。自20世紀90年代以來,空氣污染與健康關(guān)系的時間序列方法研究比較多,比較權(quán)威的有歐洲的APHEA項目、美國的NMMAPS項目、亞洲的PAPA項目等。與此同時,R軟件作為一種開源免費的軟件,已逐漸成為環(huán)境流行病學(xué)研究學(xué)者的常用工具之一。本文在回顧常用統(tǒng)計分析模型的基礎(chǔ)上,以美國芝加哥市1987-2000年大氣污染與死亡數(shù)據(jù)為例,就時間序列分析在空氣污染與健康領(lǐng)域的應(yīng)用與R軟件實現(xiàn)做一介紹。

時間序列方法介紹

時間序列分析是根據(jù)系統(tǒng)觀測得到的時間序列數(shù)據(jù),通過曲線擬合和參數(shù)估計來建立數(shù)學(xué)模型的理論和方法。在空氣污染健康效應(yīng)研究中,通常采用泊松回歸模型來估計每日污染物濃度變化與死亡人數(shù)、就診人數(shù)等健康結(jié)局指標之間的關(guān)系。目前,國際上通常使用廣義線性模型(generalized linear models,GLM)與廣義相加模型(generalized additive model,GAM)來實現(xiàn)泊松回歸過程。由于大氣污染、氣象因素等對人群健康的影響通常存在滯后性,因此要求時間序列模型不僅要反映效應(yīng)在自變量維度的變化,還要反映效應(yīng)的時間結(jié)構(gòu)(即暴露-反應(yīng)關(guān)系的時間過程)。因此,有學(xué)者將分布滯后模型與廣義線性模型或廣義相加模型結(jié)合起來用于環(huán)境因素健康效應(yīng)的定量化評估。

1.廣義線性模型(GLM) GLM是1972年由Nelder和Wedderburn提出的[2]。GLM是一般線性模型的擴展,它通過聯(lián)接函數(shù)將因變量和線性預(yù)測值關(guān)聯(lián)起來克服了線性模型要求因變量服從正態(tài)分布的限制。模型一般形式為:

(1)

在評估空氣污染與健康效應(yīng)指標的關(guān)系時,由于對于總?cè)丝趤碚f,每日死亡數(shù)、發(fā)病數(shù)是小概率事件,作為一種時間序列資料,其實際分布近似泊松分布,因此自變量和因變量之間的聯(lián)接函數(shù)通常選取對數(shù)函數(shù)來進行分析。

(2)

式(2)中，g(.)為聯(lián)接函數(shù)，f(.)為任意單變量平滑函數(shù)。在實際應(yīng)用中,模型中除了擬合普通的線性項,如大氣污染濃度外,還可以將一些與因變量之間存在復(fù)雜的非線性關(guān)系的變量,如長期趨勢、日歷效應(yīng)、氣象等混雜因素,以不同函數(shù)加和的形式擬合模型。通?？梢詫懗扇缦滦问?

(3)

3.分布滯后模型(DLM) DLM由Almon于1965年提出,并應(yīng)用于經(jīng)濟學(xué)研究[5]。DLM暴露-滯后-效應(yīng)的模型思路更加符合空氣污染導(dǎo)致健康效應(yīng)的過程,可用于分析空氣污染對人群健康急性影響的滯后效應(yīng)分布。模型一般形式為:

(4)

式(4)中,l為滯后期,q是需定義的最長滯后時間。假設(shè)空氣污染暴露的效應(yīng)存在于某一特定時間內(nèi),通過對參數(shù)q設(shè)置不同的值估計不同滯后時間的效應(yīng)。但這種簡單地將每個滯后時間與其設(shè)定的相應(yīng)參數(shù)乘積累加的方式,往往會產(chǎn)生很高的共線性而導(dǎo)致模型系數(shù)估計的不穩(wěn)定。通常采用的方法是給滯后分布強加某些限制(如假設(shè)滯后一定區(qū)間內(nèi)有相同的固定效應(yīng),或者使用多項式函數(shù)或樣條函數(shù)等連續(xù)函數(shù)來描述平滑曲線等),通過適當?shù)幕瘮?shù)轉(zhuǎn)換來提高估計結(jié)果的精確性[6]。

此外,2006年Armstrong將分布滯后非線性模型(distributed lag non-linear models,DLNM)引入氣溫健康效應(yīng)研究中[7],2010年Gasparrini等人進一步以GLM和GAM等傳統(tǒng)模型思想為基礎(chǔ),利用交叉基過程,闡述了DLNM的理論[8]。該模型同時考慮暴露因素的滯后效應(yīng)和暴露-反應(yīng)的非線性關(guān)系,突破了DLM 對暴露-反應(yīng)關(guān)系呈線性的要求。模型的基本結(jié)構(gòu)為:

(5)

式(5)中,自變量xj通常為逐日空氣污染物濃度、溫度、相對濕度等環(huán)境因素;fj表示自變量的各種基函數(shù)。Uk為其他混雜因素的線性效應(yīng),γk為其系數(shù)。雖然國際上普遍認為線性無閾值模型能夠較好地估計大氣污染與人群健康效應(yīng)之間的關(guān)系[9-10],但在實際應(yīng)用中,常常采用DLNM的思想來控制氣象因素的混雜作用。

時間序列分析在R軟件中的實現(xiàn)

1.數(shù)據(jù)簡介與分析策略

(1)數(shù)據(jù)簡介 美國芝加哥空氣污染與死亡數(shù)據(jù)是R軟件自帶的資料,該數(shù)據(jù)庫是美國NMMAPS項目的一個子集,匯總了1987-2000年間芝加哥市逐日空氣污染、氣象因素及死亡人數(shù)資料。該資料包括5114天共590252個死亡病例,變量date為日期,變量time為時間變量(取值為1,2,…,5114),變量dow代表星期幾,變量death代表總死亡人數(shù),變量temp和rhum分別代表溫度與濕度,變量pm10代表可吸入顆粒物(PM10)濃度(表1)。

(2)分析策略 分別采用廣義線性模型、廣義相加模型,并結(jié)合分布滯后模型,使用R3.4.1軟件中的tsModel、splines 、mgcv以及dlnm等程序包,以人群總死亡數(shù)為響應(yīng)變量,將大氣PM10濃度作為直線變量引入模型,同時控制氣象因素(溫度、相對濕度)、星期幾效應(yīng)、長期趨勢與季節(jié)趨勢等影響,分析單污染物模型條件下,過去兩天大氣PM10濃度的移動均值(Lag0-1)與人群日死亡數(shù)之間的關(guān)系,結(jié)果以PM10濃度每升高10μg/m3,人群總死亡增加的百分比及其95%可信區(qū)間(confidence interval,CI)表示。

表1 1987-2000年美國芝加哥市空氣污染與死亡數(shù)據(jù)庫摘錄

2.數(shù)據(jù)讀取與變量預(yù)處理

假定表1的數(shù)據(jù)以CSV格式存儲在E盤的R文件中,文件名為“chicago.csv”。通常需要以下步驟進行模型擬合前的準備:

> setwd("E:/R") #設(shè)定工作目錄

> data<-read.csv("chicago.csv ") #導(dǎo)入數(shù)據(jù)集

> data $ date<-as.Date(data $ date,"%Y/%m/%d") #將date變量轉(zhuǎn)化成日期格式

3.回歸模型的擬合

(1)廣義線性模型的實現(xiàn)

>library(tsModel);library(splines) #加載tsModel與splines包

>pm10lag01<-runMean(data $ pm10,0:1)#產(chǎn)生pm10lag01變量,賦值為過去兩天大氣PM10濃度的移動均值

> fit1<-glm(death～pm10lag01+ns(temp,3)+ns(rhum,3)+ns(time,7*14)+ dow,family=quasipoisson(),data) #GLM模型擬合

pm10lag01為過去兩天大氣PM10濃度的移動均值(下同);以自然立方樣條函數(shù)ns( )控制時間序列資料中溫度、相對濕度、時間等非線性自變量的混雜。溫度、相對濕度等氣象因素的自由度(degree of freedom,df)取值為3,時間變量的自由度取7/年。

(2)廣義相加模型的實現(xiàn)

> library(mgcv);library(tsModel) #加載tsModel與splines包

> fit2<-gam(death～ pm10lag01+ s(temp,k=3+1) +s(time,k=7 * 14+1) +s(rhum,k=3+1) +dow,family=quasipoisson(),data) #GAM模型擬合

S( )表示平滑函數(shù),常采用懲罰樣條方法進行參數(shù)的平滑;k等于df+1。此處溫度、相對濕度k值均為4,時間變量的k值取7/年×14年+1。

(3)分布滯后模型的實現(xiàn)

氣溫對健康影響的滯后性已得到公認。傳統(tǒng)的GLM與GAM模型在分析空氣污染與健康效應(yīng)之間的關(guān)系時,只考慮到當天氣溫的影響,沒有考慮其他滯后時間氣溫的混雜作用。R軟件可通過以下命令首先建立氣溫的交叉基矩陣,然后進行回歸模型的擬合:

>library(dlnm);library(tsModel);library(splines) #加載dlnm、tsModel與splines包

>cb.temp <-crossbasis(data$temp,lag=15,argvar=list(fun="ns",df=3,cen=21),arglag=list(fun="strata")) #建立氣溫的交叉基

> fit3<-glm(death～ pm10lag01+cb.temp+ns(rhum,3)+ns(time,7*14)+ dow,family=quasipoisson(),data) #模型擬合

lag設(shè)置變量的最長滯后時間,argvar=list( )、arglag=list( )分別為針對變量及其滯后時間的參數(shù)設(shè)置,fun表示調(diào)用的函數(shù),常用的函數(shù)包括ns(自然立方樣條函數(shù))、bs(B樣條函數(shù))、thr(閾值函數(shù))等;cen設(shè)置參考值。此處假定氣溫的最長滯后時間lag為15天,參考值取21℃,且氣溫在每個滯后時間的效應(yīng)一致。

需要注意的是,上述三種模型中沒有統(tǒng)一的參數(shù)設(shè)置要求,研究者應(yīng)結(jié)合自己的數(shù)據(jù)特征選擇相應(yīng)的平滑函數(shù)、自由度取值等。目前可通過經(jīng)驗法、偏自相關(guān)函數(shù)(PACF)、Akaike信息標準(Akaike’s information criterion,AIC)以及廣義交叉檢驗(general cross validation,GCV)等方法確立最佳擬合模型。有學(xué)者通過改變模型中的參數(shù)設(shè)置進行敏感性分析,觀察其對回歸系數(shù)估計的影響,以評估模型的穩(wěn)定性。

4.結(jié)果的輸出與解釋

上述模型建立后,可通過summary( )語句查看模型結(jié)果,以分布滯后模型為例,其模型配合結(jié)果如表2。

> summary(fit3) #查看模型配合結(jié)果

表2 分布滯后模型的配合結(jié)果

*：P<0.05,**:P<0.001。

從表2的模型配合結(jié)果可以看出,大氣PM10濃度(Lag01時)的回歸系數(shù)β=0.0007490,標準誤SE=0.0001278,結(jié)果具有統(tǒng)計學(xué)意義(P<0.001)。上述結(jié)果也可通過R程序直接提取。

> b<-summary(fit3) $ coeff[2,1] #提取PM10效應(yīng)的回歸系數(shù)

> se<-summary(fit3) $ coeff[2,2] #提取PM10效應(yīng)回歸系數(shù)的標準誤

> ER<-(exp(b*10)-1)*100 #計算PM10每升高10μg/m3,死亡的超額危險度ER

> ERlower<-(exp((b-1.96*se)*10)-1)*100 #計算ER的95%CI下限

>ERupper<-(exp((b+1.96*se)*10)-1)*100 #計算ER的95%CI上限

上述命令通過提取fit3模型配合結(jié)果中的回歸系數(shù)與標準誤,可計算出PM10濃度每升高10μg/m3,人群總死亡增加的百分比及其95%CI,即為0.75%(95%CI:0.50%～1.00%)。fit1模型、fit2模型結(jié)果輸出與此類似,此處不再贅述。

5.其他常用命令

除進行回歸模型的擬合外,分析過程中還經(jīng)常需要繪制各變量的時間序列圖、殘差自相關(guān)圖與偏自相關(guān)圖以及暴露-反應(yīng)關(guān)系圖等。下面列舉了一些常用的基本命令:

> plot(data $ death) #繪制時間序列圖,以死亡人數(shù)為例

> acf(resid(fit1)) #繪制殘差自相關(guān)圖,以fit1模型為例

>pacf(resid(fit1)) #繪制偏自相關(guān)圖,以fit1模型為例

>model<-gam(death～s(pm10lag01,k=4,fx=T)+ ns(temp,3) +ns(rhum,3) + ns(time,7*14)+ dow,family=quasipoisson(),data)

>plot (model) #繪制模型model條件下的暴露-反應(yīng)關(guān)系圖

小 結(jié)

近年來,R軟件在空氣污染與健康領(lǐng)域的統(tǒng)計軟件包已經(jīng)比較成熟。對于使用者而言,只要具備了相應(yīng)的方法學(xué)基礎(chǔ),即可方便地進行數(shù)據(jù)分析與圖形制作。本文采用廣義線性模型、廣義相加模型結(jié)合分布滯后模型,利用逐日大氣污染物、氣象因素及死亡人數(shù)的時間序列數(shù)據(jù),簡要介紹了時間序列分析在空氣污染與健康領(lǐng)域的應(yīng)用,并給出了R軟件實現(xiàn)的詳細過程,可為我國空氣污染健康影響研究工作提供實用的操作方法參考。但時間序列分析并不是一個固定的模型,模型結(jié)構(gòu)是由研究時所采用數(shù)據(jù)的內(nèi)在聯(lián)系決定的,建議使用者結(jié)合自己的數(shù)據(jù)特點選擇合適的模型和參數(shù)靈活分析。