国外流行病学或临床研究越来越多将交互作用纳入标准化分析程序。国内常用logistic回归主效应模型中纳入因素的乘积交互项，并根据乘积项的意义判断有无交互作用。此类分析误区是无法准确理解统计建模交互项，错误将相乘交互等于相加交互。Rothman等^[1]详细评价了相加交互与相乘交互的差异，认为统计建模中一般线性模型交互项反映的是因素间相加交互作用（additive interactions，INTA），而logistic等广义线性模型则反映因素间统计学上相乘交互作用（multiplicative interaction，INTM），即logistic乘积项仅反映统计学上的交互，而只有在生物学机制上病因因素间存在相加交互作用才可解释为相加交互作用。定量评价流行病学研究中暴露因素间及暴露因素与基因间相加交互作用需要根据3项指标，即交互对比度（interaction contrast ratio，ICR）又称相对超危险度比（relative excess risk due to interaction，RERI），归因比（attributable proportion due to interaction，AP）和交互作用指数（the synergy index，S）来计算。其中ICR的点估计较为简单，而95%CI较为复杂。当前已有文献提供计算ICR程序，但其算法冗繁且实际分析需要大幅修改原程序使推广受限^{[2, 3, 4, 5]}。Andersson等^[6]制作Excel版本工具极大推广了ICR 95%CI运用，邱宏等^[7]亦分别对分类交互Delta法进行了阐述。但Andersson等的方法不适合大数据多变量程序化分析，且SPSS进行条件logistic运用Cox分析模块无法输出协方差矩阵，Delta法估计置信区间易造成假阴性结果。Zou^[8]亦阐述了Wald置信区间法需基于风险比对称的前提，故Richardson推荐采用轮廓似然置信区间（profile likelihood confidence intervals，PL），利用Wald的ICR点估计重构ICR 95%CI。logistic回归是最常见的报道交互作用模型^[9]，我们认为同时利用Delta法、Wald法、PL法可更全面评价相互作用的效果，为此本文提供SAS 9.4宏程序，实现非条件/条件logistic和Cox比例风险模型ICR 95%CI的估计，旨在为临床与遗传流行病大数据分析因素间交互作用提供简便的方法。

Rothman等多数流行病学者认为，生物学交互作用的评价应该基于相加尺度的INTA而非相乘尺度的INTM，并以此提出构建了ICR、AP、S 3个指标，用于评价因素间是否有区别于相乘交互作用的相加交互作用，进而为生物学交互作用的评价提供依据。如果两因素有相加交互作用，则ICR 95%CI、AP 95%CI应不包含0，S 95%CI应不包含l。假设研究多风险因素中交互作用的两暴露因素为AB，则OR00表示AB均无暴露，即OR₀₀＝1；OR10表示A暴露、B无暴露，OR₀₁表示A无暴露、B暴露，OR₁₁表示A、B均暴露。Rothman评价相乘交互作用INT_M＝OR_A×B＝OR₁₁/（OR₁₀×OR₀₁），即logistic回归乘积项95%CI不包含1，表明有相乘交互作用。

相加交互作用3项指标定义为

如果两因素无相加交互作用，则ICR和AP的可信区间应包含0，S的可信区间应包含1。则二分类logistic回归模型中

以2004-2014年广东省先天性心脏病（先心）监测网数据库病例对照资料为例^{[10, 11]}，分析孕母是否被动吸烟（自报围孕期家庭、工作任一环境中接触吸烟平均时长＞15 min/d且持续1周以上）与家庭人均月收入有无相乘和相加交互作用（表 1），并比对采用Andersson分析方法的结果。

表 1 广东省先心监测网病例对照研究中家庭人均月收入

表选项

如表 2，计算得到OR₁₀＝1.383（95%CI：1.211～1.579），OR₀₁＝1.268（95%CI：1.160～1.387），OR₁₁＝2.443（95%CI：2.048～2.914），INT_M＝1.393（95%CI：1.116～1.739）。采用Delta法、Wald法和PL法分别估计ICR、AP、S及其95%CI，其交互作用见图 1和表 3。本例INT_M 95%CI＞1有相乘交互作用，孕母被动吸烟与家庭人均月收入低有正向相乘交互作用；Delta法、Wald法、PL法估计ICR 95%CI＞0、AP 95%CI＞0，S 95%CI＞1，孕母被动吸烟与家庭人均月收入低有协同相加交互作用，即孕母家庭人均月收入低同时存在被动吸烟者，其胎儿患先心的风险显著增加。ICR和S的意义相同，AP表示全部病例中可归因于两因素交互作用的病例所占比例，本例PL法中AP＝0.324，（95%CI：0.153～0.519），说明全部先心病例中归因于被动吸烟与家庭人均月收入低的相加交互作用所引起的病例占32.4%。

