既往我们对logistic回归交互项（交互作用）进行了可视化分析，很多人对交互作用影响还不是很明白，都在后台发私信询问。其实交互作用就是一个变量通过对另一个变量的影响，导致结局变量改变,最常见的就是我们药物之间的交互作用，比如β受体阻滞剂对血糖有影响，同时和降糖药使用可能会加强或者减弱降糖药的降糖效果，这就是交互作用。
通常分类和连续变量之间的交互作用可以解释为分类变量斜率的变化，
继续使用我们的不孕症数据（公众号回复：不孕症可以获得该数据）来演示。
先导入数据

library(rms)
bc<-read.csv("E:/r/test/buyunzheng.csv",sep=',',header=TRUE)

数据有8个指标，最后两个是PSM匹配结果，我们不用理他，其余六个为：
Education：教育程度，age:年龄，parity产次，induced：人流次数，case：是否不孕，这是结局指标，spontaneous：自然流产次数。
有一些变量是分类变量，我们需要把它转换一下

bc$education<-ifelse(bc$education=="0-5yrs",0,ifelse(bc$education=="6-11yrs",1,2))
bc$spontaneous<-as.factor(bc$spontaneous)
bc$case<-as.factor(bc$case)
bc$induced<-as.factor(bc$induced)
bc$education<-as.factor(bc$education)

建立模型，假如我们想知道流产和年龄这两个指标是不是存在交互影响

f1<- glm(case ~ age+spontaneous+age:spontaneous,
         family = binomial(link = logit), data = bc)
summary(f1)

通过拟合模型我们可以把方程式写出来，我们在这里使用冒号（：）来表示交互项
：
Y=-0.92912+age×-0.01455
+（-0.70606×spontaneous1）
+（-0.45443×spontaneous1）
+（0.06040×age:spontaneous1）
+（0.08028×age:spontaneous2）
我们得出方程后就可以根据不同类型的变量，写出拟合模型
对于没有流产的患者:
Y=-0.92912+age×-0.01455
对于流产1次的患者
Y=-0.92912+（age×-0.01455）-1.63517989 +0.06040×age
= -1.635 +0.04584×age
对于流产2次以上的患者
Y=-0.92912+（age×-0.01455）-0.45443 + 0.08028×age
=-1.384+0.0657×age
由此我们得出了每个类型流产对应的不孕症患者的Y的系数，我们也可以通过R来计算这些参数

dia<-coef(f1)
dia

没有流产的患者

small<-c(dia[1],dia[2])
small

流产1次的患者

med<-c(dia[1]+dia[3],dia[2]+dia[5])
med

流产2次以上的患者

lar<-c(dia[1]+dia[4],dia[2]+dia[6])
lar

得出了系数，我们就可以进行数据作图，显示交互作用
先生成预测值

bb<-predict(f1,newdata = bc,type = "link",se =T)

然后进行绘图

cols<-c("black","red","green")###设定颜色
plot(bb$fit~bc$age,col=cols[bc$spontaneous],cex=0.5) ###绘图
abline(coef = small,lwd=2,col="black")##添加线条
abline(coef = med,lwd=2,col="red")##添加线条
abline(coef = lar,lwd=2,col="green")##添加线条
legend("topright",legend = c("spontaneous0","spontaneous1","spontaneous2"),
       lty = 1,lwd = 3,col = cols) ##图例

图片还可以进一步修饰，如调整线条和图例位置等，我这里就不弄了，如图所示随着年龄增长，自然流产2次以上的患者不孕的概率比没有自然流产的患者明显增高。

交互作用效应(p for Interaction)在SCI文章中可以算是一个必杀技，几乎在高分的SCI中必出现，因为把人群分为亚组后再进行统计可以增强文章结果的可靠性，进行可视化后可以清晰的表明变量之间的关系。不仅如此，交互作用还可以使用来进行数据挖掘。在既往文章中，我们已经介绍了怎么使用R语言和SPSS对logistic回归亚组交互效应（交互作用）进行可视化分析（见下图）


后台有粉丝问能不能进行COX回归交互效应可视化分析和对连续变量进行交互效应可视化分析，都是可以的，下面我们来一一进行演示。
COX回归交互效应可视化分析和logistic回归交互项的可视化分析的步骤基本一样的，继续使用我们的乳腺癌数据（公众号回复：乳腺癌，可以获得数据）我们先导入数据和R包

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

我们先来看看数据：
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)

建立cox回归方程，我们假设age和histgrad有交互

f1<- coxph(Surv(time, status) ~ age + histgrad+
           pathsize+ er + pr + ln_yesno+age*histgrad, data=bc)
summary(f1)

这里发现age：histgrad的P值大于0.05，两者之间没有交互效应，我们不管他，这里主要是演示怎么做图

visreg(f1, "age", "histgrad", ylab="log(Hazard ratio)")

可以把3张图合并起来

visreg(f1, "age", "histgrad", ylab="log(Hazard ratio)", overlay = T)

还可以进一步修改

plot(visreg(f1,xvar = "age",by="histgrad",plot=F),xlab="age",ylab="log(Hazard ratio)",
     overlay = T,partial = F,rug=F,,line=list(lty=1:6))
