應用R軟件進行logistic回歸模型的交互作用分析*

許敏銳 強德仁 周義紅 石素逸 秦 晶 陶 源

應用R軟件進行logistic回歸模型的交互作用分析*

許敏銳1強德仁1周義紅1石素逸1秦 晶1陶 源2△

目的 應用R軟件進行logistic回歸模型的交互作用分析，為探討交互作用提供依據。方法 使用R軟件，編寫程序實現logistic或Cox回歸模型三個評價相加交互作用的指標及其可信區間的計算。結果 生物學交互作用的評價應該基于是否有相加交互作用，而流行病學研究中常運用logistic回歸模型，并納入乘積項分析因素間交互作用，其是否有意義僅反映相乘交互作用，并不能反映兩因素間相加或生物學交互作用的有無。本文通過實例分析，調用基于R軟件編寫的interact程序，可以直接計算出logistic或Cox回歸模型的三個交互作用評價指標（RERI、AP、SI）及其可信區間；并將結果與運用Andersson編制的Excel計算結果相比較，驗證了本程序的科學性和準確性。結論 應用R軟件編制程序，可實現logistic回歸模型因素間交互作用和可信區間的計算，為流行病學研究人員分析生物學交互作用提供依據。

logistic回歸 交互作用 R軟件

在統計分析中交互作用是指某因素的作用隨其他因素水平變化而變化，兩因素共同作用不等于兩因素單獨作用之和（相加交互作用）或之積（相乘交互作用）［1］。目前流行病學研究在分析因素間交互作用時，常采用納入因素乘積項到回歸方程中的方法。一般認為，線性回歸模型為相加模型，反映因素間是否有相加交互作用，而logistic回歸或Cox回歸模型為相乘模型，反映因素間是否有相乘交互作用［2］。Rothman指出logistic或Cox回歸模型中乘積項無統計學意義，并不表示兩因素無相加交互作用，也不表示無生物學交互作用，并從理論上探討了用于評價因素間是否有區別于相乘交互作用的相加交互作用，以及三個評價指標：相對超危險度比（the relative excess risk due to interaction，RERI）、歸因比（the attributable proportion due to interaction，AP）和交互作用指數（the synergy index，SI）的構造和計算方法［3］。本研究以logistic回歸分析為例，利用R軟件（http：／／www．r-project．org／）編寫計算程序，可無需計算變量間的方差和協方差矩陣，直接給出交互作用和可信區間，并將結果同Andersson等［4］編制的Excel計算結果進行比較，以期為流行病學研究中評價相加交互作用提供便捷的方法。

基本原理和方法

以最簡單的兩因素兩水平為例。假設兩暴露因子分別為A、B，1表示因素存在，0表示因素不存在，因變量為疾病的發生與否。logistic回歸模型得到的OR值作為相對危險度（RR）的估計值，OR_A0B0表示A、B都不存在時發病的OR值，分析時作為參照組；OR_A1B0表示僅A存在、B不存在時發病的OR值；OR_A0B1表示A不存在、僅B存在時發病的OR值；OR_A1B1表示A、B共同存在時發病的OR值。

Rothman用于評價相加交互作用的三個指標，①相對超危險度比：RERI＝OR_A1B1-OR_A0B1-OR_A1B0+ 1；②歸因比：AP＝RERI／OR_A1B1；③交互作用指數SI＝（OR_A1B1-1）／［（OR_A0B1-1）+（OR_A1B0-1）］。如果兩因素無相加交互作用，則RERI和AP的可信區間應包含0，SI的可信區間應包含1。

1．交互作用指標的點估計：可通過以下兩種方法，建立logistic回歸模型計算OR_A1B1、OR_A0B1和OR_A1B0，代入交互作用指標的計算公式。

（1）用兩因素A、B及乘積項A×B構建logistic回歸模型1。則有

