時(shí)間分層病例交叉研究的R軟件實(shí)現(xiàn)*

張彩霞 劉志東 張斐斐 丁國(guó)永 姜寶法Δ

·計(jì)算機(jī)應(yīng)用·

時(shí)間分層病例交叉研究的R軟件實(shí)現(xiàn)*

張彩霞1，2劉志東1，2張斐斐1，2丁國(guó)永3姜寶法1，2Δ

時(shí)間分層病例交叉研究是病例交叉研究的一種，適用于多變量的時(shí)間序列資料，現(xiàn)已廣泛應(yīng)用于環(huán)境污染、氣象因素、極端天氣事件對(duì)人群健康影響的研究。經(jīng)典的分析方法是條件 logistic回歸，常采用SPSS中的COX回歸模塊以及SAS的proc phreg、proc logistic過程實(shí)現(xiàn)，最近一些新的分析方法如Poisson回歸、條件Poisson回歸也被逐漸應(yīng)用到病例交叉研究中，但目前國(guó)內(nèi)尚無(wú)此類方法如何用軟件實(shí)現(xiàn)的相關(guān)文獻(xiàn)報(bào)道。此外，病例交叉研究的資料整理比較復(fù)雜，雖已有彭曉武［1］等有關(guān)應(yīng)用SAS整理雙向?qū)ΨQ病例交叉研究資料的方法介紹，但時(shí)間分層病例交叉資料的整理方法尚未見報(bào)道。因此，本文借助R軟件通過實(shí)例分析來(lái)實(shí)現(xiàn)基于時(shí)間序列的時(shí)間分層病例交叉資料的整理及統(tǒng)計(jì)分析，包括利用“season”包中的casecross函數(shù)實(shí)現(xiàn)時(shí)間分層病例交叉資料的條件logistic回歸分析，利用glm函數(shù)實(shí)現(xiàn)Poisson回歸分析，利用“gnm”包中的gnm函數(shù)實(shí)現(xiàn)條件Poisson回歸，最后利用山東省濟(jì)南市2005年8～9月期間洪水對(duì)細(xì)菌性痢疾每日發(fā)病例數(shù)影響的時(shí)間分層病例交叉研究數(shù)據(jù)集進(jìn)行驗(yàn)證。

時(shí)間分層病例交叉研究方法簡(jiǎn)介

病例交叉研究（case-crossover study）是由Maclure在1991年首次提出，它是一種用于研究短暫暴露對(duì)罕見急性病的瞬間影響的流行病學(xué)方法［2］。病例交叉研究可以看作病例對(duì)照研究的配對(duì)設(shè)計(jì)，它以病例自身作為對(duì)照不僅避免了選擇對(duì)照所引起的偏倚，而且可以避免各病例之間一些不可控制的因素（如年齡、智力、遺傳等）所引起的偏倚。當(dāng)對(duì)照數(shù)據(jù)是暴露人時(shí)時(shí)，病例交叉研究也可以看作回顧性的隊(duì)列研究，其分析可以按許多樣本只含有一個(gè)個(gè)體的隊(duì)列研究的薈萃分析來(lái)處理。病例交叉研究中交叉的思想來(lái)自于實(shí)驗(yàn)性交叉研究，后者實(shí)驗(yàn)對(duì)象經(jīng)歷兩個(gè)處理階段，而且各自作為自身對(duì)照以控制個(gè)體差異的影響；前者研究對(duì)象同時(shí)經(jīng)歷不同的時(shí)期即危險(xiǎn)期和對(duì)照期，通過比較個(gè)體在不同時(shí)期的暴露情況得出疾病與暴露的關(guān)聯(lián)［3］。時(shí)間分層病例交叉研究的基本原理是將時(shí)間進(jìn)行分層，病例期和對(duì)照期處于同一年、同一個(gè)月和同一個(gè)星期幾（the day of week），而且在同一時(shí)間層內(nèi)，幾個(gè)對(duì)照期是隨機(jī)分布的，病例期并非固定在某一位置。例如，假設(shè)病例期發(fā)生在2014年12月12日（星期五），則2014年12月其他的星期五均被選為對(duì)照期（見圖1）。

