如何建立非线性回归预测模型

专题合集更多教程

作者:李健民;审稿:龚志忠

 

在统计分析中,根据变量的不同类型可以建立不同的预测模型,如果因变量是连续型变量,最常见的是建立线性回归模型。但是,建立线性回归模型有很多前提条件(可以参考:SPSS操作:简单线性回归)。

 

由于实际的临床研究中,变量之间关系复杂,因变量和自变量之间并非呈现线性关系,如果强行建立线性回归模型,就会影响模型的预测准确性。对于此类数据,应该如何处理呢?之前医咖会发布过的《R语言课程》,王九谊老师在“多项式回归、分段回归、限制性立方样条...”视频课程中已做了详细介绍。本文以临床医生的角度,通过案例分析,结合R软件来讲解如何建立非线性回归模型,也对之前的视频教程内容作了延伸。

 

案例说明(模拟数据)

 

临床中心衰、肝硬化的病人,常伴有体液潴留和低钠血症,医生会选择使用托伐普坦进行超滤治疗,但是目前这个药物价格昂贵,未能广泛使用。

 

假设有一种新的利尿剂上市,价格便宜,且具有类似作用。为了探究新利尿剂的治疗效果,研究人员开展了一项临床试验,共入组149人(数据库名称为urinetest),因变量为患者每日尿量(变量名为urine),自变量为每日新利尿剂使用剂量(变量名为dosage)。

 

研究目的是为两者建立最合适的回归模型,分析步骤如下:

 

1、初步探索数据

2、建立简单线性回归

3、建立曲线方程

4、建立分段回归

5、建立样条回归

6、构建局部加权回归

7、建立广义可加模型

8、总结

 

分析步骤

 

分析数据前的准备工作

 

1、点击impordataset导入数据urinetest 

 

2、数据预览,View(urinetest)

 

 

3、加载相关的包,请加载前用install.packages()命令安装好

 

library(ggplot2)

library(segmented)

library(splines)

library(Hmisc)

library(rms)

library(mgcv)

library(caret)

 

一、数据探索

 

ggplot(urinetest, aes(dosage, urine) )+geom_point()#绘制散点图

 

从图形可以看出,当利尿剂使用剂量<25ml时,病人的尿量在2000-2300ml之间波动。当利尿剂剂量为25-30ml时,两者成线性关系。当>30ml时,随着利尿剂剂量的增加,尿量不再出现明显的变化。

 

由此看出,两者呈现出一种非线性的变化关系,存在阈值效应和饱和效应,在不同药物剂量范围内,剂量-反应关系函数差别很大,如果强行用单一的线性回归来建立预测建模,不符合临床实际,模型预测的准确性将会大打折扣。下面我们先用线性回归来分析一下。

 

二、建立线性回归模型

 

model.lm <- lm(urine ~ dosage, data = urinetest)#构建线性回归模型

 

summary(model.lm)#查看回归模型结果

 

模型结果如下:

 

(1)残差的最大值、最小值、中位数等,描述的是预测值和实际值之差的分布;

 

(2)回归方程的系数和统计学检验结果;

 

(3)模型的拟合情况。其中Residual standard error为残差标准误,是模型用自变量预测因变量的平均误差,该值越小说明模型拟合越好;Adjusted R-squared为调整R2,可理解为模型对数据集的解释程度,该值越大模型拟合程度越好。本研究中线性回归模型的残差标准误的值为159.8;调整R2为0.5902。

 

接下来看看线性回归的拟合效果

 

ggplot(urinetest, aes(dosage, urine) ) +  

  geom_point() +  

  stat_smooth(method = lm, formula = y ~ x)

 

从图形可以直观看出拟合直线与数据点存在一定的偏离,拟合效果不佳。

 

三、建立曲线方程

 

下面尝试用曲线模型去拟合,例如对数曲线型、指数曲线型、S曲线型等。我们以对数曲线为例。

 

model.log<- lm(urine ~ log(dosage), data = urinetest)#建立对数曲线方程

 

summary(model.log)#查看模型概况

 

对数曲线模型的残差标准误的值为151.5,调整R2为0.6318,两个指标比简单线性回归模型略有提高。

 

#拟合曲线

ggplot(urinetest, aes(dosage, urine)) +  

  geom_point() +  

  stat_smooth(method = lm, formula = y ~ log(x))

 

从图形可以看出,拟合曲线的效果较直线有所改善。

 

四、建立分段回归模型

 

在数据探索时我们发现,药物剂量和尿量的散点图分布呈现三段式变化特征,我们以此为依据,建立一个分段回归模型。在R中我们可以使用segmented这个包。

 

model.segmented <- segmented(model.lm)#构建分段回归模型

 

summary(model.segmented)#查看模型概况

 

分段回归结果显示,软件自动将模型分成了两段,拐点为dosage=32.534,残差标准误为124.9,调整R2为0.7499,两个指标较曲线模型得到了进一步提升。

 

#查看拟合效果

plot(dosage,urine, pch=1, cex=1.5)

abline(a=coef(model.lm)[1],b=coef(model.lm)[2],col="red",lwd= 2.5)

plot(model.segmented, col='blue', lwd= 2.5 ,add=T)

 

 

在构建的上述模型中,函数自动将模型分成了两段。但根据对散点图的分析,我们认为将模型分为三段更为合适,此时可以手动设置25和30两个剂量拐点,软件会自动寻找附近的点做为最佳拐点。

 

#手动设置拐点,分三段回归

 

model.segmented2 <- segmented(model.lm,psi=c(25,30))#构建分段回归模型

 

summary(model.segmented2)#查看模型概况

 

