武漢大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計學(xué)系(430071)

張云權(quán) 朱耀輝 李存祿 馮仁杰 馬 露△

廣義相加模型在R軟件中的實現(xiàn)

武漢大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計學(xué)系(430071)

張云權(quán) 朱耀輝 李存祿 馮仁杰 馬 露△

目的 通過R軟件實現(xiàn)廣義相加模型。方法 通過空間污染流行病學(xué)的一個實例研究介紹利用R軟件mgcv包實現(xiàn)廣義相關(guān)模型的具體步驟和評價方法。結(jié)果 廣義相加模型可在R軟件中方便實現(xiàn)。結(jié)論 R軟件作為一款自由、免費、開源的統(tǒng)計分析軟件，可靈活方便地構(gòu)建廣義相加模型，在實際研究中值得推廣。

廣義相加模型 空氣污染流行病學(xué) R軟件

廣義相加模型(generalized additional model，GAM)是對傳統(tǒng)廣義線性模型的非參數(shù)拓展，可有效處理解釋變量與效應(yīng)變量間復(fù)雜的非線性關(guān)系[1]。GAM目前已廣泛應(yīng)用于空氣污染流行病學(xué)研究中，主要用于分析空氣污染或氣象因素對人群健康事件(如發(fā)病、住院和死亡)的急性損害效應(yīng)。目前，國內(nèi)學(xué)者構(gòu)建GAM模型主要采用SAS軟件中的PROC GAM模塊實現(xiàn)，但由于SAS軟件價格昂貴，大大阻礙了GAM模型的應(yīng)用。R作為一款自由、免費、開源的統(tǒng)計分析軟件，近年來已逐漸受到越來越多的科研工作者的重視和青睞。R軟件自帶默認包中的glm函數(shù)以及mgcv包的gam函數(shù)，均可用于構(gòu)建GAM模型[2]，在國外空氣污染流行病學(xué)研究中已得到廣泛應(yīng)用。本文以R3.1.1中的mgcv包為例，通過一個研究實例簡要介紹GAM在R軟件中的實現(xiàn)方法。

研究實例

1.數(shù)據(jù)資料

為研究某地區(qū)大氣污染物對居民呼吸系統(tǒng)疾病入院人次的影響，研究人員收集了該地區(qū)2009年1月1日至2010年12月31日的日均大氣污染物濃度(PM10、SO2和NO2)、日均濕度、相對濕度以及每日的呼吸系統(tǒng)疾病入院人數(shù)等資料。具體數(shù)據(jù)形式見表1(僅顯示部分數(shù)據(jù))。

時間序列(time)是2009-2010年每日呼吸系統(tǒng)疾病入院這一事件發(fā)生次序的一列數(shù)，故取值為1，2，3，…，730；星期變量(dow)用于控制短期波動；本研究據(jù)以往文獻報道簡化為二分類，dow=0表示工作日，dow=1則表示周末(周六和周日)。

2.R實現(xiàn)

(1)構(gòu)建基礎(chǔ)模型

本研究旨在控制時間長期趨勢和季節(jié)趨勢、“星期幾”效應(yīng)、氣象因素等混雜影響的基礎(chǔ)上，評價大氣污染物對人群呼吸系統(tǒng)疾病入院率的影響。因而，結(jié)合本研究實際數(shù)據(jù)，構(gòu)建如下基礎(chǔ)模型：

Yt～Poisson(μt)

Log(μt)=s(time,dft)+as.factor(dow)+α

其中，Yt為第t日實際入院人次(服從Poisson分布)；μt為第t日入院人次的期望值；s表示非參數(shù)平滑函數(shù)；dft為非參數(shù)平滑函數(shù)中控制時間長期趨勢和季節(jié)趨勢的自由度。以下為構(gòu)建基礎(chǔ)模型的R軟件代碼(其中df待確定，#后為注釋語句，不在程序中運行)：

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ù)的自由度對模型的參數(shù)估計和模型穩(wěn)定性有一定影響。因而，選擇合適的自由度對模型構(gòu)建有重大意義，通常根據(jù)以下評判準則[2]進行設(shè)定：

①基于生物學(xué)知識和專家經(jīng)驗(包括敏感性分析)設(shè)置固定的自由度；

②赤池信息準則(Akaike information criterion，AIC)，依據(jù)AIC最小選擇自由度。

③依據(jù)殘差獨立原則，通過最小化模型殘差自相關(guān)來選擇自由度。實際工作中，我們根據(jù)基礎(chǔ)模型殘差的偏自相關(guān)(PACF)絕對值之和最小選取自由度。

④依據(jù)廣義交叉驗證(generalized cross-validation，GCV)預(yù)測污染物濃度的最佳模型(GCV-PM10)選擇自由度，這種方法是最小化誤差均方過程的一種簡化。