圖1 時(shí)間分層病例交叉設(shè)計(jì)對(duì)照期選擇示意圖

時(shí)間分層病例交叉設(shè)計(jì)相當(dāng)于按時(shí)間層匹配的病例對(duì)照研究，Janes［4］等通過統(tǒng)計(jì)模擬的方法發(fā)現(xiàn)，時(shí)間分層病例交叉研究可以同時(shí)控制季節(jié)性與星期幾效應(yīng)等混雜因素、消除時(shí)間趨勢(shì)偏倚，并能得到參數(shù)的無(wú)偏估計(jì)（條件logistic回歸）。而Poisson回歸和條件Poisson回歸通過設(shè)置啞變量（年、月、星期幾）的方法，同樣能達(dá)到按時(shí)間層匹配的目的。

實(shí)例與R軟件分析

1.資料概述

現(xiàn)以倫敦的一個(gè)空氣污染與人群死亡數(shù)據(jù)集為例演示時(shí)間分層病例交叉設(shè)計(jì)的軟件操作［5］。數(shù)據(jù)集包含從2002年1月1日到2006年12月31日期間，倫敦市每天的臭氧濃度（ozone）、溫度（temperature）、相對(duì)濕度（relative＿h(yuǎn)umidity）以及總死亡數(shù)（numdeaths）。所研究的暴露因素是臭氧濃度，健康結(jié)局是總死亡數(shù)，兩個(gè)潛在的混雜因素是溫度和相對(duì)濕度。數(shù)據(jù)集可以在 http://ije.oxfordjournals.org/content/suppl/2013/05/30/dyt092.DC1下載。

2.數(shù)據(jù)導(dǎo)入及變量預(yù)處理

library（foreign）＃加載 foreign包

data＜-read.dta（“l(fā)ondondataset.dta”）＃導(dǎo)入數(shù)據(jù)集

表1 倫敦空氣污染與人群死亡數(shù)據(jù)集（部分）

data＄temp＜-data＄temperature＃重命名溫度

data＄rh＜-data＄relative＿h(yuǎn)umidity＃重命名相對(duì)濕度

data＄dow＜-as.factor（weekdays（data＄date））＃添加星期元啞變量

data＄month＜-as.factor（months（data＄date））＃添加月份變量

data＄year＜-as.factor（format（data＄date，format＝“%Y”））＃添加年份變量

data＄stratum ＜-as.factor（data ＄year:data＄month:data＄dow）＃添加時(shí)間層變量

3.利用“season”包中的casecross函數(shù)實(shí)現(xiàn)條件logistic回歸

library（season）＃加載 season包

model1＝casecross（numdeaths～ozone＋temp＋rh，matchdow＝TRUE，stratamonth＝TRUE，data＝data）

＃用 matchdow＝TRUE，stratamonth＝TRUE匹配同一年同一月同一個(gè)星期幾

summary（model1）＃匯總結(jié)果

4.利用glm函數(shù)實(shí)現(xiàn)泊松回歸

model2＜-glm（numdeaths～ozone＋temp＋rh＋factor（stratum），data＝data，family＝poisson）

summary（model2）

5.利用“gnm”包中的gnm函數(shù)實(shí)現(xiàn)條件Poisson回歸

library（gnm）＃加載 gnm包

model3＜-gnm（numdeaths～ozone＋temp＋rh，data＝data，family＝poisson，eliminate＝factor（stratum））

summary（model3）

表2 R軟件casecross、glm、gnm函數(shù)過程參數(shù)估計(jì)結(jié)果

上述三種方法的結(jié)果列于表2，可見條件logistic回歸、Poisson回歸和條件Poisson回歸3種方法的結(jié)果是完全一致的，其中條件logistic回歸得到的結(jié)果為比值比（odds ratio，OR），Poisson回歸和條件 Poisson回歸得到的是危險(xiǎn)度比（risk ratio，RR）。