软件找到的两个最佳拐点分别为24.075和30.166,此时分段回归模型的残差标准误为99.01,调整R2为0.8427,预测效果比曲线模型明显提升。

 

#查看拟合效果

 

plot(dosage,urine, pch=1, cex=1.5)

abline(a=coef(model.lm)[1],b=coef(model.lm)[2],col="red",lwd= 2.5)

plot(model.segmented2, col='blue', lwd= 2.5 ,add=T)


 

五、样条回归

 

上述提到的曲线方程和分段回归两种方法都有一定的缺点。曲线方程是非局部的,当某一个因变量的值发生变化时,即使距离很远的点也会受到影响。如果采用多项式建立曲线方程,当多项式的幂较高时,自变量的一个微小变化,就会引起因变量很大的变化,得出的模型不适合外推到其他数据样本。而在分段回归模型中,每一段都是基于线性回归而建立的,拐点之间的连接显得比较生硬。

 

那么有没有办法建立一个既具有分段回归模型的优点,又可以拟合比较平滑的模型呢?样条回归则兼具曲线方程和分段回归的优点,可以灵活的分段展示自变量与因变量之间的关系。样条回归把数据集划分成一个个连续的区间,划分的点称为节点,每个节点之间用单独的模型(线性函数或者低阶多项式函数)来拟合。节点越多,模型就越灵活。但是过多的节点也会导致过拟合问题,所以一般先尝试设置3个节点为宜。

 

样条回归很多种,我们主要讲限制性立方样条回

 

model.spline <- lm(urine ~ rcs(dosage, c(20,30,35)))#建立样条回归,设置3个节点

 

summary(model.spline)#查看模型概况

 

样条回归模型的残差标准误为139.6,调整R2为0.6872。比线性回归和曲线回归好,但不如分段回归。

 

#样条回归拟合效果

 

ggplot(urinetest, aes(dosage, urine) ) +  

  geom_point() +  

  stat_smooth(method = lm, formula = y ~ rcs(x, c(20,30,35)) )  

 

六、Lowess函数建立局部加权回归

 

以上介绍的模型都是参数模型,选择什么样的曲线,设置多少个拐点,这些步骤都需要进行尝试,但也会容易出现过拟合现象。于是有学者提出了Lowess非参数回归,它没有回归系数可估计,只是在寻找一条拟合效果相对更好的曲线。

 

model.lowess<- lowess(urine ~ dosage)#建立局部加权回归

 

summary(model.lowess)#查看概况

 

#查看拟合

 

ggplot(urinetest, aes(dosage, urine)) +  

  geom_point() +  

  stat_smooth()

 

局部加权给一般只做数据探索,stat_smooth()就是默认用lowess画拟合图

 

七、广义可加模型

 

和lowess函数一样,广义可加模型也无法给出明确的系数,但它的适用范围更广,可以执行因变量与多个自变量之间的各种非参数拟合。

 

它可以是任意的单变量函数的叠加,这些函数既可以是线性,也可以是非线性。它的因变量可以服从二项分布、Poisson分布、Gamma分布等更广义的范畴。它的任务就是根据目前的数据,找出一条最贴合的曲线。

 

model.gam <- gam(urine ~ s(dosage), data = urinetest)#建立gam模型

 

summary(model.gam)#查看模型概况

 

广义可加模型的调整R2为0.837,但没有给出残差标准误的结果,所以我们需要利用模型生成预测值,用预测值和真实值进行比较,得出残差标准误为98.5,是上述众多模型中表现最优秀的。

 

pr.gam <- predict(model.gam,urinetest)#生成预测值

 

#计算RSME和R方

 

data.frame(RMSE = RMSE(pr.gam, urinetest$urine),

 

R2 = R2(pr.gam, urinetest$urine))

 

#查看模型拟合情况

 

ggplot(urinetest, aes(dosage, urine) ) +  

  geom_point() +  

  stat_smooth(method = gam, formula = y ~ s(x))

 

从图形可以看出,广义可加模型的曲线拟合效果非常好。虽然模型在本数据集中表现良好,但仍需要注意过拟合的情况。

 

八、总结

 

各个模型的拟合指标比较

 

 

通过比较模型指标,虽然广义可加模型表现较好,可是它并不能提供系数,无法解释变量之间的内在联系。而结合了专业背景而建立的分段回归模型表现相对更为优异。

 

分段回归的结果汇总如下:

 

 

为方便临床应用,我们将剂量节点取整,分别为A=24和C=30,分段回归方程的书写格式为:

 

Y=β01X12(X1-A)X23(X1-C)X3

 

(1)当X1≤A时,即X≤24时,X2=0,X3=0,

Y=β01X1=2168.470-0.692*X

 

(2)当A<X1≤C,即24<X≤30时,X2=1,X3=0,

Y=β01X12(X1-A)X2=2168.470-0.692*X+92.703*(X-24)=-56.402+92.011*X

 

(3)当X1>C,即X>30时,X2=1,X3=1,

Y=β01X12(X1-A)X23(X1-C)X3

=2168.470-0.692*X+92.703*(X-24)-96.839*(X-30)=2848.768-4.828*X

 

相关内容推荐:

 

R视频课程:多项式回归、分段回归、限制性立方样条...

 

扫码关注“医咖会”公众号,及时获取最新统计教程!

描述问题
选择一个标签 (请选择一个与您问题最相符的标签)
提交问题
我要提问
描述问题
选择一个标签 (请选择一个与您问题最相符的标签)
提交问题
描述问题
选择一个标签 (请选择一个与您问题最相符的标签)
    提交问题