時間分層病例交叉研究的R軟件實現*

張彩霞 劉志東 張斐斐 丁國永 姜寶法Δ

·計算機應用·

時間分層病例交叉研究的R軟件實現*

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

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

時間分層病例交叉研究方法簡介

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

圖1 時間分層病例交叉設計對照期選擇示意圖

時間分層病例交叉設計相當于按時間層匹配的病例對照研究，Janes［4］等通過統計模擬的方法發現，時間分層病例交叉研究可以同時控制季節性與星期幾效應等混雜因素、消除時間趨勢偏倚，并能得到參數的無偏估計（條件logistic回歸）。而Poisson回歸和條件Poisson回歸通過設置啞變量（年、月、星期幾）的方法，同樣能達到按時間層匹配的目的。

實例與R軟件分析

1.資料概述

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

2.數據導入及變量預處理

library（foreign）＃加載 foreign包

data＜-read.dta（“londondataset.dta”）＃導入數據集

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

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

data＄rh＜-data＄relative＿humidity＃重命名相對濕度

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）＃添加時間層變量

3.利用“season”包中的casecross函數實現條件logistic回歸

library（season）＃加載 season包

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

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

summary（model1）＃匯總結果

4.利用glm函數實現泊松回歸

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

summary（model2）

5.利用“gnm”包中的gnm函數實現條件Poisson回歸

library（gnm）＃加載 gnm包

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

summary（model3）

表2 R軟件casecross、glm、gnm函數過程參數估計結果

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

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

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

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

表4 洪水對細菌性痢疾發病數影響的參數估計結果

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

（1）控制過離散:

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

summary（model4）

用quasipoisson代替poisson分布族控制過離散現象

（2）控制自相關:

library（tsModel）＃加載 tsModel包

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

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

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＃建立滯后變量，并依次估計不同滯后的參數，結果放在tablag1表格

（4）分布滯后模型:

library（dlnm）＃加載 dlnm包

cbozone＜-crossbasis（data＄ozone，lag＝c（0，7），argvar＝list（type＝“lin”，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指相對于0μg/m3，臭氧濃度為10ug/m3的效應

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

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

tablag2

pred＄allRRfit；pred＄allRRlow；pred＄allRRhigh＃不同滯后的累積效應

討論與小結

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

［1］彭曉武，余松林，相紅，等.用SAS程序整理病例交叉設計資料.中國衛生統計，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.

國家重大科學研究計劃（973計劃）項目（2012CB955502）

1.山東大學公共衛生學院流行病學系（250012）

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

3.泰山醫學院公共衛生學院

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

（責任編輯:鄧 妍）