R绘制logistic回归限制性立方样条图

限制性立方样条用于研究X和y的非线性关系，通常的结果是一条曲线。通过曲线可以判断转折点。

限制性立方样条包括线性模型、logistic模型、Cox回归模型三种情况，今天介绍基于logistic回归模型的限制性立方样条图。

变量要求：x为连续性变量，y为2分类变量。

今天展示的是BMI对心血管风险的影响，

自变量	因变量	协变量
BMI	心血管	地区、sex、age

老规矩，直接看数据形式和最终的效果图：

数据情况：这里展示前几行

结果图：

同样，注意两个##的地方代码就是需要更改的地方

.libPaths()#查看R包位置

##这里改成自己电脑的路径

setwd("C:/Users/12974/Desktop/百度经验/简书/R绘制限制性立方样条图")#设置工作空间

getwd()#加载工作空间

#包安装

install.packages("foreign")

install.packages("ggplot2")

install.packages("rms")

install.packages("survival")

install.packages("Hmisc")

install.packages("splines")

#批量包加载

ps <- c("foreign","ggplot2","rms","survival","Hmisc","splines")

for(i in ps){library(i, character.only = T)}; rm(i)

#导入数据

mydata <- read.csv("cc1.csv",as.is = TRUE,header = T,sep = ",", fileEncoding='utf-8')

names(mydata)#查看所有变量名字

attach(mydata)

#变量因子化，意思就是把分类变量变为真正的分类变量

##地区、sex需要

mydata$地区<-as.factor(mydata$地区)

mydata$sex<-as.factor(mydata$sex)

mydata$心血管<-as.factor(mydata$心血管)

##设置分类变量的参照组

mydata$地区<-relevel(mydata$地区, ref="1")

mydata$sex<-relevel(mydata$sex, ref="1")

mydata$心血管<-relevel(mydata$心血管, ref="0")

#注意：这里我说明一下，结局变量心血管：值是1/0；

#其实不设置R默认以“0“为参照

#接着为后续程序设定数据环境，也就是打包数据，这一步在预测模型中也常做

dd <- datadist(mydata)

options(datadist='dd')

#拟合logistic回归的限制性立方样条

##（心血管 为y，bmi 为x ，3是拟合曲线的时候采用三个节点）

##后面的+地区+sex+age 是一些协变量

fit<-lrm(心血管~rcs(bmi,3)+地区+sex+age,data=mydata)

#这里先看看结果

fit

##########图片1

#可以看到，地区展示了地区=2和地区=3的结果，因为前面我们设置地区参照组是 1

#性别展示sex=2的结果，因为前面我们设置sex参照组是 1

##由于限制性立方样条推荐拟合3~5个节点，这里分别拟合3个模型

fit3<-lrm(心血管~rcs(bmi,3)+地区+sex+age,data=mydata)

fit4<-lrm(心血管~rcs(bmi,4)+地区+sex+age,data=mydata)#模型出错

fit5<-lrm(心血管~rcs(bmi,5)+地区+sex+age,data=mydata)#模型出错

#然后查看AIC,选择AIC最小的

AIC(fit3)

AIC(fit4)

AIC(fit5)

#由于上面模型出现

Warning message:

  In lrm(心血管 ~ rcs(bmi, 5) + 地区 + sex + age, data = mydata) :

  Unable to fit model using “lrm.fit”

#这样的错误，后面我猜测应该是bmi对心血管影响只有1个转折点，所以4和5个节点无法拟合

#这个也是我猜的，有大神的话，评论区可以留言

#最后我选择fit3，3个节点

fit<-update(fit3)#更新模型

#使用anova()可以看p值,这里是为后续图上放 卡方值和P值 做准备

an<-anova(fit)

#这个结果主要是看bmi对心血管非线性关系的p值，可以发现二者存在非线性关系，如下图

an

########图片2

##生成预测值，并作图 （exp相当于把预测值转换成了OR）

plot(Predict(fit, bmi,fun=exp), anova=an, pval=T)

OR<-Predict(fit, bmi,fun=exp,ref.zero = TRUE)

#ggplot画图

ggplot(OR,anova=an, pval=T)