表 2 logistic相乘交互作用评价指标

表选项

图 1 相加交互作用示意图

图选项

表 3 logistic相加交互作用评价指标及其95%CI

表选项

本研究采用ICR等指标评价围孕期被动吸烟与家庭人均月收入既存在正向相乘交互作用又存在生物学上协同相加交互作用，表明Andersson法、Delta法、Wald法、PL法均能较好解释生物学相加作用。

统计学交互作用的线性模型一般为加法模型，乘积项表示有无相加交互作用；而广义线性属于乘法模型，乘积项表示有无相乘交互作用。生物学交互作用是多风险因素在发病的生物机制上的定性概念，两因素皆为病因前提下生物学机制的相互联系，包括协同和拮抗，生物学交互不同于统计模型中乘积项的分析。故相乘交互项OR值95%CI＜1为负相乘，OR值95%CI＞1为正相乘，相加交互项ICR值95%CI＞0为协同作用，＜0为拮抗作用。有时相乘与相加交互方向相悖，此时更宜运用生物学交互。在置信区间估计方法的选择上，Andersson等使用Delta法估计，显不够稳健易造成假阴性结果，如分析广东省先心监测网数据库中孕妇被动吸烟与围孕期饮酒致新生儿先心的关系，采用Delta法，则ICR＝7.635（95%CI：-1.777～17.046），而PL法为ICR＝7.635（95%CI：1.389～32.707），由于前者假阴性导致未正确评估出二者实际存在的协同效应。在风险比不对称时，应采用PL法，因Wald法的结果不稳定。

本研究可有效弥补其他非条件logistic模型的不足，稍加改写就能计算条件logistic和Cox比例风险模型交互作用参数和P值。Andersson等未尝试计算条件logistic的ICR值，需计算协方差矩阵再手工输入数据至Excel表格，在基因研究及大型临床流行病学课题中很容易遗漏数据或难以质控。SAS宏在大数据分析中重复操作性方面则优势明显。国外已将logistic模型交互作用分析纳入标准化流程中^[12]，本文宏程序较高效的运算了二分类条件/非条件logistic和Cox交互效应。源程序主要展示非条件logistic交互，临床资料结构复杂常常需要多种统计学方法建模，例如若想进行生存资料Cox交互作用和条件logistic分析，则仅需修改过程步骤以及协方差矩阵的相关位置参数即可便捷完成计算（宏程序见附录）。

附录：SAS 9.4宏程序

%macro additive2dum_LOGISTIC（dat，a，b，statues，cov1…covN）；

data P1；

set &dat；XX＝&a*&b；

PROC logistic OUTEST＝EST DATA＝P1 COVOUT；

/*生存资料Cox比例风险模型则修改为PROCPHREG；modeltime*&statues（0）＝&a &b XX且BETA3＝t（BETAP[，2：4]）；V3＝COV [1：3，2：4]；*/

class &A（ref＝first） &b（ref＝first）/param＝ref；

MODEL &statues（EVENT＝“1”）＝&a &b XX &cov1…&covN/ risklimits COVB include＝3 selection＝stepwise；

/*条件logistic则strata duizihao；stepwise后加ties＝discrete；且BETA3＝t（BETAP[，2：4]）；V3＝COV [1：3，2：4]；*/

contrast “OR10” &a 1 &b 0 XX 0/estimate＝exp；

contrast “OR01” &a 0 &b 1 XX 0/estimate＝exp；

contrast “OR11” &a 1 &b 1 XX 1/estimate＝exp；

RUN；

PROC EXPORT DATA＝estOUTFILE＝ “D：\&b” DBMS＝excel replace；RUN；

%mend additive2dum_LOGISTIC；

%let ICR_LOGISTIC＝%str（

PROC PRINT DATA＝EST；RUN；

data parmscovb；set EST；drop _lnlike_；

if_type_＝‘PARMS’ then output parms；elseif_ type_＝‘COV’ then output covb；run；

PROC IML；USE parms；READ ALL var _all_ INTO BETAP；