6.現(xiàn)應(yīng)用以上R程序研究山東省濟(jì)南市2005年8～9月期間洪水對(duì)細(xì)菌性痢疾每日發(fā)病例數(shù)的影響，暴露因素為洪水（flood，用二分類變量0、1表示，研究期間發(fā)生了兩次洪水分別為8月1～2號(hào)和9月19～21號(hào)），健康結(jié)局是濟(jì)南市細(xì)菌性痢疾每日發(fā)病例數(shù)（N），欲控制的兩個(gè)潛在混雜因素為平均溫度（AT）和平均相對(duì)濕度（ARH）。具體數(shù)據(jù)見表3。

表3 濟(jì)南市2005年8－9月洪水狀況、氣象因素與菌痢發(fā)病數(shù)數(shù)據(jù)集（部分）

將上文中介紹的條件logistic回歸、Poisson回歸和條件Poisson回歸3種方法應(yīng)用于研究山東省濟(jì)南市2005年8～9月期間洪水對(duì)細(xì)菌性痢疾每日發(fā)病例數(shù)的影響，得到的結(jié)果是一致的（見表4）。

表4 洪水對(duì)細(xì)菌性痢疾發(fā)病數(shù)影響的參數(shù)估計(jì)結(jié)果

7.過離散、自相關(guān)、滯后等相關(guān)問題，以Poisson回歸為例

（1）控制過離散:

model4＜-glm（numdeaths～ozone＋temp＋rh＋factor（stratum），data＝data，family＝quasipoisson）

summary（model4）

用quasipoisson代替poisson分布族控制過離散現(xiàn)象

（2）控制自相關(guān):

library（tsModel）＃加載 tsModel包

reslag1 ＜-Lag（residuals（model4，type＝“deviance”），1）＃提取模型殘差項(xiàng)

model5＜-glm（numdeaths～ozone＋temp＋rh＋reslag1，data＝data，family＝quasipoisson）＃通過加入殘差項(xiàng)控制1階自相關(guān)

summary（model5）

（3）單滯后模型:

library（Epi）＃加載 Epi包

library（tsModel）＃加載 tsModel包

tablag1＜-matrix（NA，7＋1，3，dimnames＝list（paste（“Lag”，0:7），c（“RR”，“ci.low”，“ci.hi”）））＃建立7行3列的表格

for（i in 0:7）｛numdeathslag＜-Lag（data＄numdeaths，-i）

model6＜-glm（numdeathslag～ozone＋temp＋rh＋factor（stratum），data＝data，family＝quasipoisson）

tablag1［i＋1，］＜-ci.lin（model6，subset＝“ozone”，Exp＝T）［5:7］｝

tablag1＃建立滯后變量，并依次估計(jì)不同滯后的參數(shù)，結(jié)果放在tablag1表格

（4）分布滯后模型:

library（dlnm）＃加載 dlnm包

cbozone＜-crossbasis（data＄ozone，lag＝c（0，7），argvar＝list（type＝“l(fā)in”，cen＝0），arglag＝list（type＝“integer”））＃建立臭氧與滯后的交叉基

model7＜-glm（numdeaths～cbozone＋temp＋rh＋factor（stratum），data＝data，family＝quasipoisson）

pred＜-crosspred（cbozone，model7，at＝10）

＃cen＝0指臭氧濃度參照水平為0 ug/m3，at＝10指相對(duì)于0μg/m3，臭氧濃度為10ug/m3的效應(yīng)

tablag2 ＜-with（pred，t（rbind（matRRfit，matRRlow，matRRhigh）））＃提取 RR值及可信區(qū)間

colnames（tablag2） ＜-c（“RR”，“ci.low”，“ci.hi”）

tablag2

pred＄allRRfit；pred＄allRRlow；pred＄allRRhigh＃不同滯后的累積效應(yīng)

討論與小結(jié)

