R 做交互作用分析本人菜鸟从网上找的，运行到最后没出现错误提示，可结果不知道在哪输出输出结果应该是交互作用的三个指标 RERI AP S，求指教。代码如下：
test<-read.table (file="C:/test2.dat",header=T,fill=TRUE)
name(test)
confounders<-cbind(test$dummyag1,test$dummyag2,test$dummyag3)
test$A<-test$ASE
test$B<-test$APTPN22
test$case<-test$KIND

interaction<-function(inter,case,A,B,confounders){
inter$A0B0[inter$A==0 & inter$B==0]<-1
inter$A0B0[inter$A==1 & inter$B==0 | inter$A==0 & inter$B==1 | inter$A==1 & inter$B==1]<-0
inter$A1B0[inter$A==1 & inter$B==0]<-1
inter$A1B0[inter$A==0 & inter$B==0 | inter$A==0 & inter$B==1 | inter$A==1 & inter$B==1]<-0
inter$A0B1[inter$A==0 & inter$B==1]<-1
inter$A0B1[inter$A==1 & inter$B==0 | inter$A==0 & inter$B==0 | inter$A==1 & inter$B==1]<-0
inter$A1B1[inter$A==1 & inter$B==1]<-1
inter$A1B1[inter$A==1 & inter$B==0 | inter$A==0 & inter$B==1 | inter$A==0 & inter$B==0]<-0
glmdata<-cbind("A1B1"=inter$A1B1,"A1B0"=inter$A1B0,"A0B1"=inter$A0B1,"A0B0"=inter$A0B0,
as.data.frame(confounders,row.names=NULL,optional=TRUE))

interaction.Iglm<-glm(case~.,family=binomial(link="logit"),
data=glmdata)
mod<-summary(interaction.Iglm)

cor<-vcor(interaction.Iglm)
A1B1var<-cor[2,2]
A1B0var<-cor[3,3]
A0B1var<-cor[4,4]
cov23<-cor[4,3]
cov13<-cor[4,2]
cov12<-cor[3,2]

x<-summary(interaction.Iglm)$coefficients
OR<-round(exp(c(x[,1])),digits=3)
OR.LowLi<-round(exp(c(x[,1]-1.96*x[,2])),digits=3)
OR.UppLi<-round(exp(c(x[,1]+1.96*x[,2])),digits=3)
outf<-data.frame(OR,OR.LowLi,OR.UppLi)
out<-data.matrix(outf)
ORii<-out[2,1]
ORio<-out[3,1]
ORoi<-out[4,1]
ORoo<-1.0
ORg<-c(ORoo,ORio,ORoi,ORii)
ORn<-c("OR_A0B0","OR_A1B0","OR_A0B1","OR_A1B1")

RERI<-round(out[2,1]-out[4,1]-out[3,1]+1,digit=3)
hr1<-(-out[3,1])
hr2<-(-out[4,1])
hr3<--out[2,1]
SeRERI<-sqrt(hr1**2*A1B0var+hr2**2*A0B1var+hr3**2*A1B1var+2*hr1*hr2*cov23+2*hr1*hr3*cov12+2*hr3*hr2*cov12)
RERI.95Low<-round(RERI-1.96*SeRERI,digits=3)
RERI.95High<-round(RERI+1.96*SeRERI,digits=3)

ha1<-(-(ORio/ORii))
ha2<-(-(ORoi/ORii))
ha3<-(out[3,1]+out[4,1]-1)/out[2,1]
SeAP<-sqrt(ha1**2*A1B0var+ha2**2*A0B1var+ha3**2*A1B1var+2*ha1*ha2*cov23+2*ha1*ha3*cov12+2*ha3*ha2*cov12)
AP<-round(RERI/ORii,digits=3)
AP.95Low<-round(AP-1.96*SeAP,digits=3)
AP.95High<-round(AP+1.96*SeAP,digits=3)

SI<-round((out[2,1]-1)/(out[4,1]+out[3,1]-2),digits=3)
hs1<-(-out[3,1]/(out[3,1]+out[4,1]-2))
hs2<-(-out[4,1]/(out[3,1]+out[4,1]-2))
hs3<-(-out[2,1]/(out[2,1]-1))
SeSI<-sqrt(hs1**2*A1B0var+hs2**2*A0B1var+hs3**2*A1B1var+2*hs1*hs2*cov23+2*hs1*hs3*cov12+2*hs3*hs2*cov12)
SI.95Low<-round(SI-1.96*SeSI,digits=3)
SI.95High<-round(SI+1.96*SeSI,digits=3)

IM<-data.frame(RERI,RERI.95Low,RERI.95High,AP,AP.95Low,AP.95High)
SyI<-data.frame(SI,SI.95Low,SI.95.High)
ILST<-list(barplot((ORg),names.arg=(ORn)),mod,out,SyI,IM)
return(ILST)
}