新版亚组交互效应函数（P for interaction）newscitb5 1.3版本发布--用于一键生成交互效应表

在SCI文章中，交互效应表格（通常是表五）能为文章锦上添花，增加文章的信服力，增加结果的可信程度，还能进行数据挖掘。什么是亚组，通常就是特殊类型人群，比如男女，种族等，就是说你的数据放入特殊人群中结果还可靠吗？如果在各个特殊人群中，你的结果很稳定，说明你的结论很可靠。如果亚组的结论和你的数据数据结论相反，你可以拿来做个新论题。还可以比较不同亚组之间有无区别，比如做了心脏支架和没做支架的区别，可以发现很多新思路，易于数据挖掘。

目前不少粉丝私信问我，上面这种亚组交互的表格怎么做，
在既往咱们已经在文章《scitb5函数2.1版本（交互效应函数P for interaction）发布----用于一键生成交互效应表、森林图》中已经介绍了可以生成亚组交互表格的scitb5函数，但是在X为分类变量的时候，生成的亚组交互函数是如下图所示这样的

这两个表格是有区别的，主要是参照物不同，就拿第一个表格的SEX这栏来说，第一个表格主要是两个亚组男和女的情况，在男的这个亚组中，Q2-4是和Q1进行比较，女的亚组中也是一样，Q2-4是和Q1进行比较。因此会有两个Ref，因为分成了两个层。
在第二个表格进行了另外一种方式对数据的处理，等于是把数据分层了8个亚组，女性Q1组，女性Q2组，女性Q3组，女性Q4组，男性Q1组，男性Q2组，男性Q3组，男性Q4组，以女性女性Q1组为参考，其他数据和它进行了参考。
两种方法各有所长，都有文章进行发表，既往的scitb5函数生成的是第二种表格。
应粉丝的要求，我又编写了新的亚组交互函数用于生成下图也就是第一张图的这种表格，就叫newscitb5函数，想要取个响亮的名字真费劲

下面我来给大家演示一下，这种表格我们Y只能是二分类的，X可以是连续的和分类的，但最后都按分类的来处理。
咱们先导入数据和我写的函数，需要survival包和lmtest支持，这两个包要先安装

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)
bc$smoke<-factor(bc$smoke)

导入我写的函数，导入成功的话应该左下角出现13个函数

source("E:/r/myfit/newscisub13.R")

接下来咱们需要定义协变量和分层因子，cov1就是你模型中的协变量，Interaction就是你的分层变量

cov1<-c("age")
Interaction<-c("ht","ui","smoke")

生成数据，data就是你的数据，x就是你研究的变量，y就是你的结局变量，结局变量必须是分类变量，Interaction就是你的分层变量，cov就是你的协变量，family="glm"这里是固定的

out<-newscitb5(data=bc,x="race",y="low",Interaction=Interaction,cov = cov1,family="glm")

这样表格就可以轻松生成了。下面咱们来手动检验一下咱们的结果对不对，以ht为例子，如果你购买了既往的scitb5函数，也可以这样做出上面表格
先把数据分层ht的两个亚组数据

dat1<-subset(bc,bc$ht==0)
dat2<-subset(bc,bc$ht==1)

生成两个亚组的数据模型

fit1<-glm(low~race+age,family = "binomial",data = dat1)
fit2<-glm(low~race+age,family = "binomial",data = dat2)

生成两个模型结果

咱们可以看到手动手动算出来和咱们的表格一模一样，如果你购买了是scitb5，再把做出来的交互值填进去就做好啦。

下面咱们来做一个下图这样带有 p for trend的

做出这样p for trend的必须有两个条件，第一个是我们观察的变量X必须是连续变量，2是转换后的分类必须大于3个等级，就是像上图一样必须有G0,G1和G2一共3个，只有2个是做不出来趋势的。你看第一张图的时候它的X变量也是一个连续变量，做了分位数后再行亚组比较的

cov1<-c("ui")
Interaction<-c("ht","race","smoke")
out<-newscitb5(data=bc,x="age",y="low",Interaction=Interaction,cov = cov1,family="glm",contain=F)

出现下面红字并不是报错，而是提醒你观察变量X是个连续变量

