[程序分享] R 做交互作用分析有些疑问求指点 [推广有奖]

1关注
0粉丝

已卖：22份资源

博士生

52%

还不是VIP/贵宾

-

0%

威望: 0 级
论坛币: 200 个
通用积分: 3.8500
学术水平: 0 点
热心指数: 1 点
信用等级: 0 点
经验: 17834 点
帖子: 98
精华: 0
在线时间: 469 小时
注册时间: 2012-8-6
最后登录: 2024-12-3

楼主

怪峩、冭天眞 发表于 2014-11-29 09:28:41 |AI写论文

是否 +2 论坛币

k人参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

R软件做交互作用分析代码是从文献上弄的，可是运行到最后没有提示错误，可是结果在哪呢没有结果还是没输出出来....结果应该是交互作用的三个指标吧？？

菜鸟一个就会复制粘贴，代码如下求高手

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)
}