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

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

病例交叉研究由美國學(xué)者麥克盧爾（McClure）于1991年提出，是一種通過比較研究對(duì)象在急性事件發(fā)生前一段時(shí)間的暴露情況與未發(fā)生事件的某段時(shí)間內(nèi)的暴露情況，來研究短暫暴露對(duì)罕見急性病的瞬間影響的流行病學(xué)方法，目前已成為環(huán)境污染相關(guān)健康效應(yīng)研究中應(yīng)用最廣泛的設(shè)計(jì)類型之一。病例交叉研究的數(shù)據(jù)中，通常只有病例沒有對(duì)照，為了研究暴露與研究對(duì)象疾病罹患的關(guān)系，通常以該患者健康效應(yīng)出現(xiàn)時(shí)間點(diǎn)前（和/或后）的某幾個(gè)時(shí)間點(diǎn)該患者的個(gè)體暴露狀態(tài)作為其自身對(duì)照，形成1 ∶ 1 或者1 ∶ M 的配比[1]。因此在統(tǒng)計(jì)方法上通常采用條件Logistic回歸對(duì)數(shù)據(jù)進(jìn)行分析。目前常規(guī)統(tǒng)計(jì)軟件如R、SAS 等均可完成條件Logistic 回歸分析[2]。相較于傳統(tǒng)數(shù)據(jù)統(tǒng)計(jì)工具，Python 作為一款流行的計(jì)算機(jī)語言，具有強(qiáng)大的通用性與可拓展性特點(diǎn)，特別是在控制其他軟件實(shí)現(xiàn)自動(dòng)化處理，智能化完成數(shù)據(jù)的采集、清洗、預(yù)處理以及數(shù)據(jù)挖掘等方面擁有明顯優(yōu)勢(shì)。但目前將Python 應(yīng)用于流行病學(xué)研究的案例較為少見，因此本文將應(yīng)用病例交叉研究的實(shí)例，探討Python 實(shí)現(xiàn)條件Logistic 回歸的過程，并比較其與R 和SAS 統(tǒng)計(jì)軟件在建模以及參數(shù)估計(jì)結(jié)果上的異同，以拓展Python 在流行病學(xué)領(lǐng)域中的應(yīng)用。

1 資料與方法

1.1 資料來源

案例資料來源于某地某年的住院首頁資料，根據(jù)國際疾病分類第10 版（ICD-10）對(duì)疾病進(jìn)行編碼，選擇肺部感染（ICD-10 代碼：J98.414）的患者作為研究對(duì)象。在這項(xiàng)研究中，共有3 216例肺部感染患者納入研究。

氣象數(shù)據(jù)來自中國氣象數(shù)據(jù)網(wǎng)，包含該地研究期間的每日平均溫度（℃）和日平均相對(duì)濕度（%），NO2濃度（μg/m3）資料來源于當(dāng)?shù)丨h(huán)境監(jiān)測(cè)中心。研究對(duì)象病例日（部分）溫度、濕度以及NO2濃度數(shù)據(jù)如表1所示。

表1 某地污染物信息Table 1.Pollutant information of a place

1.2 模型構(gòu)建

本研究Python 使用Cox 回歸對(duì)條件Logistic回歸進(jìn)行擬合，Cox 比例風(fēng)險(xiǎn)模型的基本形式為：