有些细心的粉丝可能发现了ht等于1的时候有个p for trend的值没有算出来，是缺失的，为什么会这样呢？下面我来手动算下给大家看看，方法是参照了文章《手把手教你R语言计算校正混杂因素后的P for trend》，我这里就不解释了，直接算了
咱们先要把年龄进行分段然后进行标识

fenwei<-quantile(bc$age,c(0,.25,.50,.75,1))
Q<-cut(bc$age,breaks=c(fenwei),right = F,labels = c(0,1,2,3))#把 age 划分为 4 个等分到区间

但是传统方法算出来有个小缺点，它的最右侧不算进去，所以咱们可以看到有一个缺失，这里我使用我自编小程序转换一下

Q<-rankvar(bc$age,4)
bc$Q<-Q

大家看下我这里是没有缺失值的，其他都一样，然后和刚才一样把数据分组

dat1<-subset(bc,bc$ht==0)
dat2<-subset(bc,bc$ht==1)

再dat1中把中位数加进去

library(dplyr)
dat1<-dat1 %>% group_by(Q) %>%mutate(Median= median(age,na.rm=TRUE))%>%ungroup()

上图可以看到中位数已经生成，带入模型就可以了

f<-glm(low~Median+ui,family = "binomial",data = dat1)
summary(f)

最后咱们可以看到P值是0.361和咱们计算的一模一样，如果你是手动计算有缺失值的话可能会和我算的差一点，但是也应该差不多
接下来咱们来看看为什么 ht==1这组跑不出来，主要再数据结构上ui这里只有一个分类

我们可以看到只有一个分类是组成不了模型的

dat2<-dat2 %>% group_by(Q) %>%mutate(Median= median(age,na.rm=TRUE))%>%ungroup()
f<-glm(low~Median+ui,family = "binomial",data = dat2)

所以有时候不是数据多少的问题，结构也很重要。接着这个问题，如果使用函数过程中遇到这个报错
Error in getpfi(xlv, xxname, yvname, Interaction, flmadj, tmp0, xvname_, :
找不到对象’mdli’

也是一模一样的情况，你的代码没错，函数也没错，就是你的交互变量的亚组数据太少，生成不了模型，处理就是：可以删除莫格数据过少的亚组变量，或者增加数据量。

最后介绍一下怎么调节分位数回归的个数，默认是4个，我换个大点的数据来演示一下

library(foreign)
bc <- read.spss("E:/r/Breast cancer survival agec.sav",use.value.labels=F, to.data.frame=T)
bc <- na.omit(bc)

我们先来看看数据（公众号回复：乳腺癌，可以获得数据）：
age表示年龄，pathsize表示病理肿瘤大小（厘米），lnpos表示腋窝淋巴结阳性，histgrad表示病理组织学等级，er表示雌激素受体状态，pr表示孕激素受体状态，status结局事件是否死亡，pathscat表示病理肿瘤大小类别（分组变量），ln_yesno表示是否有淋巴结肿大，time是生存时间，后面的agec是我们自己设定的，不用管它。
分类变量转成因子

bc <- na.omit(bc)
bc$er<-as.factor(bc$er)
bc$pr<-as.factor(bc$pr)
bc$ln_yesno<-as.factor(bc$ln_yesno)
bc$histgrad<-as.factor(bc$histgrad)
bc$pathscat<-as.factor(bc$pathscat)

定义分层变量和协变量

cov1<-c("er")
Interaction<-c("ln_yesno","histgrad")

生成数据，默认是4个分组

out<-newscitb5(data=bc,x="age",y="status",Interaction=Interaction,cov = cov1,family="glm")

想要调节成5个只需调节num参数

out<-newscitb5(data=bc,x="age",y="status",Interaction=Interaction,cov = cov1,family="glm",num=5)

这样就轻松生成了。这个过程还生成了每个分位数的P值，但是不知道怎么做表格，所以没有导出来。COX回归版本正在写，还要等一等。

如果还是不会可以看下面视频操作。

有问题可以私信我。获取新版亚组交互函数，可以看下面这篇文章

新版亚组交互效应函数（P for interaction）newscitb5 1.3版本发布–用于一键生成交互效应表