讯享网

 <p id="20PGNB4Q"><strong>作者：李健民；审稿：龚志忠</strong></p><p id="20PGNB4S">在统计分析中，根据变量的不同类型可以建立不同的预测模型，如果因变量是连续型变量，最常见的是建立线性回归模型。但是，建立线性回归模型有很多前提条件（可以参考：）。</p><p id="20PGNB4U">由于实际的临床研究中，变量之间关系复杂，因变量和自变量之间并非呈现线性关系，如果强行建立线性回归模型，就会影响模型的预测准确性。对于此类数据，应该如何处理呢？之前医咖会发布过的《R语言课程》，王九谊老师在“”视频课程中已做了详细介绍。本文以临床医生的角度，通过案例分析，结合R软件来讲解如何建立<strong>非线性回归模</strong><strong>型</strong>，也对之前的视频教程内容作了延伸。</p><p id="20PGNB52"><strong>案例说明（模拟数据）</strong></p><p id="20PGNB55">临床中心衰、肝硬化的病人，常伴有体液潴留和低钠血症，医生会选择使用托伐普坦进行超滤治疗，但是目前这个药物价格昂贵，未能广泛使用。</p><p id="20PGNB57">假设有一种新的利尿剂上市，价格便宜，且具有类似作用。为了探究新利尿剂的治疗效果，研究人员开展了一项临床试验，共入组149人（数据库名称为urinetest），因变量为患者每日尿量（变量名为urine），自变量为每日新利尿剂使用剂量（变量名为dosage）。</p><p id="20PGNB59">研究目的是为两者建立最合适的回归模型，<strong>分析步骤</strong>如下：</p><p id="20PGNB5B">1、初步探索数据</p><p id="20PGNB5C">2、建立简单线性回归</p><p id="20PGNB5D">3、建立曲线方程</p><p id="20PGNB5E">4、建立分段回归</p><p id="20PGNB5F">5、建立样条回归</p><p id="20PGNB5G">6、构建局部加权回归</p><p id="20PGNB5H">7、建立广义可加模型</p><p id="20PGNB5I">8、总结</p><p id="20PGNB5L"><strong>分析步骤</strong></p><p id="20PGNB5O">分析数据前的<strong>准备工作</strong></p><p id="20PGNB5Q">1、点击impordataset导入数据urinetest</p><p id="20PGNB5S">2、数据预览，View（urinetest）</p><p class="f_center"><img src="https://nimg.ws.126.net/?url=http%3A%2F%2Fdingyue.ws.126.net%2F2023%2F0810%2F66e20b25j00rz6ayx000adgg0089008g.jpg&thumbnail=660x&quality=80&type=jpg"/><br/></p><p id="20PGNB5V">3、加载相关的包，请加载前用install.packages()命令安装好</p><p id="20PGNB61">library(ggplot2)</p><p id="20PGNB62">library(segmented)</p><p id="20PGNB63">library(splines)</p><p id="20PGNB64">library(Hmisc)</p><p id="20PGNB65">library(rms)</p><p id="20PGNB66">library(mgcv)</p><p id="20PGNB67">library(caret)</p><p id="20PGNB6A"><strong>一、数据探索</strong></p><p id="20PGNB6D">ggplot(urinetest, aes(dosage, urine) )+geom_point()#绘制散点图<br/></p><p class="f_center"><img src="https://nimg.ws.126.net/?url=http%3A%2F%2Fdingyue.ws.126.net%2F2023%2F0810%2Fed8b43c8j00rz6ayy000ld200bj0070g00bj0070.jpg&thumbnail=660x&quality=80&type=jpg"/><br/></p><p id="20PGNB6G">从图形可以看出，当利尿剂使用剂量&lt;25ml时，病人的尿量在2000-2300ml之间波动。当利尿剂剂量为25-30ml时，两者成线性关系。当&gt;30ml时，随着利尿剂剂量的增加，尿量不再出现明显的变化。</p><p id="20PGNB6I">由此看出，两者呈现出一种<strong>非线性的变化关系</strong>，存在阈值效应和饱和效应，在不同药物剂量范围内，剂量-反应关系函数差别很大，如果强行用单一的线性回归来建立预测建模，不符合临床实际，模型预测的准确性将会大打折扣。下面我们先用线性回归来分析一下。</p><p id="20PGNB6L"><strong>二、建立线性回归模型</strong></p><p id="20PGNB6O">model.lm &lt;- lm(urine ~ dosage, data = urinetest)#构建线性回归模型</p><p id="20PGNB6Q">summary(model.lm)#查看回归模型结果</p><p class="f_center"><img src="https://nimg.ws.126.net/?url=http%3A%2F%2Fdingyue.ws.126.net%2F2023%2F0810%2Fb20d5355j00rz6ayy002sd200bj006bg00bj006b.jpg&thumbnail=660x&quality=80&type=jpg"/><br/></p><p id="20PGNB6T"><strong>模型结果</strong>如下：</p><p id="20PGNB6V">（1）残差的最大值、最小值、中位数等，描述的是预测值和实际值之差的分布；</p><p id="20PGNB71">（2）回归方程的系数和统计学检验结果；</p><p id="20PGNB73">（3）模型的拟合情况。其中Residual standard error为残差标准误，是模型用自变量预测因变量的平均误差，该值越小说明模型拟合越好；Adjusted R-squared为调整R2，可理解为模型对数据集的解释程度，该值越大模型拟合程度越好。本研究中线性回归模型的残差标准误的值为159.8；调整R2为0.5902。</p><p id="20PGNB74">接下来看看<strong>线性回归的拟合效果</strong></p><p id="20PGNB76">ggplot(urinetest, aes(dosage, urine) ) +</p><p id="20PGNB77">geom_point() +</p><p id="20PGNB78">stat_smooth(method = lm, formula = y ~ x)</p><p class="f_center"><img src="https://nimg.ws.126.net/?url=http%3A%2F%2Fdingyue.ws.126.net%2F2023%2F0810%2F242c2eabj00rz6ayz000vd200bj006zg00bj006z.jpg&thumbnail=660x&quality=80&type=jpg"/><br/></p><p id="20PGNB7B">从图形可以直观看出拟合直线与数据点存在一定的偏离，拟合效果不佳。</p><p id="20PGNB7E"><strong>三、建立曲线方程</strong></p><p id="20PGNB7H">下面尝试用曲线模型去拟合，例如对数曲线型、指数曲线型、S曲线型等。我们以对数曲线为例。</p><p id="20PGNB7J">model.log&lt;- lm(urine ~ log(dosage), data = urinetest)#建立对数曲线方程</p><p id="20PGNB7L">summary(model.log)#查看模型概况</p><p class="f_center"><img src="https://nimg.ws.126.net/?url=http%3A%2F%2Fdingyue.ws.126.net%2F2023%2F0810%2F144ab4b2j00rz6ayz002sd200b9006kg00b9006k.jpg&thumbnail=660x&quality=80&type=jpg"/><br/></p><p id="20PGNB7O">对数曲线模型的残差标准误的值为151.5，调整R2为0.6318，两个指标比简单线性回归模型略有提高。<br/></p><p id="20PGNB7Q">#拟合曲线</p><p id="20PGNB7R">ggplot(urinetest, aes(dosage, urine)) +</p><p id="20PGNB7S">geom_point() +</p><p id="20PGNB7T">stat_smooth(method = lm, formula = y ~ log(x))</p><p class="f_center"><img src="https://nimg.ws.126.net/?url=http%3A%2F%2Fdingyue.ws.126.net%2F2023%2F0810%2F3fa03057j00rz6az0000wd200bj0070g00bj0070.jpg&thumbnail=660x&quality=80&type=jpg"/><br/></p><p id="20PGNB80">从图形可以看出，拟合曲线的效果较直线有所改善。</p><p id="20PGNB83"><strong>四、建立分段回归模型</strong></p><p id="20PGNB86">在数据探索时我们发现，药物剂量和尿量的散点图分布呈现三段式变化特征，我们以此为依据，建立一个分段回归模型。在R中我们可以使用segmented这个包。</p><p id="20PGNB88">model.segmented &lt;- segmented(model.lm)#构建分段回归模型</p><p id="20PGNB8A">summary(model.segmented)#查看模型概况</p><p class="f_center"><img src="https://nimg.ws.126.net/?url=http%3A%2F%2Fdingyue.ws.126.net%2F2023%2F0810%2F3d980a35j00rz6az00037d200bj0076g00bj0076.jpg&thumbnail=660x&quality=80&type=jpg"/><br/></p><p id="20PGNB8D"><strong>分段回归结果</strong>显示，软件自动将模型分成了两段，拐点为dosage=32.534，残差标准误为124.9，调整R2为0.7499，两个指标较曲线模型得到了进一步提升。</p><p id="20PGNB8F">#查看拟合效果</p><p id="20PGNB8G">plot(dosage,urine, pch=1, cex=1.5)</p><p id="20PGNB8H">abline(a=coef(model.lm)[1],b=coef(model.lm)[2],col="red",lwd= 2.5)</p><p id="20PGNB8I">plot(model.segmented, col='blue', lwd= 2.5 ,add=T)</p><p class="f_center"><img src="https://nimg.ws.126.net/?url=http%3A%2F%2Fdingyue.ws.126.net%2F2023%2F0810%2F465d18b7j00rz6az1002fd200u000f0g00it009e.jpg&thumbnail=660x&quality=80&type=jpg"/><br/></p><p id="20PGNB8M">在构建的上述模型中，函数自动将模型分成了两段。但根据对散点图的分析，我们认为将模型分为三段更为合适，此时可以手动设置25和30两个剂量拐点，软件会自动寻找附近的点作为**拐点。</p><p id="20PGNB8O">#手动设置拐点，分三段回归</p><p id="20PGNB8Q">model.segmented2 &lt;- segmented(model.lm,psi=c(25,30))#构建分段回归模型</p><p id="20PGNB8S">summary(model.segmented2)#查看模型概况</p><p class="f_center"><img src="https://nimg.ws.126.net/?url=http%3A%2F%2Fdingyue.ws.126.net%2F2023%2F0810%2Fa0e6f2caj00rz6az1003nd200bj007tg00bj007t.jpg&thumbnail=660x&quality=80&type=jpg"/><br/></p><p id="20PGNB8V">软件找到的两个**拐点分别为24.075和30.166，此时分段回归模型的残差标准误为99.01，调整R2为0.8427，预测效果比曲线模型明显提升。</p><p id="20PGNB91">#查看拟合效果</p><p id="20PGNB93">plot(dosage,urine, pch=1, cex=1.5)</p><p id="20PGNB94">abline(a=coef(model.lm)[1],b=coef(model.lm)[2],col="red",lwd= 2.5)</p><p id="20PGNB95">plot(model.segmented2, col='blue', lwd= 2.5 ,add=T)</p><p class="f_center"><img src="https://nimg.ws.126.net/?url=http%3A%2F%2Fdingyue.ws.126.net%2F2023%2F0810%2F155fdd38j00rz6az2002ad200u000elg00it0095.jpg&thumbnail=660x&quality=80&type=jpg"/><br/></p><p id="20PGNB99"><strong>五、样条回归</strong></p><p id="20PGNB9C">上述提到的曲线方程和分段回归两种方法都有一定的缺点。曲线方程是非局部的，当某一个因变量的值发生变化时，即使距离很远的点也会受到影响。如果采用多项式建立曲线方程，当多项式的幂较高时，自变量的一个微小变化，就会引起因变量很大的变化，得出的模型不适合外推到其他数据样本。而在分段回归模型中，每一段都是基于线性回归而建立的，拐点之间的连接显得比较生硬。</p><p id="20PGNB9E">那么有没有办法建立一个既具有分段回归模型的优点，又可以拟合比较平滑的模型呢？<strong>样条回归则兼具曲线方程和分段回归的优点，可以灵活的分段展示自变量与因变量之间的关系</strong>。样条回归把数据集划分成一个个连续的区间，划分的点称为节点，每个节点之间用单独的模型（线性函数或者低阶多项式函数）来拟合。节点越多，模型就越灵活。但是过多的节点也会导致过拟合问题，所以一般先尝试设置<strong>3个节点</strong>为宜。</p><p id="20PGNB9G">样条回归很多种，我们主要讲<strong>限制性立方样条回</strong><strong>归</strong>。<br/></p><p id="20PGNB9I">model.spline &lt;- lm(urine ~ rcs(dosage, c(20,30,35)))#建立样条回归，设置3个节点</p><p id="20PGNB9K">summary(model.spline)#查看模型概况</p><p class="f_center"><img src="https://nimg.ws.126.net/?url=http%3A%2F%2Fdingyue.ws.126.net%2F2023%2F0810%2F1d1a7b8cj00rz6az3002rd200bj005tg00bj005t.jpg&thumbnail=660x&quality=80&type=jpg"/><br/></p><p id="20PGNB9N">样条回归模型的残差标准误为139.6，调整R2为0.6872。比线性回归和曲线回归好，但不如分段回归。<br/></p><p id="20PGNB9P">#样条回归拟合效果</p><p id="20PGNB9R">ggplot(urinetest, aes(dosage, urine) ) +</p><p id="20PGNB9S">geom_point() +</p><p id="20PGNB9T">stat_smooth(method = lm, formula = y ~ rcs(x, c(20,30,35)) )</p><p class="f_center"><img src="https://nimg.ws.126.net/?url=http%3A%2F%2Fdingyue.ws.126.net%2F2023%2F0810%2Fef7b0166j00rz6az3000td200bj006wg00bj006w.jpg&thumbnail=660x&quality=80&type=jpg"/><br/></p><p id="20PGNBA2"><strong>六、Lowess函数建立局部加权回归</strong><br/></p><p id="20PGNBA5">以上介绍的模型都是参数模型，选择什么样的曲线，设置多少个拐点，这些步骤都需要进行尝试，但也会容易出现过拟合现象。于是有学者提出了Lowess非参数回归，它没有回归系数可估计，只是在寻找一条拟合效果相对更好的曲线。</p><p id="20PGNBA7">model.lowess&lt;- lowess(urine ~ dosage)#建立局部加权回归</p><p id="20PGNBA9">summary(model.lowess)#查看概况</p><p id="20PGNBAB">#查看拟合</p><p id="20PGNBAD">ggplot(urinetest, aes(dosage, urine)) +</p><p id="20PGNBAE">geom_point() +</p><p id="20PGNBAF">stat_smooth()</p><p class="f_center"><img src="https://nimg.ws.126.net/?url=http%3A%2F%2Fdingyue.ws.126.net%2F2023%2F0810%2F31e1119fp00rz6az30002d2009q002ig009q002i.png&thumbnail=660x&quality=80&type=jpg"/><br/></p><p class="f_center"><img src="https://nimg.ws.126.net/?url=http%3A%2F%2Fdingyue.ws.126.net%2F2023%2F0810%2F75f473eej00rz6az4000vd200bj006zg00bj006z.jpg&thumbnail=660x&quality=80&type=jpg"/><br/></p><p id="20PGNBAJ">局部加权给一般只做数据探索，stat_smooth()就是默认用lowess画拟合图</p><p id="20PGNBAM"><strong>七、广义可加模型</strong></p><p id="20PGNBAP">和lowess函数一样，广义可加模型也无法给出明确的系数，但它的适用范围更广，可以执行因变量与多个自变量之间的各种非参数拟合。</p><p id="20PGNBAR">它可以是任意的单变量函数的叠加，这些函数既可以是线性，也可以是非线性。它的因变量可以服从二项分布、Poisson分布、Gamma分布等更广义的范畴。它的任务就是根据目前的数据，找出一条最贴合的曲线。</p><p id="20PGNBAT">model.gam &lt;- gam(urine ~ s(dosage), data = urinetest)#建立gam模型</p><p id="20PGNBAV">summary(model.gam)#查看模型概况</p><p class="f_center"><img src="https://nimg.ws.126.net/?url=http%3A%2F%2Fdingyue.ws.126.net%2F2023%2F0810%2Fcbd17c85j00rz6az4003md200g300abg00g300ab.jpg&thumbnail=660x&quality=80&type=jpg"/><br/></p><p id="20PGNBB2">广义可加模型的调整R2为0.837，但没有给出残差标准误的结果，所以我们需要利用模型生成预测值，用预测值和真实值进行比较，得出残差标准误为98.5，是上述众多模型中表现最优秀的。</p><p id="20PGNBB4">pr.gam &lt;- predict(model.gam,urinetest)#生成预测值</p><p id="20PGNBB6">#计算RSME和R方</p><p id="20PGNBB8">data.frame(RMSE = RMSE(pr.gam, urinetest$urine),</p><p id="20PGNBBA">R2 = R2(pr.gam, urinetest$urine))</p><p class="f_center"><img src="https://nimg.ws.126.net/?url=http%3A%2F%2Fdingyue.ws.126.net%2F2023%2F0810%2F0p00rz6az50003d200cq0030g00cq0030.png&thumbnail=660x&quality=80&type=jpg"/><br/></p><p id="20PGNBBD">#查看模型拟合情况</p><p id="20PGNBBF">ggplot(urinetest, aes(dosage, urine) ) +</p><p id="20PGNBBG">geom_point() +</p><p id="20PGNBBH">stat_smooth(method = gam, formula = y ~ s(x))</p><p class="f_center"><img src="https://nimg.ws.126.net/?url=http%3A%2F%2Fdingyue.ws.126.net%2F2023%2F0810%2Faadf6031j00rz6az5001vd200nq00eag00it00bb.jpg&thumbnail=660x&quality=80&type=jpg"/><br/></p><p id="20PGNBBK">从图形可以看出，广义可加模型的曲线拟合效果非常好。虽然模型在本数据集中表现良好，但仍需要注意过拟合的情况。</p><p id="20PGNBBN"><strong>八、总结</strong></p><p id="20PGNBBQ">各个模型的拟合指标比较<br/></p><p class="f_center"><img src="https://nimg.ws.126.net/?url=http%3A%2F%2Fdingyue.ws.126.net%2F2023%2F0810%2Fe6621b3bj00rz6az60013d200d400azg00d400az.jpg&thumbnail=660x&quality=80&type=jpg"/><br/></p><p id="20PGNBBU">通过比较模型指标，虽然广义可加模型表现较好，可是它并不能提供系数，无法解释变量之间的内在联系。而结合了专业背景而建立的分段回归模型表现相对更为优异。</p><p id="20PGNBC0"><strong>分段回归的结果</strong>汇总如下：</p><p class="f_center"><img src="https://nimg.ws.126.net/?url=http%3A%2F%2Fdingyue.ws.126.net%2F2023%2F0810%2F58b6a4edj00rz6az6000ld200bv0060g00bv0060.jpg&thumbnail=660x&quality=80&type=jpg"/><br/></p><p id="20PGNBC3">为方便临床应用，我们将剂量节点取整，分别为A=24和C=30，<strong>分段回归方程的书写格式</strong>为：</p><p id="20PGNBC5">Y=β0+β1X1+β2（X1-A）X2+β3（X1-C）X3</p><p id="20PGNBC7">（1）当X1≤A时，即X≤24时，X2=0，X3=0，</p><p id="20PGNBC8">Y=β0+β1X1=2168.470-0.692*X</p><p id="20PGNBCA">（2）当A＜X1≤C，即24＜X≤30时，X2=1，X3=0，</p><p id="20PGNBCB">Y=β0+β1X1+β2（X1-A）X2=2168.470-0.692*X+92.703*(X-24)=-56.402+92.011*X</p><p id="20PGNBCD">（3）当X1＞C，即X＞30时，X2=1，X3=1，</p><p id="20PGNBCE">Y=β0+β1X1+β2（X1-A）X2+β3（X1-C）X3</p><p id="20PGNBCF">=2168.470-0.692*X+92.703*(X-24)-96.839*(X-30)=2848.768-4.828*X</p><p id="20PGNBCH">在医咖会7周年之际，我们推出了会员，包含课程、直播、专属科研工作台、科研工具、积分等12项权益。在会员正式上线之前，有一个<strong>为期15天的会员预约</strong>活动，花费1元可享50元现金抵扣。<br/></p>