临床上,因变量和临床的结局有时候不是线性关系,而回归模型有一个重要的假设就是自变量和因变量呈线性关联,因此非线性关系模型用回归分析来拟合受到限制。因此,一个更好的解决方法是拟合自变量与因变量之间的非线性关系,限制性立方(Restricted cubic spline,RCS)就是分析非线性关系的最常见的方法之一。
既往教程中我们介绍了使用R语言在COX回归模型基础上绘制限制立方条图,后台有不少粉丝说道不会制作logistic回归和线性回归的限制立方条图,我们今天分别来讲讲。
先讲logistic回归绘制立方条图:
logistic回归主要用于结果是二分类变量,没有相关时间变量的模型。我们继续使用我们原来的乳腺癌数据凑合一下,先导入数据,
library(foreign)
library(rms)
be <- read.spss("E:/r/test/Breast cancer survival agec.sav",
use.value.labels=F, to.data.frame=T)
be <- na.omit(be)
我们先来看看数据:
age表示年龄,pathsize表示病理肿瘤大小(厘米),lnpos表示腋窝淋巴结阳性,histgrad表示病理组织学等级,er表示雌激素受体状态,pr表示孕激素受体状态,status结局事件是否死亡,pathscat表示病理肿瘤大小类别(分组变量),ln_yesno表示是否有淋巴结肿大,time是生存时间,后面的agec是我们自己设定的,不用管它。
转换分类变量并且抽取一部分变量等下来拟合模型,这回我们就不需要时间变量time了
attach(be)
be<-data.frame(age,status,ln_yesno)
be$ln_yesno<-as.factor(be$ln_yesno)
对数据进行打包,整理
dd <- datadist(be) #为后续程序设定数据环境
options(datadist='dd') #为后续程序设定数据环境
拟合模型
fit <-lrm(status ~ rcs(age, 4)+ln_yesno,data=be)
an<-anova(fit)
绘图
plot(Predict(fit, age,fun=exp), anova=an, pval=T)
后面的步骤其实和原教程基本一样,我这里就简单点,直接上代码了,需要详细点的可以看上一篇教程
dd$limits$age[2] <-50###设定标准
fit1=update(fit)###更新模型
OR<-Predict(fit, age,fun=exp,ref.zero = TRUE) ##生成预测值
ggplot()+geom_line(data=OR, aes(age,yhat),linetype=1,size=1,alpha = 0.9,colour="red")+
geom_ribbon(data=OR, aes(age,ymin = lower, ymax = upper),alpha = 0.3,fill="red")+
geom_hline(yintercept=1, linetype=2,size=1)+theme_classic()+
labs(title = "RCS", x="age", y="OR(95%CI)")###绘图
接下来讲怎么绘制线性回归模型,这里我们使用我们原来的臭氧的数据,主要描述的是臭氧浓度和大气一些相关指标的情况,因为有些数据是非线性的,使用Logistic回归不合适我们先导入数据看一下
library(foreign)
library(rms)
be <- read.spss("E:/r/test/ozone.sav",
use.value.labels=F, to.data.frame=T)
names(be)
数据中有七个变量,ozon每日臭氧水平为结局变量,Inversion base height(ibh)反转基准高度,Pressure gradient (mm Hg) 压力梯度(mm Hg),Visibility (miles) 能见度(英里),Temperature (degrees F) 温度(华氏度),Day of the year日期,vh我也不知道是什么,反正就是一参数,这里所有的变量都是连续的。
对数据进行打包,整理
dd <- datadist(be) #为后续程序设定数据环境
options(datadist='dd') #为后续程序设定数据环境
拟合模型
fit<-ols(ozon ~rcs(ibh, 4)+temp,data=be)
summary(fit)
an<-anova(fit)
生成预测值并且做图,我们要注意一下,这里的我们使用线性回归,所以我们是使用β,不再使用OR或者HR了,所以fun=exp这步不需要了
Predict(fit,ibh)
plot(Predict(fit,ibh),anova=an, pval=T)
也可以像上面教程一样加一条参考线,我这里就不加了文章来源:https://www.toymoban.com/news/detail-412976.html
OLS1<-Predict(fit, ibh)
ggplot()+geom_line(data=OLS1, aes(ibh,yhat),linetype=1,size=1,alpha = 0.9,colour="red")+
geom_ribbon(data=OLS1, aes(ibh,ymin = lower, ymax = upper),alpha = 0.3,fill="red")+
theme_classic()+
labs(title = "RCS", x="ibh", y="ozon")
由上图看出,臭氧浓度和高度关系并不是越高越浓,而是在一定高度内臭氧浓度最高。
更多精彩文章请关注公众号:零基础说科研
文章来源地址https://www.toymoban.com/news/detail-412976.html
到了这里,关于R语言绘制限制立方条图2(基于logistic回归和线性回归)的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!