h(t,X)=h0(β'X)=h0(t)exp(β1X1+β2X2+…+βmXm)

h(t,X)是具有協(xié)變量X 的個(gè)體在時(shí)刻t 時(shí)的風(fēng)險(xiǎn)函數(shù)，t 為生存時(shí)間，X=(X1,X2,…,Xm)'是可能影響生存時(shí)間的有關(guān)因素。h0(t)是所有協(xié)變量取值為0 時(shí)的風(fēng)險(xiǎn)函數(shù)，稱為基線風(fēng)險(xiǎn)函數(shù)。β=(β1,β2,…,βm)為Cox 模型的回歸系數(shù)，是待估的回歸參數(shù)[3]。

2 結(jié)果

2.1 數(shù)據(jù)預(yù)處理

根據(jù)病例交叉研究的原理，在原始數(shù)據(jù)的基礎(chǔ)上，需要為每個(gè)病例日匹配3 至4 個(gè)對(duì)照日，使得對(duì)照組的特征與病例組的特征相似，以減少潛在混雜因素對(duì)研究結(jié)果的干擾。方法為選擇病例日和對(duì)照日為同年、同月的同一個(gè)星期幾，本研究對(duì)應(yīng)選擇了10 995 個(gè)對(duì)照日。SAS 擬合條件Logistic 回歸有兩種方法，分別為直接使用Logistic 回歸和借用Cox 回歸并定義分層變量后實(shí)現(xiàn)，兩者運(yùn)行結(jié)果相同。本文Python 和SAS 均采用分層Cox 風(fēng)險(xiǎn)比例模型進(jìn)行擬合，此法需新增一個(gè)時(shí)間變量（time），令time=1-case（病例日case編碼為1，對(duì)照日case編碼為0），設(shè)置原則為：病例日對(duì)應(yīng)的值小于對(duì)照日對(duì)應(yīng)的值即可。新設(shè)置的變量time 可作為Cox 風(fēng)險(xiǎn)比例模型中的生存時(shí)間變量，case 相當(dāng)于終檢變量[4]。匹配后數(shù)據(jù)信息（部分）見表2。

表2 預(yù)處理后某地污染物信息Table 2.Pollutant information of a place after pretreatment

2.2 模型擬合

條件Logistic 回歸在Python 中的實(shí)現(xiàn)首先需調(diào)用pandas 庫導(dǎo)入病例交叉數(shù)據(jù)并命名為“wb”，“columns.tolist()”為定義新列所用的函數(shù)，“col_name.insert(7,'time')”確定新列所在位置以及名稱，“df['time']=1-df['case']”是新列“time”生成的計(jì)算方式，本文中原則是time=1-case，最后生成新的數(shù)據(jù)集命名為“base_data”。具體命令如下：

import pandas as pd

wb=pd.read_excel(r"D: case-crossover.xlsx ",sheet_name="Sheet1")

col_name=wb.columns.tolist()

col_name.insert(7,'time')

df=wb.reindex(columns=col_name)

df['time']=1-df['case']

print(df)

base_data=df

print(base_data)

然后調(diào)用lifelines 庫中的CoxPHFitter 函數(shù)，“binglimerge.fit(base_data,'time',event_col='case',strata=['ID'])”，依次在括號(hào)中放入數(shù)據(jù)集、生存時(shí)間變量、終檢變量、分層變量。具體命令如下：

from lifelines import CoxPHFitter

binglimerge=CoxPHFitter()

binglimerge.fit(base_data,'time',event_col='case',strata=['ID'])

binglimerge.print_summary()

reults=binglimerge.summary

round(reults,7)

2.3 模型驗(yàn)證

R 4.2.1 軟件采用survival 包中的clogit 函數(shù)對(duì)條件Logistic 回歸模型進(jìn)行擬合，對(duì)實(shí)例資料的分析過程為：

library(readxl)

library(survival)

base_data<- read_xlsx("D:/ case-crossover.xlsx ")

options(digits = 8)

mod<-clogit(case ～ no2+temperature+humidity+strata(ID),base_data)

summary(mod)

AIC(mod)

SAS 9.1 版本采用phreg 過程步對(duì)條件Logistic回歸模型進(jìn)行擬合。其與Python 類似，在導(dǎo)入數(shù)據(jù)后，首先需對(duì)數(shù)據(jù)進(jìn)行預(yù)處理，增加新變量time(time=1-case)。對(duì)實(shí)例資料的分析過程為：

libname orange "D：sas data";

data base_data;set orange.binglijiaocha;time=1-case;run;

proc phreg data=base_data;

model time*case(0)= temperature humidity no2/ties=discrete risklimits;strata ID;

quit;

2.4 結(jié)果比較

Python 和SAS 在調(diào)用數(shù)據(jù)后，需要通過特定代碼運(yùn)行新增變量“time”，以方便采用分層Cox風(fēng)險(xiǎn)比例模型進(jìn)行擬合。R 則無需進(jìn)行上述操作，可直接通過clogit 函數(shù)實(shí)現(xiàn)模型擬合。與Python和SAS 相比，R 沒有默認(rèn)輸出AIC 值，需另運(yùn)行“AIC()”函數(shù)實(shí)現(xiàn)其結(jié)果輸出。

三款軟件輸出的主要結(jié)果基本相同（表3）。針對(duì)P 值的檢驗(yàn)方法上，R 與Python 輸出參數(shù)為z 值，SAS 輸出參數(shù)為χ2值，兩種檢驗(yàn)也是完全等價(jià)的（z 值的平方與χ2值相等）。

3 討論

在使用Python 和SAS 這兩款軟件時(shí)，本研究均用分層Cox 風(fēng)險(xiǎn)比例模型的運(yùn)行代碼來擬合條件Logistic 回歸模型，而R 語言則直接運(yùn)用survival 包中的clogit 函數(shù)進(jìn)行擬合，不用另對(duì)始變量進(jìn)行處理。其擬合原理是在分層Cox 模型中，各層的基線風(fēng)險(xiǎn)函數(shù)之間完全無關(guān)，而且Cox 風(fēng)險(xiǎn)比例模型在擬合時(shí)并沒有估計(jì)基線風(fēng)險(xiǎn)函數(shù)，只對(duì)各協(xié)變量的系數(shù)值β 進(jìn)行了估計(jì)，這和條件Logistic 回歸模型只求出系數(shù)值β 的思路一致[5]。有研究對(duì)Cox 比例風(fēng)險(xiǎn)模型總偏似然函數(shù)和條件Logistic 回歸分析的似然函數(shù)理論公式進(jìn)行推導(dǎo)后，發(fā)現(xiàn)它們完全等同[4]。本研究中三款軟件均采用極大似然估計(jì)法對(duì)參數(shù)進(jìn)行估計(jì)，其運(yùn)行結(jié)果完全相同，證實(shí)了擬合結(jié)果的可靠性。在衡量最優(yōu)模型的標(biāo)準(zhǔn)中，Python 以及SAS 軟件均自動(dòng)輸出AIC 值，但R 未自動(dòng)輸出該值，原因是R 調(diào)用的clogit 函數(shù)中不含衡量最優(yōu)模型標(biāo)準(zhǔn)的相關(guān)值的運(yùn)算代碼。另外，三款軟件輸出參數(shù)有z 值與χ2值的差異，其原因是不同軟件的開發(fā)人員在統(tǒng)計(jì)檢驗(yàn)傾向上不同，但z 值的平方等于χ2值，可以認(rèn)為Waldχ2檢驗(yàn)是等價(jià)于Z 檢驗(yàn)的[6]。

Python 作為一款面向?qū)ο蟮母呒?jí)編程語言，已經(jīng)成為最受歡迎的程序設(shè)計(jì)語言之一，在各行各業(yè)都發(fā)揮著重要的作用，常用于Web 應(yīng)用開發(fā)、人工智能、自動(dòng)化運(yùn)維、游戲開發(fā)等領(lǐng)域，其價(jià)值不可估量[7]。但在統(tǒng)計(jì)分析方面，Python的統(tǒng)計(jì)功能相對(duì)R 來說還比較薄弱，其自帶的處理功能和函數(shù)模型不及R 齊全，本研究中Python得采用Cox 風(fēng)險(xiǎn)比例模型去擬合條件Logistic 回歸，并且整個(gè)運(yùn)行過程相較另外兩款統(tǒng)計(jì)軟件都更復(fù)雜。在增加新變量方面，SAS 的操作步驟比Python 簡(jiǎn)潔很多。在可視化方面，Python 擁有Matplotlib 及Numpy 等繪圖庫[8]，可滿足可視化需求。R 作為一款為統(tǒng)計(jì)分析而設(shè)計(jì)的軟件，其可視化功能更為強(qiáng)大，它采用簡(jiǎn)潔的函數(shù)就能構(gòu)建各類圖形，并且在默認(rèn)條件下的繪圖品質(zhì)就能達(dá)到出版要求，但是R 在智能化方面以及非統(tǒng)計(jì)分析領(lǐng)域的應(yīng)用遠(yuǎn)不及Python。

綜上所述，將Python 應(yīng)用于統(tǒng)計(jì)分析領(lǐng)域，憑借其豐富的第三方庫以及快速運(yùn)算大數(shù)據(jù)的優(yōu)勢(shì)，能大大提高數(shù)據(jù)的智能化處理與分析效率。本研究使用Python 軟件實(shí)現(xiàn)了條件Logistic 回歸的統(tǒng)計(jì)建模，在實(shí)際研究中有一定的參考價(jià)值。