環境流行病學研究中廣義可加模型在Python中的實現

李湘瑩，李培政，王 靜，羅晨曦，張清煜，馬 露

武漢大學公共衛生學院（武漢 430071）

環境流行病學領域研究工作中往往涉及大量時間序列數據，其影響因素的多樣性和復雜性使得在研究時往往難以確定回歸關系的基本形式，廣義可加模型（generalized additive models，GAM）是解決這一問題的一種途徑。GAM 是在廣義線性模型（generalized linear mode，GLM）基礎上發展而來的，它提供了更為靈活的建?？蚣?，其反應變量可以為非正態的指數族分布，解釋變量與反應變量間可存在復雜的非線性關系，由指定的協變量平滑函數加上線性項的常規參數分量之和得出。伴隨環境流行病學領域研究資料的數字化高速發展，多源數據融合、智能化管理和分析已成為必然趨勢。Python 作為一個擁有強大第三方庫的開源軟件，其規范且相對簡潔的編程語言，在同應用程序整合、數據源連接與讀取、調用其他語言，以及實現人工智能等方面展現出了顯著優勢[1]。然而，在環境流行病學領域中，使用Python 進行GAM 分析的研究尚較少。本文利用環境流行病學研究中的一組時間序列資料，分別進行Python、R 和SAS 的GAM 建模，比較它們在運算邏輯、參數、結果方面的差異，以拓寬Python 軟件在環境流行病學領域的應用場景。

1 資料與方法

1.1 資料來源

因呼吸系統疾病入院人次數據主要來源于某地區衛生信息中心2017 年1 月1 日至2018 年12 月31 日記錄的該地醫院病案首頁。每日大氣中PM2.5濃度、日均溫度、平均濕度等數據來源于當地環境監測部門官方網站。PM2.5濃度和研究對象日入院人次逐日變化圖如圖1 所示。

圖1 2017—2018年某地PM2.5濃度和研究對象日入院人次逐日變化圖Figure 1. Daily variation of PM2.5 concentration and the number of hospital admissions of research objects from 2017 to 2018

1.2 模型構建

為了評價大氣污染物與因呼吸系統疾病入院人次的關系，需要對長期趨勢和季節趨勢、星期效應、平均溫度、相對濕度等產生的影響進行控制。此處長期趨勢指連續數年觀察中入院人數呈現的某種總的變動趨勢；季節趨勢指入院人數在每年隨季節產生的周期變化。這里，我們對數據庫進行了調整，為了控制長期趨勢和季節趨勢，將事件發生的日期調整為這一事件發生的先后次序，引入時間序列變量（命名為time），記作1，2，3…。為了控制星期效應的影響引入星期變量（命名為dow），工作日（星期一至星期五）入院為dow=0、周末（周六至周日）入院為dow=1。平均溫度的單位為℃，相對濕度的單位為%，PM2.5濃度單位為μg/m3，數據形式如表1（僅展示部分數據）。

表1 顆粒物、氣象及研究對象數據資料Table 1. Data of particulate matter, meteorology and research objects

本研究首先在入院前三天的單日和移動平均暴露的多個滯后結構下估計了空氣污染物對全因住院的短期影響(lag0-3, lag01-03)。然后根據模型的廣義交叉驗證（generalized cross-validation，GCV）和R2，最終選擇以PM2.5單污染物滯后0天的暴露模型作為先驗滯后結構，利用Python 中的statsmodels 庫構建GAM 模型，并與R 軟件、SAS 軟件的輸出結果進行比較，驗證Python 結果。

GAM 的基本形式為：

其中，g(μ)為連接函數，g(μ)的選擇由反應變量的分布形式決定：①當反應變量服從正態分布時，等同于一般可加模型，g(μ)=μ。②當反應變量服從二項分布時，選用logit 連接。③當反應變量服從Gamma 分布時，選用Identity 連接[2]。④當反應變量服從Poisson分布時，選用Log連接[3]。Si(·)為非參數光滑樣條函數，ε 為誤差項。模型構建過程如下：

（1）調用相關參數并載入數據

加載statsmodles 庫與pandas 庫，從statmodels庫中加載GLMGam 語句，CyclicCubicSplines 語句，以便于讀取數據和擬合模型。使用pd.read 語句載入數據保存在base_data 中。命令如下：

（2）對解釋變量進行樣條變換

在Python statsmodles 庫中，使用GLMGam 函數進行擬合。分類變量被視為線性項，對非線性解釋變量通過樣條函數進行擬合。GLMGam 函數中主要選擇的平滑基函數有兩類，分別是B 樣條函數與循環立方樣條函數。在本例中我們選用循環立方樣條函數對長期趨勢、溫度、濕度這三個變量進行樣條變換。命令如下：

（3）確定模型自由度