close parms；print betap；use covb；read all var _all_ into cov；

close covb；print cov；

BETA3＝t（BETAP[，2：4]）；V3＝COV [2：4，2：4]；PRINT BETA3 V3；

C＝BETA3[1，]＋BETA3[2，]＋BETA3[3，]；D＝EXP（C）-1；

I＝EXP（BETA3[1，]）＋EXP（BETA3[2，]）-1；F＝EXP（BETA3[1，]）＋EXP（BETA3[2，]）-2；H＝EXP（C）/D；

ICR＝EXP（C）-EXP（BETA3[1，]）-EXP（BETA3[2，]）＋1；AP＝ICR/EXP（C）；SI＝D/F；

A＝J（3，1，0）；A[1，]＝EXP（C）-EXP（BETA3[1，]）；A[2，]＝EXP（C）-EXP（BETA3[2，]）；

A[3，]＝EXP（C）；VARICR＝t（A）*V3*A；Z＝ICR/SQRT（VARICR）；

PVAL＝2*（1-PROBNORM（ABS（Z）））；L_ICR＝ICR-1.96* SQRT（VARICR）；

H_ICR＝ICR＋1.96*SQRT（VARICR）；B＝J（3，1，0）；

B[1，]＝（EXP（BETA3[2，]-1））/exp（C）；

B[2，]＝（EXP（BETA3[1，]-1））/exp（C）；

B[3，]＝I/exp（C）；VARAP＝t（B）*V3*B；

L_AP＝AP-1.96*SQRT（VARAP）；H_AP＝AP＋1.96*SQRT（VARAP）；S＝J（3，1，0）；

S[1，]＝H-EXP（BETA3[1，]）/F；S[2，]＝H-EXP（BETA3[2，]）/F；S[3，]＝H；

VARS＝t（S）*V3*S；L_LSI＝log（SI）-1.96*SQRT（VARS）；

H_LSI＝log（SI）＋1.96*SQRT（VARS）；L_SI＝exp（L_LSI）；

H_SI＝exp（H_LSI）；PRINT ICR L_ICR H_ICR AP L_AP H_AP SI L_SI H_SI Z PVAL；

run；quit；

option mcompilenote＝all；

%macro addit_sign（dat，a，b，y， beta0，beta1，beta2， beta3）；

PROC NLMIXED DATA＝&dat DF＝10000；

PARMS b0＝&beta0 b1＝&beta1 b2＝&beta2 b3＝&beta3；

odds＝exp（b0＋b1*&a＋b2*&b＋b3*&a*&b）；pi＝odds/（1＋odds）；MODEL &y～Binary（pi）；

estimate ‘ICR＝ICR using OR’ exp（b1＋b2＋b3）-exp（b1）-exp（b2）＋1；

estimate ‘AP using OR’（exp（b1＋b2＋b3）-exp（b1）-exp（b2）＋1）/（exp（b1＋b2＋b3））；

estimate ‘SI＝ICR using OR’（exp（b1＋b2＋b3）-1）/（exp（b1）＋exp（b2）-2）；run；

%mend addit_sign；option mcompilenote＝all；

%macro addit_sign1（dat，a，b，y，beta0，beta1，beta2，WaldICR）；

PROC NLP DATA＝&dat VARDEF＝N technique＝dbldogpcov；

PARMS b0＝&beta0，b1＝&beta1，b2＝&beta2，ICR＝&WaldICR；

b_interaction＝log（（ICR＋exp（b1）＋exp（b2）-1）/（exp（b1） * exp（b2）））；

if &y＝1 then p＝exp（b0＋（b1*&a）＋（b2*&b）＋（b_ interaction*&a*&b））/（1＋exp（b0＋（b1*&a）＋（b2*&b）＋（b_interaction * &a*&b）））；

else p＝1/（1＋exp（b0＋（b1*&a）＋（b2*&b）＋ （b_interaction * &a*&b）））；

loglike＝log（p）； max loglike ；PROFILE ICR /ALPHA＝0.05；RUN；

%mend；

%additive2dum_LOGISTIC（logistic，smoke，salary，y，cov1…covN）；&ICR_Logistic；

%addit_sign （logistic，smoke，salary，y， -0.1882， 0.3241， 0.2377， 0.3312）；

%addit_sign1 （logistic，smoke，salary，y， -0.1882， 0.3241， 0.2377， 0.7917）；