legend("bottomright",
       c("histgrad1","histgrad2","histgrad3"),
       lty=c(1,1,1),
       col=c("red","green","blue"),
       lwd=c(2,1,1),
       bty="n")

也可以按箱线图来表示

plot(visreg(f1, "histgrad", by="age",plot=F),xlab="age",ylab="log(Hazard ratio)",
     overlay = T,partial = F,rug=F,,line=list(lty=1:6))
legend("bottomright",
       c("age38","age55","age74"),
       lty=c(1,1,1),
       col=c("red","green","blue"),
       lwd=c(2,1,1),
       bty="n")

从图上中数据的可信区间严重重合也可以看出来两者并没有交互。
下面来演示连续变量的交互可视化分析，使用我们的臭氧数据（公众号回复：臭氧可以获得数据）
我们先导入R包和臭氧数据

library(splines)
be <- read.spss("E:/r/test/ozone.sav",
                use.value.labels=F, to.data.frame=T)

数据中有七个变量，ozon每日臭氧水平为结局变量，Inversion base height（ibh）反转基准高度，Pressure gradient (mm Hg) 压力梯度（mm Hg），Visibility (miles) 能见度（英里），Temperature (degrees F) 温度（华氏度），Day of the year日期，vh我也不知道是什么，反正就是一参数，这里所有的变量都是连续的。
假设我们想知道反转基准高度和温度对臭氧浓度有无存在交互影响

fit <- lm(ozon ~ vh +ns(ibh, df=2)*ns(temp, df=2), data=be)
summary(fit)

visreg2d(fit, "ibh", "temp")

也可以自己更改想要的颜色

visreg2d(fit, "ibh", "temp", color.palette=colorRampPalette(c("black", "white", "purple")))

也可以绘制3D立体图

visreg2d(fit,"ibh", "temp", plot.type="persp")

也可以绘制成这种

plot(visreg(fit,xvar = "ibh",by="temp",plot=F),xlab="age",ylab="log(Probability)",
     overlay = T,partial = F,rug=F,,line=list(lty=1:6))
legend("topright",
       c("temp41","temp62","temp81"),
       lty=c(1,1,1),
       col=c("red","green","blue"),
       lwd=c(2,1,1),
       bty="n")

交互作用效应(p for Interaction)在SCI文章中可以算是一个必杀技，几乎在高分的SCI中必出现，因为把人群分为亚组后再进行统计可以增强文章结果的可靠性，不仅如此，交互作用还可以使用来进行数据挖掘。在既往文章中，我们已经介绍了怎么使用R语言可视化visreg包对交互作用进行可视化分析（见下图），

但是使用visreg包后我们对图形的细节掌控力不够，今天我们使用重新绘制手动图形，这在SCI论文中非常实用，几乎可以做出所有的交互可视化，也有利于我们了解制图原理。
废话不多说，我们实操一下。我们使用人流后导致不孕的数据集(关注公众号后回复：不孕症，可以获得数据)，我们先导入看一下

library(rms)
bc<-read.csv("E:/r/test/buyunzheng.csv",sep=',',header=TRUE)

数据有8个指标，最后两个是PSM匹配结果，我们不用理他，其余六个为：
Education：教育程度，age:年龄，parity产次，induced：人流次数，case：是否不孕，这是结局指标，spontaneous：自然流产次数。
有一些变量是分类变量，我们需要把它转换一下

bc$education<-ifelse(bc$education=="0-5yrs",0,ifelse(bc$education=="6-11yrs",1,2))
bc$spontaneous<-as.factor(bc$spontaneous)
bc$case<-as.factor(bc$case)
bc$induced<-as.factor(bc$induced)
bc$education<-as.factor(bc$education)

建立模型，假如我们想知道流产（包括人流和自然流产）和年龄这两个指标是不是存在交互影响

f1<- glm(case ~ age + education + parity + induced+spontaneous+age*induced*spontaneous,
           family = binomial(link = logit), data = bc)
summary(f1)

我们看到自然流产和年龄之间可能存在交互（下图）

重新建立模型

f1<- glm(case ~ age  + parity +spontaneous+age:spontaneous,
         family = binomial(link = logit), data = bc)

建立新数据集

attach(bc)
newdata2<-data.frame(age,spontaneous,parity=mean(parity))

把模型带入新数据集

newdata3<-data.frame(newdata2,predict(f1,newdata = newdata2,type = "link",se =T))

生成预测值

newdata4<-within(newdata3,{
  pre<-plogis(fit)
  ll<-plogis(fit-(1.96*se.fit))
  ul<-plogis(fit+(1.96*se.fit))
})

最后绘图

ggplot(newdata4,aes(x=age,y=pre))+
  geom_line(aes(col=spontaneous),size=1)

最后对图形进行修饰

最后得出和visreg包同样的结果，随着年龄增长，自然流产2次以上的患者不孕的概率比没有自然流产的患者明显增高。
我们还可以在细节和年龄范围进一步修饰，这里就不多说了，总之，数据交互可视化是一项非常实用的技能，不仅可以用于数据分析，还可以用于数据挖掘。