焦奎壯，馬煦晰，馬小茜，劉朝屹，張 青，馬 露

武漢大學(xué)公共衛(wèi)生學(xué)院（武漢 430071）

縱向數(shù)據(jù)是流行病學(xué)研究中最常見的資料類型之一，常見于隊(duì)列研究、定群研究等研究設(shè)計(jì)中。這些研究設(shè)計(jì)中暴露因素或健康結(jié)局指標(biāo)均具有隨時(shí)間變化而變化的特點(diǎn)，為探索暴露因素與健康結(jié)局間的關(guān)系，往往需對(duì)研究對(duì)象進(jìn)行隨訪或重復(fù)測(cè)量。針對(duì)同一研究對(duì)象的多次測(cè)量，研究結(jié)果間通常存在相關(guān)關(guān)系，若不滿足獨(dú)立性的條件，不適用一般線性模型或廣義線性模型，且當(dāng)觀測(cè)值存在缺失時(shí)，重復(fù)測(cè)量方差分析也不適用。因此，廣義估計(jì)方程（generalized estimating equations，GEE）和混合線性模型（mixed linear model，MLM）被廣泛應(yīng)用于縱向數(shù)據(jù)的統(tǒng)計(jì)分析。目前，主流的統(tǒng)計(jì)分析軟件如SAS、SPSS和R等均能實(shí)現(xiàn)GEE和MLM的建模分析[1-3]。Python作為一款開源免費(fèi)軟件，因其強(qiáng)大的大數(shù)據(jù)處理第三方庫(kù)（Pandas、Numpy、Scipy等）、內(nèi)存優(yōu)化系統(tǒng)和豐富的應(yīng)用場(chǎng)景（爬蟲等），可方便快速地實(shí)現(xiàn)數(shù)據(jù)的獲取、清洗、管理和分析，顯著縮短數(shù)據(jù)分析時(shí)間[4-5]，近年來(lái)越來(lái)越受到國(guó)內(nèi)科研工作者的歡迎。隨著我國(guó)醫(yī)療系統(tǒng)信息化建設(shè)的快速推進(jìn)，醫(yī)療大數(shù)據(jù)的智能化統(tǒng)計(jì)分析是必然的發(fā)展趨勢(shì)[6]。運(yùn)用Python軟件實(shí)現(xiàn)流行病學(xué)研究中的統(tǒng)計(jì)分析，目前尚不多見。本研究以Python 3.8.5中的statsmodels庫(kù)為例，通過(guò)研究實(shí)例介紹GEE和MLM在Python軟件中的實(shí)現(xiàn)方法。同時(shí)，采用R 4.0.5軟件中的geepack包和lmerTest包構(gòu)建GEE與MLM模型[7]，作為本次Python結(jié)果的對(duì)照，驗(yàn)證Python輸出結(jié)果是否正確。

1 資料與方法

1.1 資料來(lái)源

為研究某地區(qū)大氣顆粒物PM2.5對(duì)肺功能的影響，收集該地區(qū)連續(xù)十日的日均PM2.5濃度、溫濕度和研究對(duì)象的肺功能數(shù)據(jù)。其中，PM2.5濃度、溫濕度數(shù)據(jù)分別來(lái)自武漢市生態(tài)環(huán)境保護(hù)局（http://hbj.wuhan.gov.cn/）和湖北省氣象局（http://hb.cma.gov.cn/），均采用Python爬蟲收集。

采用方便抽樣方法抽取研究對(duì)象，使用第1秒用力呼氣容積（forced expiratory volume in one second，F(xiàn)EV1）作為肺功能評(píng)估指標(biāo)，使用肺功能測(cè)試儀于每日7:30 AM由研究對(duì)象自測(cè)，測(cè)量10日后通過(guò)軟件將數(shù)據(jù)讀入計(jì)算機(jī)。以某位研究對(duì)象（20歲）為例，顯示十天連續(xù)測(cè)量數(shù)據(jù)，詳見表1。

表1 連續(xù)10日顆粒物、氣象及研究對(duì)象數(shù)據(jù)資料Table 1. Data of particulate matter, meteorology and research objects in ten consecutive days

1.2 模型構(gòu)建

本研究通過(guò)控制研究對(duì)象年齡、BMI和大氣溫濕度，評(píng)價(jià)PM2.5對(duì)研究對(duì)象肺功能的影響，以PM2.5單污染物滯后2天的暴露模型為例展示GEE與MLM在Python軟件中的實(shí)現(xiàn)。

1.2.1 GEE模型構(gòu)建