（2）根據兩因素A、B，建立新的交互作用啞變量A _B（A0B0表示A＝0且B＝0，分析時作為參照組，A0B1表示A＝0且B＝1，A0B1表示A＝0且B＝1，A1B1表示A＝1且B＝1），構建logistic回歸模型2。

模型1和2中的β1、β2相同，而模型2中的β3等于模型1中的β1+β2+β3。

2．交互作用指標的區間估計：運用Hosmer［5］介紹的delta方法估計可信區間，利用R軟件編寫計算程序，計算交互作用和可信區間。同時介紹使用Andersson等［4］編制的Excel計算交互作用和可信區間的方法，并將兩者計算的結果進行比較，判定計算方法的科學性。計算所需的因素的方差和協方差可由R軟件建立logistic回歸模型得到。

程序簡介

R軟件作為免費開源的軟件，能夠通過編寫自定義程序實現一些功能，得到越來越多人的認可和使用。通過R軟件的自定義function，根據交互作用的基本原理和方法，編寫程序語句，調用程序可以快捷方便地得出交互作用的結果。我們將交互作用計算的程序（interact）已經在文后的附錄中給出，供有興趣的研究者參考，對于編程基礎較為薄弱者，可以直接按照附錄的編程語句，在R軟件中編寫并運行。根據需要探討交互作用變量建立logistic回歸模型，并調用編寫的交互作用程序interact，即可得出OR值、RERI、AP和SI值，同時給出交互作用示意圖。

實例分析

1．模擬數據庫 以模擬的inter數據庫為例，設置兩個分類變量A（0為無暴露，1為暴露）、B（0為無暴露，1為暴露），一個結局變量case（0為無結局，1為結局），兩個混雜調整變量（x1，x2），具體見表1。

表1 模擬數據庫inter基本情況

2．設立新的啞變量，建立logistic回歸模型 根據兩個變量A、B設置一個新的啞變量A_B（A0B0表示A＝0且B＝0，分析時作為參照組，A1B0表示A＝1且B＝0，A0B1表示A＝0且B＝1，A1B1表示A＝1且B＝1），以新設置的啞變量A_B建立logistic回歸模型Iglm，可在模型中放入需要調整的變量（x1，x2）。

inter＄A_B＜-ifelse（inter＄A＝＝0＆inter＄B＝＝0，“A1B0”，ifelse（inter＄A＝＝0＆inter＄B＝＝1，“A0B1”，ifelse（inter＄A＝＝1＆inter＄B＝＝0，“A1B0”，“A1B1”）））＃建立新的啞變量A_B＃

Iglm＜-glm（case～as．factor（A_B）+x1+x2，family＝binomial，data＝inter）＃以A_B構建logistic回歸模型（注意：模型中A_B變量在前，調整變量在后）＃

3．計算交互作用和可信區間 在R軟件中運行編寫的interact程序并調用，即可得出OR值、RERI、AP和SI值，同時給出交互作用示意圖（圖1）。運行結果顯示，調整了x1和x2因素后，以A0B0組為參照組，OR_A0B1＝1．828，OR_A1B0＝2．912，OR_A1B1＝8．290；RERI（95%CI）為4．550（0．361，8．739），AP（95%CI）為0．549（0．303，0．794），SI（95%CI）為2．660（1．382，5．121）；RERI和AP的可信區間應大于0，SI的可信區間應大于1，則說明A、B之間存在交互作用。RERI和SI意義相同，AP表示全部病例中可歸因于兩因素交互作用的病例所占的比例，本例AP（95%CI）＝0．549（0．303，0．794），說明全部病例中歸因于A和B的交互作用所引起的病例占54．9%。

圖1 R軟件調用interact程序計算交互作用結果和交互作用示意圖

4．用Andersson編制的Excel計算交互作用和可信區間 將新設置的啞變量A_B建立logistic回歸模型Iglm，可在模型中放入需要調整的變量（x1，x2），并將啞變量A_B的回歸系數β1、β2、β3，以及方差和協方差矩陣，輸入Andersson編制的Excel中，可得到RERI、AP和SI的點估計、95%CI（表2）。對比可見，兩種計算的結果完全一致。