#进一步美化

#anova=an, pval=T：增加卡方值和P值

p1<-ggplot(anova=an, pval=T)+

  #画曲线

  geom_line(data=OR, aes(bmi,yhat),linetype=1,size=1,alpha = 0.9,colour="red")+

  #画置信区间

  geom_ribbon(data=OR, aes(bmi,ymin = lower, ymax = upper),alpha = 0.3,fill="red")+

  #x轴任意刻度：增加一条竖线

  geom_vline(aes(xintercept=23.9), colour="#BB0000", linetype="dashed")+

  #y轴任意刻度：增加一条横线

  geom_hline(yintercept=1, linetype=2,size=1)+

  #去除背景

  theme_classic()+

  #增加标签

  labs(title = "RCS", x="bmi", y="OR (95%CI)")+

  #x轴范围

  scale_x_continuous(limits = c(10, 50),

                     #x轴刻度

                   breaks = c(10,15,20,25,30,35,40,45,50))+

  #y轴范围

  scale_y_continuous(limits = c(0, 1.5),

                     #y轴刻度

                     breaks = c(0,0.5,1,1.5,2,2.5,3))+

  #手动给图上增加标签

  geom_text(aes(x=15,y=1,label='bmi=23.9'),

            vjust=1.5,hjust=0,size=2.5)

p1

########结果图1

#结果解释：

随着bmi增加，心血管风险增加，当BMI>23.9,心血管风险变化不大；

由于是模拟的数据，可能和实际不符，这里仅供参考

#那么接下来，绘制性别分层的图

##计算不同性别的OR值

OR1 <- Predict(fit, bmi, sex=c('1','2'),

               fun=exp,type="predictions",

               ref.zero=TRUE,conf.int = 0.95,digits=2)

#美化

p2<-ggplot()+

  #画曲线，多color = sex

  geom_line(data=OR1, aes(bmi,yhat, color = sex),

            linetype="solid",size=1,alpha = 0.9)+

  #画置信区间，多color = sex

  geom_ribbon(data=OR1,

              aes(bmi,ymin = lower, ymax = upper,fill = sex),

              alpha = 0.2)+

  #两条线的颜色

  scale_color_manual(values = c('red','blue'))+

  #两个置信区间的颜色

  scale_fill_manual(values = c("red","blue"))+

  #x轴任意刻度：增加一条竖线

  geom_vline(aes(xintercept=23.9), colour="#BB0000", linetype="dashed")+

  #x轴任意刻度：再增加一条竖线

  geom_vline(aes(xintercept=26), colour="#BB0000", linetype="dashed")+

  #y轴任意刻度：增加一条横线

  geom_hline(yintercept=1, linetype=2,size=1)+

  #去除背景

  theme_classic()+

  #增加标签

  labs(title = "RCS", x="bmi", y="OR (95%CI)")+

  #x轴范围

  scale_x_continuous(limits = c(10, 50),

                     #x轴刻度

                     breaks = c(10,15,20,25,30,35,40,45,50))+

  #y轴范围

  scale_y_continuous(limits = c(0, 1.5),

                     #y轴刻度

                     breaks = c(0,0.5,1,1.5,2,2.5,3))+

  #手动给图上增加标签

  geom_text(aes(x=15,y=1,label='bmi=23.9'),

            vjust=1.5,hjust=0,size=2.5)+

  #手动给图上增加标签

  geom_text(aes(x=27,y=0.9,label='bmi=23.9'),

            vjust=1.5,hjust=0,size=2.5)

p2

########结果图2

#可以看到男性和女性的结果与前面基本一致

#最后导出图片

ggsave(filename = "结果1.png",#命名

       plot=p1,#哪张图

       path = "C:/Users/12974/Desktop/百度经验/简书/R绘制限制性立方样条图",

       #保存路径

       units="px",

       width = 1200,#宽度

       height = 800 #高度

)

ggsave(filename = "结果2.png",#命名

       plot=p2,#哪张图

       path = "C:/Users/12974/Desktop/百度经验/简书/R绘制限制性立方样条图",

       #保存路径

       units="px",

       width = 1200,#宽度

       height = 800 #高度

)