武漢大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)學(xué)系(430071)
張?jiān)茩?quán) 朱耀輝 李存祿 馮仁杰 馬 露△
廣義相加模型在R軟件中的實(shí)現(xiàn)
武漢大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)學(xué)系(430071)
張?jiān)茩?quán) 朱耀輝 李存祿 馮仁杰 馬 露△
目的 通過(guò)R軟件實(shí)現(xiàn)廣義相加模型。方法 通過(guò)空間污染流行病學(xué)的一個(gè)實(shí)例研究介紹利用R軟件mgcv包實(shí)現(xiàn)廣義相關(guān)模型的具體步驟和評(píng)價(jià)方法。結(jié)果 廣義相加模型可在R軟件中方便實(shí)現(xiàn)。結(jié)論 R軟件作為一款自由、免費(fèi)、開(kāi)源的統(tǒng)計(jì)分析軟件,可靈活方便地構(gòu)建廣義相加模型,在實(shí)際研究中值得推廣。
廣義相加模型 空氣污染流行病學(xué) R軟件
廣義相加模型(generalized additional model,GAM)是對(duì)傳統(tǒng)廣義線性模型的非參數(shù)拓展,可有效處理解釋變量與效應(yīng)變量間復(fù)雜的非線性關(guān)系[1]。GAM目前已廣泛應(yīng)用于空氣污染流行病學(xué)研究中,主要用于分析空氣污染或氣象因素對(duì)人群健康事件(如發(fā)病、住院和死亡)的急性損害效應(yīng)。目前,國(guó)內(nèi)學(xué)者構(gòu)建GAM模型主要采用SAS軟件中的PROC GAM模塊實(shí)現(xiàn),但由于SAS軟件價(jià)格昂貴,大大阻礙了GAM模型的應(yīng)用。R作為一款自由、免費(fèi)、開(kāi)源的統(tǒng)計(jì)分析軟件,近年來(lái)已逐漸受到越來(lái)越多的科研工作者的重視和青睞。R軟件自帶默認(rèn)包中的glm函數(shù)以及mgcv包的gam函數(shù),均可用于構(gòu)建GAM模型[2],在國(guó)外空氣污染流行病學(xué)研究中已得到廣泛應(yīng)用。本文以R3.1.1中的mgcv包為例,通過(guò)一個(gè)研究實(shí)例簡(jiǎn)要介紹GAM在R軟件中的實(shí)現(xiàn)方法。
1.數(shù)據(jù)資料
為研究某地區(qū)大氣污染物對(duì)居民呼吸系統(tǒng)疾病入院人次的影響,研究人員收集了該地區(qū)2009年1月1日至2010年12月31日的日均大氣污染物濃度(PM10、SO2和NO2)、日均濕度、相對(duì)濕度以及每日的呼吸系統(tǒng)疾病入院人數(shù)等資料。具體數(shù)據(jù)形式見(jiàn)表1(僅顯示部分?jǐn)?shù)據(jù))。
時(shí)間序列(time)是2009-2010年每日呼吸系統(tǒng)疾病入院這一事件發(fā)生次序的一列數(shù),故取值為1,2,3,…,730;星期變量(dow)用于控制短期波動(dòng);本研究據(jù)以往文獻(xiàn)報(bào)道簡(jiǎn)化為二分類(lèi),dow=0表示工作日,dow=1則表示周末(周六和周日)。
2.R實(shí)現(xiàn)
(1)構(gòu)建基礎(chǔ)模型
本研究旨在控制時(shí)間長(zhǎng)期趨勢(shì)和季節(jié)趨勢(shì)、“星期幾”效應(yīng)、氣象因素等混雜影響的基礎(chǔ)上,評(píng)價(jià)大氣污染物對(duì)人群呼吸系統(tǒng)疾病入院率的影響。因而,結(jié)合本研究實(shí)際數(shù)據(jù),構(gòu)建如下基礎(chǔ)模型:
Yt~Poisson(μt)
Log(μt)=s(time,dft)+as.factor(dow)+α
其中,Yt為第t日實(shí)際入院人次(服從Poisson分布);μt為第t日入院人次的期望值;s表示非參數(shù)平滑函數(shù);dft為非參數(shù)平滑函數(shù)中控制時(shí)間長(zhǎng)期趨勢(shì)和季節(jié)趨勢(shì)的自由度。以下為構(gòu)建基礎(chǔ)模型的R軟件代碼(其中df待確定,#后為注釋語(yǔ)句,不在程序中運(yùn)行):
install.packages(“mgcv”)#安裝mgcv包
library(mgcv)#載入mgcv包
base_mod<-gam(y ~ s(time,df) +as.factor(dow),family = poisson,data = mydata)
#利用mgcv包中g(shù)am函數(shù)建立基礎(chǔ)模型base_mod,指定分布族為Poisson分布,數(shù)據(jù)集為mydata
(2)確定模型自由度
在基礎(chǔ)模型構(gòu)建之后,最重要的工作就是確定模型中非參數(shù)平滑函數(shù)的自由度df。在廣義相加模型中,由于平滑函數(shù)的自由度對(duì)模型的參數(shù)估計(jì)和模型穩(wěn)定性有一定影響。因而,選擇合適的自由度對(duì)模型構(gòu)建有重大意義,通常根據(jù)以下評(píng)判準(zhǔn)則[2]進(jìn)行設(shè)定:
①基于生物學(xué)知識(shí)和專(zhuān)家經(jīng)驗(yàn)(包括敏感性分析)設(shè)置固定的自由度;
②赤池信息準(zhǔn)則(Akaike information criterion,AIC),依據(jù)AIC最小選擇自由度。
③依據(jù)殘差獨(dú)立原則,通過(guò)最小化模型殘差自相關(guān)來(lái)選擇自由度。實(shí)際工作中,我們根據(jù)基礎(chǔ)模型殘差的偏自相關(guān)(PACF)絕對(duì)值之和最小選取自由度。
④依據(jù)廣義交叉驗(yàn)證(generalized cross-validation,GCV)預(yù)測(cè)污染物濃度的最佳模型(GCV-PM10)選擇自由度,這種方法是最小化誤差均方過(guò)程的一種簡(jiǎn)化。
本研究中,時(shí)間的自由度則通過(guò)最小化模型殘差自相關(guān)[4]來(lái)選擇。在R軟件中,筆者通過(guò)編寫(xiě)循環(huán)語(yǔ)句,設(shè)定時(shí)間的每年自由度為i(從1到20),分別構(gòu)建20個(gè)相應(yīng)自由度的模型并計(jì)算出每個(gè)模型殘差偏自相關(guān)絕對(duì)值的和stat,并利用plot作圖方式將i與stat之間的關(guān)系顯示出來(lái)(亦可通過(guò)逐一建立不同自由度的多個(gè)GAM模型,對(duì)每個(gè)模型其殘差偏自相關(guān)絕對(duì)值的和進(jìn)行比較,從而最終確定模型自由度)。以下為作者自行編寫(xiě)的R循環(huán)語(yǔ)句代碼(僅供參考):
stat<-NULL
mod<-list()
for(i in 1:20){
mod[[i]]<-gam(y~s(time,df=2*i)+as.factor(dow),family=poisson,data=mydata)
tt<-sum(abs((pacf(mod[[i]]$residuals))$acf))
stat<-append(stat,tt)
}
上述代碼中,列表對(duì)象mod中包含20個(gè)GAM模型,tt為每個(gè)模型殘差偏自相關(guān)絕對(duì)值的和,stat則是由20個(gè)GAM模型的tt值構(gòu)成的向量。
根據(jù)圖1并結(jié)合最小化模型殘差自相關(guān)原則可知,本研究中時(shí)間的非參數(shù)平滑自由度應(yīng)設(shè)定為28(14/年)。同時(shí),考慮到溫度和濕度等氣象因素也可能與人群健康有關(guān),根據(jù)既往專(zhuān)家經(jīng)驗(yàn)設(shè)定溫度和相對(duì)濕度的非參數(shù)平滑函數(shù)自由度為3[3],因而本研究模型調(diào)整為:
Log(μt)=s(time,df=14×2)+as.factor(dow)+s(temperature,df=3)+s(humidity,df=3)+α
(3)構(gòu)建污染物模型
本研究從污染物的眾多暴露模型(包括滯后和累計(jì)平均)和污染物模型(單污染物和多污染物)中,以SO2單污染物滯后3d的暴露模型為例,探討SO2對(duì)人群呼吸系統(tǒng)疾病入院率的影響,模型為:
Log(μt)=βSO2+s(time,df=14×2)+as.factor(dow)+s(temperature,df=3)+s(humidity,df=3)+α
相應(yīng)R程序?yàn)椋?/p>
final_mod<-gam(y~so2+s(time,df=14*2)+s(temperature,df=3)+s(humidity,df=3)+as.factor(dow),family=poisson,data=mydata)
summary(final_mod)
模型建立后,可通過(guò)summary(final_mod)語(yǔ)句查看模型結(jié)果,其中參數(shù)估計(jì)結(jié)果見(jiàn)表2所示。
(4)模型的比較與評(píng)價(jià)
在R中,GAM的嵌套模型可通過(guò)anova(model1,model2,test=“Chisq”)語(yǔ)句進(jìn)行比較,也可通過(guò)AIC(model1,model2)直接比較模型的AIC值。通過(guò)比較選出最優(yōu)模型后,可通過(guò)觀察模型中非參數(shù)平滑函數(shù)的自由度改變對(duì)污染物參數(shù)的影響大小來(lái)最終評(píng)價(jià)模型是否穩(wěn)健,亦稱(chēng)敏感性分析。
(5)模型參數(shù)的解釋
由于本研究實(shí)例中,分布族為Poisson分布,鏈接函數(shù)為對(duì)數(shù)函數(shù),直接對(duì)模型參數(shù)進(jìn)行解釋意義不大。對(duì)模型中的系數(shù)和置信區(qū)間進(jìn)行指數(shù)轉(zhuǎn)換后得表3。
可以看出,在控制其他因素的影響后,周末入院率是工作日的1.08倍(星期變量分為周末和工作日),即周末入院率比工作日高8.29%(95%CI:2.97%~13.87%);SO2每升高1μg/m3,入院率增加0.16%(95%CI:0.04%~2.89%)。
本文結(jié)合空氣污染流行病學(xué)中的一個(gè)研究實(shí)例,簡(jiǎn)要介紹了廣義相關(guān)模型在R軟件mgcv包中的實(shí)現(xiàn)方法。由于氣溫、相對(duì)濕度、風(fēng)速等氣象因素與健康效應(yīng)之間可能存在某種非線性關(guān)聯(lián),GAM可有效處理和識(shí)別變量間的非線性相關(guān),極大地發(fā)展了傳統(tǒng)的線性模型。
R軟件中的mgcv包可通過(guò)gam函數(shù)輕松構(gòu)建GAM模型,并靈活設(shè)定模型中的各項(xiàng)參數(shù)。在GAM模型中,由于非參數(shù)平滑函數(shù)的自由度設(shè)定對(duì)模型的擬合效果至關(guān)重要,本文介紹了自由度選擇中常用的幾種方法和準(zhǔn)則。在實(shí)例研究中,基于最小化模型殘差自相關(guān)的原則,通過(guò)R中的循環(huán)語(yǔ)句對(duì)不同自由度下的模型殘差偏自相關(guān)絕對(duì)值之和進(jìn)行計(jì)算,并通過(guò)plot函數(shù)作圖對(duì)其進(jìn)行可視化展示,可非常直觀地選擇自由度。
此外,R軟件亦可通過(guò)繪制污染物以及氣象因素的暴露-反應(yīng)曲線[4-5],更直觀、形象地揭示污染物、氣象因素對(duì)健康效應(yīng)的線性或非線性影響,從而深入探索污染物的暴露閾值和氣象因素的最適宜范圍。在探索污染物與氣象因素的交互作用時(shí),亦可通過(guò)R繪制污染物、氣象因素與健康效應(yīng)的三維透視圖直觀地展示[6]。
總之,R軟件可以實(shí)現(xiàn)廣義相加模型的統(tǒng)計(jì)建模,并可靈活設(shè)定非參數(shù)平滑參數(shù),同時(shí)能對(duì)污染物和氣象因素的暴露-反應(yīng)關(guān)系以及交互效應(yīng)進(jìn)行可視化探索。
[1]英董,趙耐青,湯軍克,等.廣義相加模型在氣溫健康效應(yīng)研究中的應(yīng)用.中國(guó)衛(wèi)生統(tǒng)計(jì),2008,25(2):144-146.
[2]Peng RD,Dominici F,Louis TA.Model choice in time series studies of air pollution and mortality .J R Stat Soc Ser A Stat Soc,2006,169(2):179-198.
[3]張衍燊,周脈耕,賈予平,等.天津市可吸入顆粒物與城區(qū)居民每日死亡關(guān)系的時(shí)間序列分析.中華流行病學(xué)雜志,2010,31(5):544-548.
[4]楊敏娟,潘小川.北京市大氣污染與居民心腦血管疾病死亡的時(shí)間序列分析.環(huán)境與健康雜志,2008,25(4):294-297.
[5]張?jiān)?闞海東,彭麗,等.日均氣溫與呼吸系統(tǒng)疾病日入院人次相關(guān)性的時(shí)間序列分析.中華預(yù)防醫(yī)學(xué)雜志,2014,48(9):795-799.
[6]Burkart K,Canario P,Breitner S,et al.Interactive short-term effects of equivalent temperature and air pollution on human mortality in Berlin and Lisbon.Environ Pollut,2013,183:54-63.
(責(zé)任編輯:郭海強(qiáng))
△通信作者:馬露,E-mail:malu@whu.edu.cn