討 論

統計學交互作用和生物學交互作用在病因學研究中有一定的區別，不能等同于統計模型中乘積項的分析。統計學交互作用是指在統計模型中納入乘積項的意義，在線性模型中是加法模型，乘積項表示有無相加交互作用，而對于logistic或Cox等乘法模型，乘積項表示有無相乘交互作用。生物學交互作用是指兩因素且同時存在時，是否具有在生物機制上聯合作用，包括協同作用和拮抗作用［1］。

Rothman［3］提出對于生物學交互作用的評價應基于相加尺度，對logistic、Cox回歸等相乘模型構建了本文介紹的三項指標，用于評價因素間是否有相加交互作用。邱宏等［6］介紹了在SPSS中運用Multinomial logistic過程構建回歸建模，將模型參數估計值和因素間的協方差矩陣帶入Andersson等編制的Excel計算表計算交互作用和可信區間，其操作過程較為復雜，在填寫協方差矩陣的時候易出錯，尤其在探討多個因素之間兩兩交互作用時，可以節省大量的時間，避免出錯。R軟件作為一種免費的軟件，應用越來越廣泛［7］，目前尚無運用R軟件進行二分類變量logistic回歸模型交互作用分析的使用介紹。本研究應用R軟件，編寫計算相加交互作用和可信區間的程序，通過調用程序即可得出三個相加交互作用指標的點估計和可信區間，為研究人員分析交互作用提供參考依據。

本方法僅適用于兩因素二分類的相加交互作用評價，在因素變量設置時，一般以風險的一類作為暴露組，尤其是在保護因素時，應當將無暴露設置為1，有暴露設置為0，以避免解釋上混亂。當因素變量為多分類或連續變量時，該計算方法以及Andersson編制的Excel法均不適用。對此Assmann等［8］提出Bootstrap法，在原始數據中做重復千次、萬次的模擬隨機抽樣，估計RERI，AP和SI及其可信區間。使用Bootstrap法在R軟件中分析兩個連續自變量或連續變量與分類變量間的交互作用的方法，邱宏等［9］已經做了詳細介紹，可供流行病學交互作用分析提供參考和借鑒。

表2 交互作用指標和可信區間Excel計算結果

［1］Ahlbom A，Alfredsson L．Interaction：A word w ith two meanings creates confusion．European Journal of Epidemiology，2005，20（7）：563-564．

［2］Knol MJ，Vand TI，Grobbee DE，et al．Estimating interaction on an additive scale between continuous determinants in a logistic regression model．International Journal of Epidemiology，2007，36（5）：1111-1118．

［3］Rothman KJ．Epidem iology：An introduction．New York：Oxford University Press，2002：168-180．

［4］Andersson T，Alfredsson L，K?llberg H，et al．Calculating measures of biological interaction．European Journal of Epidemiology，2005，20（7）：575-579．

［5］Hosmer DW，Lemeshow S．Confidence interval estimation of interaction．Epidemiology，1992，3（5）：333-338．

［6］邱宏，余德新，王曉蓉，等．logistic回歸模型中交互作用的分析及評價．中華流行病學雜志，2008，29（9）：934-937．

［7］張俊國，劉麗，李麗霞，等．懲罰廣義線性模型在遺傳關聯研究中的應用及R軟件實現．中國衛生統計，2016，33（4）：582-586．

［8］Assmann SF，Hosmer DW，Lemeshow S，et al．Confidence intervals formeasures of interaction．Epidem iology，1996，7（3）：286-90．

［9］邱宏，余德新，謝立亞，等．logistic回歸模型中連續變量交互作用的分析．中華流行病學雜志，2010，31（7）：812-814．

（責任編輯：鄧 妍）

常州市武進區科技發展計劃項目（WS201432）

1．常州市武進區疾病預防控制中心（213164）

2．常州市第一人民醫院（213000）

△通信作者：陶源，E-mail：taodazanze＠163．com