武漢大學(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)方法。

研究實(shí)例

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