R语言手把手教你制作一个交互效应表

交互作用效应(p for Interaction)在SCI文章中可以算是一个必杀技，几乎在高分的SCI中必出现，因为把人群分为亚组后再进行统计可以增强文章结果的可靠性.


今天我们使用R语言来绘制上图这样一个表格，继续使用的是我们的早产数据，我们先导入数据

bc<-read.csv("E:/r/test/zaochan.csv",sep=',',header=TRUE)
bc <- na.omit(bc)

这是一个关于早产低体重儿的数据（公众号回复：早产数据，可以获得该数据），低于2500g被认为是低体重儿。数据解释如下：low 是否是小于2500g早产低体重儿，age 母亲的年龄，lwt 末次月经体重，race 种族，smoke 孕期抽烟，ptl 早产史（计数），ht 有高血压病史，ui 子宫过敏，ftv 早孕时看医生的次数
bwt 新生儿体重数值。
我们先把分类变量转成因子

bc$race<-ifelse(bc$race=="black",1,ifelse(bc$race=="white",2,3))
bc$smoke<-ifelse(bc$smoke=="nonsmoker",0,1)
bc$low<-factor(bc$low)
bc$race<-factor(bc$race)
bc$ht<-factor(bc$ht)
bc$ui<-factor(bc$ui)

假设我们想研究母亲年龄和早产关系的影响，我们需了解不同种族间年龄和早产的交互影响。
先要把年龄分成3组数据

dat1<-subset(bc,bc$race==1)
dat2<-subset(bc,bc$race==2)
dat3<-subset(bc,bc$race==3)

分别使用这3个数据建立模型，注意一下模型是没有race这个指标的

f1<-glm(low~age + lwt + smoke + ptl + ui + ftv, data=dat1,
    family=binomial(link = "logit"))
f2<-glm(low~age + lwt + smoke + ptl + ui + ftv, data=dat2,
        family=binomial(link = "logit"))
f3<-glm(low~age + lwt + smoke + ptl + ui + ftv, data=dat3,
        family=binomial(link = "logit"))

接下来需要对每个模型提取系数，我们先对模型f1进行提取

b<-summary(f1)$coeff[2,1]#提取系数
se<-summary(f1)$coeff[2,2]#提取标准误
OR<-round(exp(b),2)
ll<-round(exp(b-1.96*se),2)
ul<-round(exp(b+1.96*se),2)

接下来需要对其他模型进行一样的提取系数，我这里写成一个循环来跑，先设定一个输出函数

jj<-function(data) {
  dat<-data
  fit<-glm(low~age + lwt + smoke + ptl + ui + ftv, data=dat,
           family=binomial(link = "logit"))
  b<-summary(fit)$coeff[2,1]#提取系数
  se<-summary(fit)$coeff[2,2]#提取标准误
  OR<-round(exp(b),3)
  ll<-round(exp(b-1.96*se),3)
  ul<-round(exp(b+1.96*se),3)
  d<-data.frame(OR,ll,ul)
}

把3个数据整成一个列表

dt<-list(dat1,dat2,dat3)

使用sapply函数来跑

out<-sapply(dt,jj)

需要对生成结果转置一下

out1<-as.data.frame(t(out))

这样结果OR值和可信区间就生成了，交互效应的P值需要另外计算，我这里就不演示了，有兴趣的看我既往文章，我这里P值是0.7358

p<-0.7358

接下来设置一下表格，现增加一行空值

out2<-rbind(c("","",""),out1)

放入我们的结果并重新命名

out2$p.for.Interaction<-c(0.075,"","","")
rownames(out2)<-c("race","black","white","other")

最后结果展示