自由度很大程度上會影響模型擬合和預測的準確性，目前主要有四種確定模型自由度的方法：①根據既往學者的經驗設置固定的自由度；②根據最小化模型殘差自相關絕對值最小來選擇；③廣義交叉驗證（generalized crossvalidation，GCV）依據GCV 值最小選擇；④赤池信息準則（Akaike information criterion,AIC），依據AIC 值最小選擇[4]。在實際操作中對于最優自由度需要進行多次評估，本例依照既往學者的經驗并結合AIC 值、GCV 值，設置固定自由度，其中長期趨勢自由度設為15/年，平均溫度自由度設為4，相對濕度自由度設為4。命令如下：

（4）構建基礎模型

由于因呼吸疾病系統入院人次服從Poisson分布，故連接函數選擇Log 連接，根據GAM 函數定義，構建如下模型：

其中β1、β2分別為PM2.5和短期波動（dow）的系數，s為非參數平滑函數，ε 為殘差項。

使用GLMgam 函數構建模型，對三個非線性解釋變量使用循環立方樣條函數擬合，定義反應變量為泊松分布，使用fit 語句擬合函數，summary 語句查看結果。具體命令如下：

1.3 模型驗證

采 用R 3.6.1 軟 件 的mgcv 包 與SAS 9.2 中的proc gam 語句構建GAM 模型作為參照，由于Python 與SAS 軟件中可選擇樣條函數的限制，為保證結果的可比較性，本文選用R、Python 軟件中的循環立方樣條（cyclic cubic splines）函數及SAS 軟件中的三次平滑樣條函數對非參數項進行擬合，并設置統一且固定的自由度，驗證輸出結果是否正確。R、SAS 相關代碼見框1。

框1 R、SAS中的GAM代碼Box 1. GAM codes in R and SAS

2 結果

2.1 GAM建模及參數估計結果

Python、R、SAS 三種軟件輸出結果與指數轉換結果如表2。Python 與R 輸出的結果主要包括：偏回歸系數、標準誤、Z值、P值等。SAS除上述結果外還會輸出迭代匯總、平滑模型分析及圖例。

表2 Python、R、SAS參數估計結果Table 2. Results of parameter estimation in Python, R and SAS

Python 輸出的截距、PM2.5、星期變量參數估計結果分別為：ε=3.937，α=0.003，Z=1298.86，P＜0.01；β1=0.087，α=0.009，Z=9.29，P＜0.01；β2=-0.299，α=0.004，Z=-68.63，P＜0.01。

指數轉換后的結果說明，在控制了長期趨勢和季節趨勢、星期效應、平均溫度、相對濕度的影響后，每日因呼吸系統疾病入院人次與PM2.5濃度變化有關聯。三種軟件輸出的污染物及星期變量的估計值和標準誤結果有所不同，但差異較小。Python 參數估計結果顯示，PM2.5每升高100 μg/m3，入院人次增加9.1%[95%CI（7.2%，11.0%）]。R 輸出結果為PM2.5每升高100 μg/m3，入院人次增加9.1%[95%CI（7.3%，11.1%）]，SAS 輸出結果為PM2.5每升高100 μg/m3，入院人次增加9.1%[95%CI（7.6%，10.6%）]。

在模型擬合結果中，R 軟件輸出R2=0.659，Python 輸出Pseudo R2=1.000，但SAS 中僅針對各非參數項輸出GCV 值，GCV 值越小表明擬合效果越好，本例中各非參數項GCV 值均較小。

3 討論

本文結合環境流行病學研究實例，簡要介紹了GAM 模型的基本形式，并比較Python、R 和SAS 在建模和分析結果上的異同。結果發現三種軟件的輸出結果基本一致，Python 輸出結果可信。

在GAM 模型構建的參數設置上，三種軟件存在差異。第一，三者內置函數擬合過程不同。Python 軟件與R 軟件是通過懲罰似然最大化來擬合平滑參數，通過對每個平滑參數添加懲罰項來避免擬合過度或擬合不良[5]，而SAS 軟件是通過雙重迭代方法估計平滑參數[6]。第二，三者構建模型命令代碼不同。Python 軟件先統一對非線性參數項使用同一樣條函數，再將非線性項與線性項通過GLMGam 函數構建模型進行擬合。R 軟件與SAS 軟件均通過相應函數直接將線性項與非線性項放在同一語句內構建模型。第三，三者可選用的樣條函數不同。Python 軟件目前已得到驗證的樣條函數僅有B樣條函數與循環立方樣條函數，SAS 軟件中只可以使用三次平滑樣條函數、局部回歸與薄板樣條函數[7]，R 軟件則能在這兩種軟件的基礎上提供如三次自然樣條函數等更多可供選擇的平滑函數[8]。這些在參數設置和選擇上的不同可能是其結果出現差異的主要原因。

盡管在本研究中，運用Python 進行GAM 建模無論是從程序語言上還是參數設置上，都沒有展現出其在統計分析方面的優勢，但Python 作為當前主流的計算機語言之一，其在網絡爬蟲、數據可視化、機器學習等方面擁有豐富的第三方庫資源支持，能夠高效實施多源數據連接與管理，并最終使實現數據收集與處理的環境統一成為可能。不斷拓展Python 在環境流行病學領域的應用，必將為促進計算機與醫學相關領域交叉融合發展提供借鑒和參考。