廣義相加模型在R軟件中的實現

武漢大學公共衛生學院流行病與衛生統計學系(430071)

張云權 朱耀輝 李存祿 馮仁杰 馬 露△

廣義相加模型在R軟件中的實現

武漢大學公共衛生學院流行病與衛生統計學系(430071)

張云權 朱耀輝 李存祿 馮仁杰 馬 露△

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

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

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

研究實例

1.數據資料

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

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

2.R實現

(1)構建基礎模型

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

Yt～Poisson(μt)

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

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

install.packages(“mgcv”)#安裝mgcv包

library(mgcv)#載入mgcv包

base_mod<-gam(y ～ s(time,df) +as.factor(dow),family = poisson,data = mydata)

#利用mgcv包中gam函數建立基礎模型base_mod，指定分布族為Poisson分布，數據集為mydata

(2)確定模型自由度

在基礎模型構建之后，最重要的工作就是確定模型中非參數平滑函數的自由度df。在廣義相加模型中，由于平滑函數的自由度對模型的參數估計和模型穩定性有一定影響。因而，選擇合適的自由度對模型構建有重大意義，通常根據以下評判準則[2]進行設定：

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

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

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

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

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

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為每個模型殘差偏自相關絕對值的和，stat則是由20個GAM模型的tt值構成的向量。

根據圖1并結合最小化模型殘差自相關原則可知，本研究中時間的非參數平滑自由度應設定為28(14/年)。同時，考慮到溫度和濕度等氣象因素也可能與人群健康有關，根據既往專家經驗設定溫度和相對濕度的非參數平滑函數自由度為3[3]，因而本研究模型調整為：

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

(3)構建污染物模型

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

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

相應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)語句查看模型結果，其中參數估計結果見表2所示。

(4)模型的比較與評價

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

(5)模型參數的解釋

由于本研究實例中，分布族為Poisson分布，鏈接函數為對數函數，直接對模型參數進行解釋意義不大。對模型中的系數和置信區間進行指數轉換后得表3。

可以看出，在控制其他因素的影響后，周末入院率是工作日的1.08倍(星期變量分為周末和工作日)，即周末入院率比工作日高8.29%(95%CI:2.97%～13.87%)；SO2每升高1μg/m3，入院率增加0.16%(95%CI:0.04%～2.89%)。

討 論

本文結合空氣污染流行病學中的一個研究實例，簡要介紹了廣義相關模型在R軟件mgcv包中的實現方法。由于氣溫、相對濕度、風速等氣象因素與健康效應之間可能存在某種非線性關聯，GAM可有效處理和識別變量間的非線性相關，極大地發展了傳統的線性模型。

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

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

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

[1]英董,趙耐青,湯軍克,等.廣義相加模型在氣溫健康效應研究中的應用.中國衛生統計,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]張衍燊,周脈耕,賈予平,等.天津市可吸入顆粒物與城區居民每日死亡關系的時間序列分析.中華流行病學雜志,2010,31(5):544-548.

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

[5]張越,闞海東,彭麗,等.日均氣溫與呼吸系統疾病日入院人次相關性的時間序列分析.中華預防醫學雜志,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