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 #高度
)