GEE模型是Liang 和 Zeger于1986年在擬似然的基礎(chǔ)上對(duì)廣義線性模型的推廣，旨在分析縱向數(shù)據(jù)中因素對(duì)總體平均水平的影響[8-10]。本研究以FEV1為因變量，記為Yij，表示有i個(gè)研究對(duì)象（1，…，n），每個(gè)研究對(duì)象有j個(gè)觀察值（1，…，p），協(xié)變量記為Xij。構(gòu)建如下模型：

建立Yij與Xij間的函數(shù)關(guān)系，鏈接函數(shù)的選擇主要有以下幾種形式：①因變量服從正態(tài)分布，鏈接函數(shù)選擇identify；②因變量服從二項(xiàng)分布，鏈接函數(shù)選擇logit；③因變量服從泊松分布，鏈接函數(shù)選擇log；④因變量服從負(fù)二項(xiàng)分布，鏈接函數(shù)選擇negativebinomial。

（2）Var(Yij) = v(μij)·φ，建立 Yij的方差與均值間的函數(shù)關(guān)系。

（3）Ri(α)是Yij的作業(yè)相關(guān)矩陣，表示因變量各次的重復(fù)測(cè)量值間的相關(guān)性大小，作業(yè)相關(guān)矩陣包括以下幾種形式：①可交換，又稱等相關(guān)，即任意兩次不同時(shí)間的觀測(cè)值間相關(guān)是相等的；②相鄰相關(guān)，即只有相鄰時(shí)間的兩次測(cè)量值之間具有相關(guān)性；③自相關(guān)，即相關(guān)與間隔次數(shù)有關(guān)，相隔次數(shù)越長(zhǎng)，相關(guān)關(guān)系越??；④不確定型相關(guān)，即相關(guān)矩陣非對(duì)角線上的元素均不等；⑤獨(dú)立，即因變量之間不相關(guān)。

根據(jù)Liang和Zeger的定義，構(gòu)建如下 GEE模型：

Vi是Yi的協(xié)方差矩陣，。根據(jù)給定的α和φ的估計(jì)值，通過(guò)迭代重復(fù)加權(quán)，采用最小二乘法，求解上述方程，最后得出β的一致性估計(jì)。

1.2.2 MLM模型構(gòu)建

MLM是一般線性模型（Y = Xβ + ε）的拓展[11-12]，其保留了傳統(tǒng)線性模型中的殘差需滿足正態(tài)性的假定，而對(duì)獨(dú)立性與方差齊性不做要求[13]，引入隨機(jī)效應(yīng)部分 Zγ，可表達(dá)為：Y = Xβ + Zγ + ε，Y表示因變量測(cè)量值的向量，X為固定效應(yīng)設(shè)計(jì)矩陣，β為固定效應(yīng)參數(shù)向量，Z為隨機(jī)效應(yīng)設(shè)計(jì)矩陣，γ為隨機(jī)效應(yīng)參數(shù)向量，服從均值向量為0、方差/協(xié)方差矩陣為G的正態(tài)分布，表示為γ～N(0，G)，隨機(jī)效應(yīng)主要有以下三種[14-15]：①隨機(jī)回歸系數(shù)帶來(lái)的隨機(jī)效應(yīng)；②隨機(jī)截距帶來(lái)的隨機(jī)效應(yīng)；③隨機(jī)回歸系數(shù)和隨機(jī)截距帶來(lái)的隨機(jī)效應(yīng)。ε為隨機(jī)誤差向量，服從均值向量為0、方差/協(xié)方差矩陣為R的正態(tài)分布，即ε～N(0，R)。

1.3 模型驗(yàn)證

采用R4.0.5軟件的geepack包、lmerTest包分別構(gòu)建GEE與MLM模型，作為參照，驗(yàn)證Python輸出結(jié)果是否正確。

2 結(jié)果

2.1 GEE建模及參數(shù)估計(jì)結(jié)果

本研究因變量FEV1為定量連續(xù)數(shù)據(jù)，服從正態(tài)分布，故鏈接函數(shù)選擇identify，根據(jù)以上定義，可構(gòu)建如下GEE模型：

β0為截距，β1、β2、β3、β4、β5為各協(xié)變量的回歸系數(shù)，CORR表示作業(yè)相關(guān)矩陣，ε表示殘差項(xiàng)，相關(guān)代碼見框1，GEE參數(shù)估計(jì)結(jié)果如表2所示。

框1 廣義估計(jì)方程代碼呈現(xiàn)Box 1. GEE codes in Python and R

表2 廣義估計(jì)方程參數(shù)估計(jì)結(jié)果Table 2. Results of GEE parameter estimation

2.2 MLM建模及參數(shù)估計(jì)結(jié)果

