应用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