許敏銳 強(qiáng)德仁 周義紅 石素逸 秦 晶 陶 源

應(yīng)用R軟件進(jìn)行l(wèi)ogistic回歸模型的交互作用分析*

許敏銳1強(qiáng)德仁1周義紅1石素逸1秦 晶1陶 源2△

目的 應(yīng)用R軟件進(jìn)行l(wèi)ogistic回歸模型的交互作用分析，為探討交互作用提供依據(jù)。方法 使用R軟件，編寫程序?qū)崿F(xiàn)logistic或Cox回歸模型三個(gè)評(píng)價(jià)相加交互作用的指標(biāo)及其可信區(qū)間的計(jì)算。結(jié)果 生物學(xué)交互作用的評(píng)價(jià)應(yīng)該基于是否有相加交互作用，而流行病學(xué)研究中常運(yùn)用logistic回歸模型，并納入乘積項(xiàng)分析因素間交互作用，其是否有意義僅反映相乘交互作用，并不能反映兩因素間相加或生物學(xué)交互作用的有無(wú)。本文通過(guò)實(shí)例分析，調(diào)用基于R軟件編寫的interact程序，可以直接計(jì)算出logistic或Cox回歸模型的三個(gè)交互作用評(píng)價(jià)指標(biāo)（RERI、AP、SI）及其可信區(qū)間；并將結(jié)果與運(yùn)用Andersson編制的Excel計(jì)算結(jié)果相比較，驗(yàn)證了本程序的科學(xué)性和準(zhǔn)確性。結(jié)論 應(yīng)用R軟件編制程序，可實(shí)現(xiàn)logistic回歸模型因素間交互作用和可信區(qū)間的計(jì)算，為流行病學(xué)研究人員分析生物學(xué)交互作用提供依據(jù)。

logistic回歸 交互作用 R軟件

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

基本原理和方法

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

Rothman用于評(píng)價(jià)相加交互作用的三個(gè)指標(biāo)，①相對(duì)超危險(xiǎn)度比：RERI＝OR_A1B1-OR_A0B1-OR_A1B0+ 1；②歸因比：AP＝RERI／OR_A1B1；③交互作用指數(shù)SI＝（OR_A1B1-1）／［（OR_A0B1-1）+（OR_A1B0-1）］。如果兩因素?zé)o相加交互作用，則RERI和AP的可信區(qū)間應(yīng)包含0，SI的可信區(qū)間應(yīng)包含1。

1．交互作用指標(biāo)的點(diǎn)估計(jì)：可通過(guò)以下兩種方法，建立logistic回歸模型計(jì)算OR_A1B1、OR_A0B1和OR_A1B0，代入交互作用指標(biāo)的計(jì)算公式。

（1）用兩因素A、B及乘積項(xiàng)A×B構(gòu)建logistic回歸模型1。則有

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

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

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

程序簡(jiǎn)介

R軟件作為免費(fèi)開(kāi)源的軟件，能夠通過(guò)編寫自定義程序?qū)崿F(xiàn)一些功能，得到越來(lái)越多人的認(rèn)可和使用。通過(guò)R軟件的自定義function，根據(jù)交互作用的基本原理和方法，編寫程序語(yǔ)句，調(diào)用程序可以快捷方便地得出交互作用的結(jié)果。我們將交互作用計(jì)算的程序（interact）已經(jīng)在文后的附錄中給出，供有興趣的研究者參考，對(duì)于編程基礎(chǔ)較為薄弱者，可以直接按照附錄的編程語(yǔ)句，在R軟件中編寫并運(yùn)行。根據(jù)需要探討交互作用變量建立logistic回歸模型，并調(diào)用編寫的交互作用程序interact，即可得出OR值、RERI、AP和SI值，同時(shí)給出交互作用示意圖。

實(shí)例分析

1．模擬數(shù)據(jù)庫(kù) 以模擬的inter數(shù)據(jù)庫(kù)為例，設(shè)置兩個(gè)分類變量A（0為無(wú)暴露，1為暴露）、B（0為無(wú)暴露，1為暴露），一個(gè)結(jié)局變量case（0為無(wú)結(jié)局，1為結(jié)局），兩個(gè)混雜調(diào)整變量（x1，x2），具體見(jiàn)表1。

表1 模擬數(shù)據(jù)庫(kù)inter基本情況

2．設(shè)立新的啞變量，建立logistic回歸模型 根據(jù)兩個(gè)變量A、B設(shè)置一個(gè)新的啞變量A_B（A0B0表示A＝0且B＝0，分析時(shí)作為參照組，A1B0表示A＝1且B＝0，A0B1表示A＝0且B＝1，A1B1表示A＝1且B＝1），以新設(shè)置的啞變量A_B建立logistic回歸模型Iglm，可在模型中放入需要調(diào)整的變量（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構(gòu)建logistic回歸模型（注意：模型中A_B變量在前，調(diào)整變量在后）＃

3．計(jì)算交互作用和可信區(qū)間 在R軟件中運(yùn)行編寫的interact程序并調(diào)用，即可得出OR值、RERI、AP和SI值，同時(shí)給出交互作用示意圖（圖1）。運(yùn)行結(jié)果顯示，調(diào)整了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的可信區(qū)間應(yīng)大于0，SI的可信區(qū)間應(yīng)大于1，則說(shuō)明A、B之間存在交互作用。RERI和SI意義相同，AP表示全部病例中可歸因于兩因素交互作用的病例所占的比例，本例AP（95%CI）＝0．549（0．303，0．794），說(shuō)明全部病例中歸因于A和B的交互作用所引起的病例占54．9%。

圖1 R軟件調(diào)用interact程序計(jì)算交互作用結(jié)果和交互作用示意圖

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

討 論

統(tǒng)計(jì)學(xué)交互作用和生物學(xué)交互作用在病因?qū)W研究中有一定的區(qū)別，不能等同于統(tǒng)計(jì)模型中乘積項(xiàng)的分析。統(tǒng)計(jì)學(xué)交互作用是指在統(tǒng)計(jì)模型中納入乘積項(xiàng)的意義，在線性模型中是加法模型，乘積項(xiàng)表示有無(wú)相加交互作用，而對(duì)于logistic或Cox等乘法模型，乘積項(xiàng)表示有無(wú)相乘交互作用。生物學(xué)交互作用是指兩因素且同時(shí)存在時(shí)，是否具有在生物機(jī)制上聯(lián)合作用，包括協(xié)同作用和拮抗作用［1］。

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

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

表2 交互作用指標(biāo)和可信區(qū)間Excel計(jì)算結(jié)果

［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回歸模型中交互作用的分析及評(píng)價(jià)．中華流行病學(xué)雜志，2008，29（9）：934-937．

［7］張俊國(guó)，劉麗，李麗霞，等．懲罰廣義線性模型在遺傳關(guān)聯(lián)研究中的應(yīng)用及R軟件實(shí)現(xiàn)．中國(guó)衛(wèi)生統(tǒng)計(jì)，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回歸模型中連續(xù)變量交互作用的分析．中華流行病學(xué)雜志，2010，31（7）：812-814．

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

常州市武進(jìn)區(qū)科技發(fā)展計(jì)劃項(xiàng)目（WS201432）

1．常州市武進(jìn)區(qū)疾病預(yù)防控制中心（213164）

2．常州市第一人民醫(yī)院（213000）

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