写在前面的话,此函数不适用于NHANES数据,也不能用于COX回归,请注意甄别。
在SCI文章中,交互效应表格(通常是表五)几乎是高分SCI必有。因为增加了亚组人群分析,增加了文章的可信度,能为文章锦上添花,增加文章的信服力,还能进行数据挖掘。
在既往文章《scitb5函数1.7版本(交互效应函数P for interaction)发布----用于一键生成交互效应表、森林图》中,本人发布了自己编写的scitb5函数,用于绘制交互效应表。效果还不错,旧版本不能生成亚组的P值和调节小数位数,新版本增加了以上功能,还增加了线性回归森林图的绘制,下面我来演示一下。
导入我们的早产数据和函数
source("E:/r/myfit/21scitb5.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)
定义协变量和分层因子
cov1<-c("lwt","smoke","ptl","ui","ftv","race")
Interaction<-c("race","smoke","ui")
生成结果
out<-scitb5(data=bc,x="age",y="low",Interaction=Interaction,cov = cov1,family="glm")
调节小数位数,dec调整OR的小数位数,pdec调整P值的小数位数,p.intervaue调整交互P值小数位数
out<-scitb5(data=bc,x="age",y="low",Interaction=Interaction,cov = cov1,family="glm",dec = 5,pdec = 5,p.intervaue = 5)
一键绘制森林图
plotforest(out)
加下来说一个包含的问题,在既往函数中,我要求协变量COV是要包含分层变量,不然就会报错,但也有些粉丝说这样有时候不够灵活,新版本中这个包含规则也是可以关掉的,我重新设置一下协变量h额分层
cov1<-c("lwt","smoke","ptl","ui","ftv")
Interaction<-c("race","smoke","ui")
以前这样是做不出来的,因为"race"没有在协变量里面
out<-scitb5(data=bc,x="age",y="low",Interaction=Interaction,cov = cov1,family="glm") #运行不了
现在可以了
out<-scitb5(data=bc,x="age",y="low",Interaction=Interaction,cov = cov1,family="glm",contain = F)
生成森林图
plotforest(out)
森林图可以调整演示和范围
plotforest(out,xticks=c(0,10,20),col = "red")
下面演示一下分组变量,这里ht是个分类数据,表示是否有高血压。
cov1<-c("lwt","smoke","ptl","ui","ftv","race")
Interaction<-c("race","smoke","ui")
out<-scitb5(data=bc,x="ht",y="low",Interaction=Interaction,cov = cov1,family="glm")
上图我要解释一下,做分类变量的时候需要设定一个参考,在race这个分层比较的时候,ht.0_race.1是被认定做参考的,什么意思呢?就是当ht等于0和race等于1这个亚组的人群是被默认为做参考比较的,其他组都是和它进行比较最后一组有个数据缺失,是因为这个分层一个数据都没有,所以就是NA。
所以说,分类变量进行亚组交互的时候,分类最好不要太多,要不数据会很大,而且有些层分不到数据。下面来演示一下粉丝的数据,主要是演示线性回归绘制森林图
我们先导入数据,数据我就不解释了
bc<-read.csv("E:/r/fensi/DII_GROUP.csv",sep=',',header=TRUE)
dput(names(bc))
bc[,c("sex_02", "nianling_02", "marriage_02", "GROUP","smoke_02")] <- lapply(bc[,c("sex_02", "nianling_02", "marriage_02","GROUP", "smoke_02")], factor)
str(bc)
########
Interaction<-c("sex_02", "nianling_02", "marriage_02", "smoke_02")
cov<-c("sex_02", "nianling_02", "marriage_02", "smoke_02","wc_02")
生成结果,这里的Y是连续变量,所以是线性回归
out<-scitb5(data=bc,x="GROUP",y="DII_total_02",Interaction=Interaction,cov = cov,family="glm")
旧版本是做不了线性回归森林图的,我看了一下有的文章确实有做线性回归森林图,所以在新版本也开发了这个线性回归森林图,方法也是一样的
plotforest(out)
我们注意一下上图有红箭头,有小数点存在,这和我的算法有关,你可以自己调整一下
plotforest(out,xticks=c(-3,0,5))
如果你还不会也没有关系,请看下面的视频操作
我暂时还没有把两个函数进行和并的打算,NHANES数据的COX回归亚组交互函数也在编写中,可能还要等一等。
需要获取scitb5.coxph 函数可以看下面这篇文章
scitb5函数2.1版本(交互效应函数P for interaction)发布----用于一键生成交互效应表、森林图