本文應(yīng)用R軟件進(jìn)行時(shí)間分層病例交叉資料的統(tǒng)計(jì)分析，簡(jiǎn)化了時(shí)間序列資料繁瑣的整理，增加了統(tǒng)計(jì)方法選擇的靈活性。R軟件開源、免費(fèi)，近年來(lái)在數(shù)據(jù)探索、統(tǒng)計(jì)分析、作圖等領(lǐng)域應(yīng)用廣泛，其幫助功能十分強(qiáng)大，本文中出現(xiàn)的所有命令、參數(shù)都可以通過幫助系統(tǒng)找到其詳細(xì)解釋。針對(duì)3種不同的統(tǒng)計(jì)方法進(jìn)行R軟件實(shí)例分析后，結(jié)果顯示三種方法參數(shù)估計(jì)的結(jié)果是一致的。利用我們的洪水發(fā)生情況與細(xì)菌性痢疾發(fā)病數(shù)數(shù)據(jù)集進(jìn)行驗(yàn)證也得到了相同的結(jié)果。國(guó)外多項(xiàng)研究亦表明，時(shí)間分層病例交叉設(shè)計(jì)的條件logistic回歸與泊松回歸分析是等價(jià)的［6－8］。條件logistic回歸是病例交叉設(shè)計(jì)最常用的統(tǒng)計(jì)分析方法，但其并不能控制時(shí)間序列數(shù)據(jù)的自相關(guān)、過離散等問題，對(duì)此我們可以利用Poisson回歸來(lái)控制這些問題，但Poisson回歸也存在估計(jì)參數(shù)較多的缺點(diǎn)。條件Poisson回歸是Poisson回歸的進(jìn)一步發(fā)展，其顯著優(yōu)點(diǎn)是減少了需要估計(jì)的參數(shù)，縮短了軟件運(yùn)算的時(shí)間，同時(shí)又不影響參數(shù)估計(jì)的準(zhǔn)確性［9］。本文同時(shí)針對(duì)相關(guān)研究中的自相關(guān)、過離散、滯后等問題進(jìn)行了R軟件操作，提供了一個(gè)簡(jiǎn)便靈活的R程序。讀者進(jìn)行時(shí)間分層病例交叉設(shè)計(jì)的統(tǒng)計(jì)分析時(shí)，建議結(jié)合數(shù)據(jù)特點(diǎn)選擇合適的方法靈活分析。

［1］彭曉武，余松林，相紅，等.用SAS程序整理病例交叉設(shè)計(jì)資料.中國(guó)衛(wèi)生統(tǒng)計(jì)，2012，29（1）:135-138.

［2］Maclure M.The case-crossover design:amethod for studying transient effects on the risk of acute events.Am J Epidemiol，1991，133（2）:144-53.

［3］胡以松.病例交叉研究.疾病控制雜志，2001，5（4）:341-343.

［4］Janes H，Sheppard L，Lum ley T.Case-crossover analyses of air pollution exposure data:referent selection strategies and their implications for bias.Epidemiology，2005，16（6）:717-726.

［5］Bhaskaran K，Gasparrini A，Hajat S，et al.Time series regression studies in environmental epidemiology.International journal of epidemiology，2013，42（4）:1187-1195.

［6］Levy D，Lum ley T，Sheppard L，etal.Referent selection in case-crossover analyses of acute health effects of air pollution.Epidemiology，2001，12（2）:186-192.

［7］Lu Y，Zeger SL.On the equivalence of case-crossover and time series methods in environmental epidemiology.Biostatistics，2007，8（2）:337-344.

［8］NavidiW.Poisson Regression and the Case-Crossover Design:Similarities and Differences.Communications in Statistics（Theory and Methods），2008，37（2）:213-220.

［9］Armstrong BG，Gasparrini A，Tobias A.Conditional Poisson models:a flexible alternative to conditional logistic case cross-over analysis.BMC Medical Research Methodology，2014，14（1）:122.

國(guó)家重大科學(xué)研究計(jì)劃（973計(jì)劃）項(xiàng)目（2012CB955502）

1.山東大學(xué)公共衛(wèi)生學(xué)院流行病學(xué)系（250012）

2.山東大學(xué)氣候變化與健康研究中心

3.泰山醫(yī)學(xué)院公共衛(wèi)生學(xué)院

△通信作者:姜寶法，E-mail:bjiang@sdu.edu.cn

（責(zé)任編輯:鄧 妍）