本研究中，時間的自由度則通過最小化模型殘差自相關(guān)[4]來選擇。在R軟件中，筆者通過編寫循環(huán)語句，設(shè)定時間的每年自由度為i(從1到20)，分別構(gòu)建20個相應(yīng)自由度的模型并計算出每個模型殘差偏自相關(guān)絕對值的和stat，并利用plot作圖方式將i與stat之間的關(guān)系顯示出來(亦可通過逐一建立不同自由度的多個GAM模型，對每個模型其殘差偏自相關(guān)絕對值的和進行比較，從而最終確定模型自由度)。以下為作者自行編寫的R循環(huán)語句代碼(僅供參考)：

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)

}

上述代碼中，列表對象mod中包含20個GAM模型，tt為每個模型殘差偏自相關(guān)絕對值的和，stat則是由20個GAM模型的tt值構(gòu)成的向量。

根據(jù)圖1并結(jié)合最小化模型殘差自相關(guān)原則可知，本研究中時間的非參數(shù)平滑自由度應(yīng)設(shè)定為28(14/年)。同時，考慮到溫度和濕度等氣象因素也可能與人群健康有關(guān)，根據(jù)既往專家經(jīng)驗設(shè)定溫度和相對濕度的非參數(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)建污染物模型

本研究從污染物的眾多暴露模型(包括滯后和累計平均)和污染物模型(單污染物和多污染物)中，以SO2單污染物滯后3d的暴露模型為例，探討SO2對人群呼吸系統(tǒng)疾病入院率的影響，模型為：

Log(μt)=βSO2+s(time,df=14×2)+as.factor(dow)+s(temperature,df=3)+s(humidity,df=3)+α

相應(yīng)R程序為：

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)

模型建立后，可通過summary(final_mod)語句查看模型結(jié)果，其中參數(shù)估計結(jié)果見表2所示。

(4)模型的比較與評價

在R中，GAM的嵌套模型可通過anova(model1,model2,test=“Chisq”)語句進行比較，也可通過AIC(model1,model2)直接比較模型的AIC值。通過比較選出最優(yōu)模型后，可通過觀察模型中非參數(shù)平滑函數(shù)的自由度改變對污染物參數(shù)的影響大小來最終評價模型是否穩(wěn)健，亦稱敏感性分析。

(5)模型參數(shù)的解釋

由于本研究實例中，分布族為Poisson分布，鏈接函數(shù)為對數(shù)函數(shù)，直接對模型參數(shù)進行解釋意義不大。對模型中的系數(shù)和置信區(qū)間進行指數(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é)中的一個研究實例，簡要介紹了廣義相關(guān)模型在R軟件mgcv包中的實現(xiàn)方法。由于氣溫、相對濕度、風速等氣象因素與健康效應(yīng)之間可能存在某種非線性關(guān)聯(lián)，GAM可有效處理和識別變量間的非線性相關(guān)，極大地發(fā)展了傳統(tǒng)的線性模型。

R軟件中的mgcv包可通過gam函數(shù)輕松構(gòu)建GAM模型，并靈活設(shè)定模型中的各項參數(shù)。在GAM模型中，由于非參數(shù)平滑函數(shù)的自由度設(shè)定對模型的擬合效果至關(guān)重要，本文介紹了自由度選擇中常用的幾種方法和準則。在實例研究中，基于最小化模型殘差自相關(guān)的原則，通過R中的循環(huán)語句對不同自由度下的模型殘差偏自相關(guān)絕對值之和進行計算，并通過plot函數(shù)作圖對其進行可視化展示，可非常直觀地選擇自由度。

此外，R軟件亦可通過繪制污染物以及氣象因素的暴露-反應(yīng)曲線[4-5]，更直觀、形象地揭示污染物、氣象因素對健康效應(yīng)的線性或非線性影響，從而深入探索污染物的暴露閾值和氣象因素的最適宜范圍。在探索污染物與氣象因素的交互作用時，亦可通過R繪制污染物、氣象因素與健康效應(yīng)的三維透視圖直觀地展示[6]。

總之，R軟件可以實現(xiàn)廣義相加模型的統(tǒng)計建模，并可靈活設(shè)定非參數(shù)平滑參數(shù)，同時能對污染物和氣象因素的暴露-反應(yīng)關(guān)系以及交互效應(yīng)進行可視化探索。

[1]英董,趙耐青,湯軍克,等.廣義相加模型在氣溫健康效應(yīng)研究中的應(yīng)用.中國衛(wèi)生統(tǒng)計,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)系的時間序列分析.中華流行病學(xué)雜志,2010,31(5):544-548.

[4]楊敏娟,潘小川.北京市大氣污染與居民心腦血管疾病死亡的時間序列分析.環(huán)境與健康雜志,2008,25(4):294-297.

[5]張越,闞海東,彭麗,等.日均氣溫與呼吸系統(tǒng)疾病日入院人次相關(guān)性的時間序列分析.中華預(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.

(責任編輯：郭海強)

△通信作者：馬露，E-mail：malu@whu.edu.cn