根據(jù)MLM定義，可構(gòu)建如下MLM模型：FEV1 = β0+ β1PM2.5+ β2Age + β3BMI +β4Temperature + β5Humidity + Ζγ + ε

β0為總截距，β1、β2、β3、β4、β5為各協(xié)變量的回歸系數(shù)，ε為殘差項(xiàng)，Zγ是隨機(jī)效應(yīng)，對(duì)應(yīng)研究對(duì)象的個(gè)體差異帶來(lái)的隨機(jī)截距。相關(guān)代碼見框2，MLM參數(shù)估計(jì)結(jié)果如表3所示。

表3 混合線性模型參數(shù)估計(jì)結(jié)果Table 3. Results of MLM parameter estimation

根據(jù)模型中自變量的檢驗(yàn)結(jié)果可知，PM2.5是肺功能降低的危險(xiǎn)因素，PM2.5每升高1 μg/m3，研究對(duì)象2天后的FEV1減少8 mL（P＜0.05）。GEE模型中，Python與R輸出的檢驗(yàn)統(tǒng)計(jì)量分別為z和χ2，同一自變量參數(shù)檢驗(yàn)結(jié)果與P值基本一致；MLM模型中，Python與R輸出的檢驗(yàn)統(tǒng)計(jì)量分別為z和t，同一自變量參數(shù)檢驗(yàn)結(jié)果與P值基本一致，表示Python輸出結(jié)果基本可信。

框2 混合線性模型代碼呈現(xiàn)Box 2. MLM codes in Python and R

3 討論

本文結(jié)合流行病學(xué)中的研究實(shí)例，簡(jiǎn)要介紹了GEE和MLM在Python中的具體操作，拓展了縱向數(shù)據(jù)分析的實(shí)現(xiàn)方法。盡管縱向數(shù)據(jù)在Python與R軟件中的實(shí)現(xiàn)代碼有些許不同，但輸出的參數(shù)和P值檢驗(yàn)結(jié)果基本一致，反映Python輸出結(jié)果可信。在不同的軟件中，相同的方法模型輸出的統(tǒng)計(jì)量有所差異，這是開發(fā)人員的統(tǒng)計(jì)檢驗(yàn)傾向?qū)е碌?。鑒于Z統(tǒng)計(jì)量值的平方等于Wald χ2統(tǒng)計(jì)量的值，因此Wald χ2檢驗(yàn)是等價(jià)于Z檢驗(yàn)的[16]；在大樣本（n＞50）的情況下，t檢驗(yàn)與Z檢驗(yàn)結(jié)果也是近似的[17]，所以軟件輸出參數(shù)統(tǒng)計(jì)量的不同不影響參數(shù)檢驗(yàn)結(jié)果。

Python在數(shù)據(jù)分析方面存在諸多優(yōu)點(diǎn)，其完善的作圖功能以及豐富的數(shù)據(jù)分析庫(kù)、機(jī)器學(xué)習(xí)庫(kù)的發(fā)展，越來(lái)越符合大數(shù)據(jù)背景下的數(shù)據(jù)分析要求。例如，陳偉等使用Python相似度分析與標(biāo)簽云分析方法進(jìn)行文本數(shù)據(jù)分析，拓展了大數(shù)據(jù)審計(jì)的研究方向[18]；楊俊秀等采用Python的matplotlib庫(kù)、numpy庫(kù)和scikit-learn庫(kù)完成了數(shù)據(jù)的可視化、整理與分析，顯著提高了高頻電子線路實(shí)驗(yàn)數(shù)據(jù)的處理分析效率[5]。但Python也存在著不足之處，其用作數(shù)據(jù)分析工具的時(shí)間較短，更多的數(shù)據(jù)分析方法程序有待進(jìn)一步開發(fā)及完善。本研究使用的statsmodels庫(kù)目前可實(shí)現(xiàn)簡(jiǎn)單線性模型、廣義線性模型、GEE、MLM等大部分模型的建模。隨著statsmodels庫(kù)的進(jìn)一步開發(fā)，相信未來(lái)statsmodels庫(kù)可滿足更多的建模需求。

綜上，將Python用于環(huán)境流行病學(xué)，在實(shí)現(xiàn)數(shù)據(jù)獲取、處理與分析的過(guò)程中，統(tǒng)一了語(yǔ)言環(huán)境，避免了數(shù)據(jù)在不同平臺(tái)間的轉(zhuǎn)換，提高了數(shù)據(jù)分析的效率。使用Python軟件可靈活實(shí)現(xiàn) GEE和MLM的統(tǒng)計(jì)建模，在實(shí)際研究中有一定參考價(jià)值。