精华内容
下载资源
问答
  • 限制性立方样条
    多人点赞
    2022-03-29 21:20:43

    R语言生存分析之限制性立方样条分析:为什么需要样条分析、样条回归、和多项式回归有什么不同?是否存在非线性关系、限制性立方样条cox回归模型

    目录

    更多相关内容
  • 限制性立方样条函数 (restricted cubic spline,RCS) 属于多项式中的一种,最大的特点在于它进行了样条插值并对趋势首尾两端进行了线性限制,即在第一个节点前的趋势和最后一个节点后的趋势强制限制为线性。...


    在上一篇文章等时替代模型的案例论文中,我们可以看到,COX比例风险回归模型中的连续变量可能会对生存时间产生非线性风险,如果不考虑这个风险,就会对结果造成偏差。非线性关系的拟合思路是在线性关系上加入高次项,如二次项和三次项,这样就构成了高次多项式模型。因此,作者对于非线性趋势的拟合采用了限制性立方样条(restricted cubic spline,RCS)的方法,其中也简单提到了,RCS是用来拟合自变量和因变量之间非线性关系的一个好方法。

    限制性立方样条函数属于多项式中的一种,最大的特点在于它进行了样条插值并对趋势首尾两端进行了线性限制,即在第一个节点前的趋势和最后一个节点后的趋势强制限制为线性。在这篇文章中,笔者就将对限制性立方样条的原理和案例进行总结。

    一、背景

    (一)什么是线性 (Linearity) ?

    维基百科的定义:

    Linearity is the property of a mathematical relationship (function) that can be graphically represented as a straight line. 1
    线性是指可以用图形表示为直线的数学关系或函数

    Additivity: f(x + y) = f(x) + f(y).
    Homogeneity of degree 1: f(αx) = α f(x) for all α.

    对于多维函数的推广,线性是指函数兼容加法和定标的性质,也称为叠加原理

    (二)为什么要做线性假设?

    线性是建立线性回归模型时的常见假设。

    因为现实世界里面许多量之间具有线性或者近似线性的依赖关系,即使不是线性关系,我们对变量进行适当的变换,得到的新的变量可能也会具有近似的线性关系。一旦问题落脚到线性问题,那在数学里面积累的大量的、丰富的处理线性关系的理论和方法就能用上了。

    而且线性模型具有一个最大的优势,就是易于解释。在线性关系中,连续变量在整个变量范围内具有相同的影响。

    例如,在传统的等时替代模型中,我们用系数b代表替代效应。系数b1就表示在保持总体活动时间不变的情况下,用其他活动去替换x1的时候,对健康结果变量y的估计效果。

    这很容易理解,线性使得估计易于解释,也让传统的模型更具有优势:模型中每增加一个单位就会给出相应系数在结果中的变化。

    但是,要知道,自然界中的规律很难会遵循一条简单的直线。对于今天的这些复杂的大型数据集,我们的模型应该允许,也必须允许非线性效应的存在。

    (三)如何面对非线性的难题?

    有许多非线性替代方案可以帮助我们去面对这些非线性的难题,从而更好地找到变量之间的联系。
    这些替代方案大多数都依赖于将单个连续变量转换为多个。其中最简单的方法就是把变量转换为多项式。

    例如,下面这个模型:

    我们可以对年龄这个变量进行扩展,让它还包括年龄的平方:
    在这里插入图片描述
    变成多项式后,所得的就是曲线了。这样就产生了一个问题,当我们添加平方项后,前面的系数就变得难以解释,并且如果我们再添加三次项,即Age³,那前面的系数就几乎不可能解释了。
    正是由于这种解释困难,我们要么使用 rms::contrast 函数,要么使用stats::predict来说明变量的行为方式。

    (四)多项式回归

    本篇文章的案例数据均来自下面随机生成的数据:

    options(max.print=25)
    set.seed(76)
    
    f <- function(x){
      f_x <- 0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10
      return(f_x)
    }
    
    values <- data_frame(
      x = seq(from=0, to=1, length=500),
      f_x = f(x),
      eps = rnorm(n=500, mean=0, sd=2),
      y = f_x + eps
    )
    
    values %>% 
      ggplot(aes(x=x)) +
      stat_function(fun = f) +
      geom_point(aes(y=y))
    

    对于函数f,在范围(0,1)上按照期望为0,方差为2的正态分布,随机产生500个噪点,结果如下图:

    在这里插入图片描述

    对噪点分别进行一次、二次、三次、四次多项式回归:

    ggplot(values,aes(x=x,y=y))+geom_point()
    
    lmp1=lm(y~poly(x,1),data = values)
    poly1=predict(lmp1,values)
    values$poly1=poly1
    
    lmp2=lm(y~poly(x,2),data = values)
    poly2=predict(lmp2,values)
    values$poly2=poly2
    
    lmp3=lm(y~poly(x,3),data = values)
    poly3=predict(lmp3,values)
    values$poly3=poly3
    
    lmp4=lm(y~poly(x,4),data = values)
    poly4=predict(lmp4,values)
    values$poly4=poly4
    
    
    valueslong=gather(values,key="model",value = "value",
                        c("poly1","poly2","poly3","poly4"))
    ggplot(valueslong)+theme_bw()+geom_point(aes(x,y))+
      geom_line(aes(x=x,y=value,linetype=model,colour=model),size=1.0)+
      theme_update(axis.title.y=element_blank())
    

    在这里插入图片描述

    多项式的使用是一种粗略的估算方法,具有可以产生平滑拟合的优点,但是它所拟合的函数也十分复杂,最大的缺点就是不灵活。

    例如,二次曲线的形状必须是“U”形(或倒“U”形)——上升的必须以相同的速度下降;多项式通常会导致曲线在末端不能很好地拟合数据。

    因此,对于多项式,我们也需要找到更好的方法去替代它,本篇文章所要介绍的方法——样条曲线 (Spline) ,就是一类灵活的替代函数。

    二、样条函数

    (一)样条 (Spline) 的定义

    1. 原始定义

    样条(spline)原本是指工匠的工具,它是一种灵活的细木条或金属条,用来绘制平滑曲线。一般会用几个重物施加在工具的不同位置,让条带根据重物的数量和位置发生弯曲。
    在这里插入图片描述

    2. 数学定义

    类似的作用,样条曲线用于统计,是为了方便在数学上重现灵活的形状。样条曲线本质是一个分段多项式函数,此函数受限于某些控制点,称为 “节点”,节点放置在数据范围内的多个位置,多项式的类型以及节点的数量和位置决定了样条曲线的类型。

    样条函数是一组平滑连接的分段多项式。 “平滑连接”意味着对于 n 次多项式,样条函数及其前 n-1 阶导数在节点处都是连续的。

    (二)节点的选择

    1. 节点 (knots)

    如果使用 k 个结,拟合 n 次多项式需要 k + n + 1 个回归参数(包括截距)。

    在实践中,通常使用三次样条(即 3 次多项式),除截距外还需要 k + 3 个系数,而拟合线性模型只需 1 个系数。 这是允许拐点的多项式的最小度数,为拟合数据提供了足够的灵活性,同时不像高阶样条那样需要那么多的自由度。 由于一阶和二阶导数(斜率和斜率的变化率)在节点处是连续的,三次样条曲线在外观上是平滑的2
    在这里插入图片描述

    2. R软件拟合示例

    R示例:

    library(tidyverse)
    library(broom)
    
    df_values <- c(2, 5, 10, 15, 25, 50)
    overall <- NULL
    
    for(df in df_values){
      overall <- smooth.spline(values$x, values$y, df=df) %>%
        augment() %>%
        mutate(df=df) %>%
        bind_rows(overall)
    }
    multiple_df <- overall %>% 
      ggplot(aes(x=x)) +
      geom_point(aes(y=y)) +
      geom_line(aes(y=.fitted)) +
      facet_wrap(~df, nrow=2) +
      labs(title="Splines fit with different degrees of freedom")
    
    multiple_df + 
      stat_function(fun = f)
    

    如下图所示,是节点分别为2、5、10、15、25、50时,利用R软件,对黑色噪点进行拟合得到的图形。(蓝色是拟合曲线,红色是真实曲线)

    在这里插入图片描述

    (三)样条的分类

    1. B-splines(B样条/多次回归样条)

    如果样条是基曲线的线性组合, 则称为B样条(B-spline)
    基函数为:
    在这里插入图片描述
    B-spline定义为:
    在这里插入图片描述

    在这里插入图片描述

    2. Natural splines(自然样条/限制性立方样条)

    在这里插入图片描述

    ‘Natural Cubic Spline’ — is a piece-wise cubic polynomial that is twice continuously differentiable. It is considerably ‘stiffer’ than a polynomial in the sense that it has less tendency to oscillate between data points.
    在这里插入图片描述

    三、限制性立方样条

    (一)定义

    回归样条(regression spline)本质上是一个分段多项式,但它一般要求每个分段点上连续并且二阶可导。

    即:设自变量数据的范围在区间[a,b],并根据需要分成k个段:a = t0 < t1 < … < tk-1 < tk = b,在每个区间 [ti-1,ti) 分别用一个多项 Si(x) 式表示,则回归样条 f(x) = Si(x) 当 x∈[ ti-1,ti ),并且 f ″(x) 在 [a,b] 存在连续3

    限制性立方样条 (restricted cubic spline,RCS) 是在回归样条的基础上附加要求:样条函数在自变量数据范围两端的两个区间 [t0,t1) 和 (tk-1,tk] 内是线性函数,满足上述性质的样条函数通常用 RCS(X) 表示,其表达式为:

    在这里插入图片描述
    更简单地一种理解:
    在这里插入图片描述

    (二)RCS节点的数量和位置

    结的数量比位置更重要。有研究表明,五个结足以很好地适应实践中可能出现的任何模式。

    Harrell (2010) 指出“对于许多数据集,k = 4 提供了模型的充分拟合,并且是在灵活性和由过度拟合小样本引起的位置损失之间的良好折衷”4

    • 如果样本量很小,则应使用三个节点,以便在节点之间有足够的观测值,以便能够拟合每个多项式。
    • 如果样本量很大并且有理由相信所研究的关系变化很快,则可以使用超过 5 个节点。

    对于一些研究,理论可能会暗示结的位置。更常见的是,节点的位置是根据连续变量的分位数预先指定的,这确保了每个区间中有足够的观测值来估计三次多项式。
    在这里插入图片描述

    (三)R软件rms程序包操作示例

    x = seq(from=0, to=1, length=500)
    y = f(x)
    model <- lm(y ~ x)
    
    # restricted cubic spline
    model.spline1 <- lm(y ~ rcs(x, 3))
    model.spline2 <- lm(y ~ rcs(x, 5))
    model.spline3 <- lm(y ~ rcs(x, 7))
    model.spline4 <- lm(y ~ rcs(x, 9))
    model.spline5 <- lm(y ~ rcs(x, 10))
    
    # compare models (to determine number of knots)
    AIC(model.spline1,model.spline2,model.spline3,model.spline4,model.spline5)
    
    

    比较不同节点数对应的AIC值:
    在这里插入图片描述

    ggplot(data=values, aes(x=x, y=y))  +
      geom_point()+geom_smooth(aes(x=x, y=y), 
                               method="lm", formula=y ~ rcs(x, 3), 
                               se=FALSE, colour="blue") +
      geom_point()+geom_smooth(aes(x=x, y=y), 
                               method="lm", formula=y ~ rcs(x, 5), 
                               se=FALSE, colour="red") +
      geom_point()+geom_smooth(aes(x=x, y=y), 
                               method="lm", formula=y ~ rcs(x, 9), 
                               se=FALSE, colour="orange") 
    

    在这里插入图片描述

    四、相关文献阅读

    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    参考文献


    1. 维基百科:https://en.wikipedia.org/wiki/Linearity ↩︎

    2. Croxford, Ruth. “Restricted Cubic Spline Regression : A Brief Introduction.” (2016). ↩︎

    3. 罗剑锋,金欢,李宝月,赵耐青.限制性立方样条在非线性回归中的应用研究[J].中国卫生统计,2010,27(03):229-232. ↩︎

    4. Harrell, F.E. (2010) Regression Modeling Strategies: with applications to linear models, logistic regression, and survival analysis. Springer-Verlag New York, Inc. New York, USA. ↩︎

    展开全文
  • 我用stata跑出来的限制性立方样条图类似"✓"型 但是用testparm指令跑出来的非线性检验Pnon-linearity值却>0.05,即认为是线性 两个对应不起来是什么原因呢? 感谢大家帮助
  • 相关介绍: 在病因推断、剂量效应研究中,时常要分析自变量和因变量的数量关系。广义线性模型,如Logistic回归、Possion回归等是应用比较广泛的方法。它的一个重要假设是通过...限制性立方样条中节点的个数(k)和...

    图片

    • 相关介绍:

    在病因推断、剂量效应研究中,时常要分析自变量和因变量的数量关系。广义线性模型,如Logistic回归、Possion回归等是应用比较广泛的方法。它的一个重要假设是通过选择合适的链接函数,因变量与自变量的关系呈线性。这个假设在某些情况下并不成立。

    此时一个常见的处理是采用百分位数等方法将连续性变量分段(P value for trend)。但是分段往往主观,而且损失信息,并有可能引入偏倚。本文介绍限制性立方样条拟合自变量和因变量之间的非线性关系。

    限制性立方样条中节点的个数(k)和位置可以由研究人员根据研究背景和数据选择。当有较强的背景知识支持,知道自变量和应变量的关系在某些特定节点转折时,可以选择这些特定的节点。但是实际上往往没有足够的背景知识足以确定节点的个数和位置。所幸绝大多数情况下,节点的位置对限制性立方样条的拟合影响不大,相对来说节点的个数是更关键的参数。节点的个数决定曲线的形状,或者说平滑程度。当节点的个数为2时,得到的拟合曲线就是一条直线。当节点个数等于样本量时,相当于将各个点用线段相连,得到的是完全拟合但是不平滑的折线。由于节点个数的选择和自由度有关,所以当样本量比较大的时候可以取较多的节点。但是,一般来说节点个数取3~7就足够了。



    • RCS案例图形:

    图片


    • R语言代码实现:

      利用library()函数加载本案例所需R包,R程序代码如下:

    library(smoothHR)
    library(survival)
    library(rms)
    library(Hmisc)

    导入数据,并查看变量名称,R程序代码如下:

    RCS <- read.csv("RCS.csv")
    names(RCS

    因变量group 、自变量VO2max、协变量Age、BMI。

    寻找最优节点数knots,R程序代码如下:

    for (i in 3:7) {
      fit <- glm(group~rcs(BMI,nk=i)+Age+VO2max,data=RCS, family = binomial())
      tmp <- extractAIC(fit)
      if(i == 3) {AIC = tmp[2]; nk = 3}
      if(tmp[2] < AIC) {AIC = tmp[2]; nk = i} 
    }

    以上代码为根据AIC准则筛选最小AIC对应的knots。结果显示最优knots为3。

    开始拟合模型,R程序代码如下:

    # 连续变量BMI被样条化
    fit <- lrm(group~rcs(BMI,nk=nk)+Age+VO2max,data=RCS, x=TRUE,y=TRUE)
    anova(fit)

    输出结果如下:

     Wald Statistics     Response: group  Factor     Chi-Square d.f.   P      BMI         97.12     4    <.0001  Nonlinear  6.35     3    0.033 Age         27.62     1    <.0001 VO2max       0.21     1    0.6444 TOTAL      111.02     6    <.0001

    结果显示非线性的P值=0.033,提示存在非线性。

    将非线性P值提取出来,方便后续作图,R程序代码如下:​​​​​​​

    ##p-non-linear
    p <-round(anova(fit)[,3],3)

    设置参考值,这是一个具有生物/医学意义的参考值。无论参考值为何,该点处的OR为1,R程序代码如下:

    refvalue <- 21.75 

    对数据进行打包,并指定参考值,R程序代码如下:​​​​​​​

    ddist<-datadist(RCS)
    ddist$limits$BMI[2]<-refvalue
    options(datadist="ddist")

    将模型预测的OR值储存在pred_OR中,R程序代码如下:

    pred_OR<-Predict(fit,BMI,ref.zero=TRUE,fun=exp)

    开始绘制图形:​​​​​​​

    # 设置密度曲线的背景色
    violet <- "#89439B"
    # 绘制左右双坐标轴baseplot
    par(mar = c(5, 4, 4, 4) + 0.3)
    par(xpd=NA)
    ylim.bot<-min(pred_OR[,"lower"])
    ylim.top<-max(pred_OR[,"upper"])
    # 先画密度图以免遮挡下面的线图
    dens <- density(RCS$BMI) # 计算密度
    plot(dens$x,dens$y, col=ggplot2::alpha(violet,0.5), type="l", xlab = "", ylab = "",xaxt="n",yaxt="n")
    polygon(dens$x,dens$y,col = ggplot2::alpha(violet,0.5),border = ggplot2::alpha(violet,0.5)) # 颜色透明化防遮盖线条
    par(new=TRUE) # 新增画布
    plot(pred_OR[,1],pred_OR[,"yhat"], 
         xlab = "BMI",ylab = paste0("OR where the refvalue for BMI is ",refvalue),
         type = "l",ylim = c(ylim.bot,ylim.top),
         col="red",lwd=2) 
    lines(pred_OR[,1],pred_OR[,"lower"],lty=2,lwd=1.5)
    lines(pred_OR[,1],pred_OR[,"upper"],lty=2,lwd=1.5)
    lines(x=range(pred_OR[,1]),y=c(1,1),lty=3,col="grey40",lwd=1.3) #
    points(refvalue,1,pch=16,cex=1.2)
    #附上refvalue的值,具体位置可以自己修改
    text(refvalue + 2, 1.1, paste0("refvalue = ",refvalue)) 
    # 绘制图例,注意非线性p值在变量p中的位置
    legend("topright",
           paste0("P-overall ",ifelse(round(p[1],3) < 0.001,"< 0.001",round(p[1],3)),
                  "\nP-non-linear = ",ifelse(round(p[2],3) < 0.001,"< 0.001",round(p[2],3))),
           bty="n",cex=0.8)
    legend("bottomleft",lty = c(1,2),col = c("red","black"),
           c("Estimation","95% CI"),
           bty="n",cex=0.8)

    输出图形如案例图形:

    图片

    展开全文
  • R语言如何确认限制性样条分析的最佳节点个数、进行方差分析通过p值确认指定连续变量和风险值HR之间是否存在非线性关系、限制性立方样条cox回归模型

    R语言如何确认限制性样条分析的最佳节点个数、进行方差分析通过p值确认指定连续变量和风险值HR之间是否存在非线性关系限制性立方样条cox回归模型

    目录

    展开全文
  • R语言使用cph函数和rcs函数构建限制性立方样条cox回归模型
  • 一、背景介绍 最近在研究怎么处理论文数据,各种分析软件都有使用,比如:SPSS、Origin、stata16、medcalc和R语言都有些研究,其中除R语言外都是收费的。不过经过一番功夫,我这边有SPSS、stata16、Origin和medcalc...
  • R语言生存分析之限制性立方样条(RCS, Restricted cubic spline)分析详解实战:拟合连续性自变量和事件风险之间的关系:基于survival包lung数据 目录 R语言生存分析之限制性立方样条(RCS, Restricted cubic spline...
  • R语言rms包生存分析之限制性立方样条(RCS, Restricted cubic spline)分析:拟合连续性自变量和事件风险之间的关系并绘制直方图、平滑曲线、双Y轴于同一个图像中 目录 R语言rms包生存分析之限制性立方样条(RCS, ...
  • R语言使用Predict函数计算指定连续变量和风险比HR值的关系、基于限制性立方样条分析方法、限制性立方样条cox回归模型
  • R语言限制性立方样条(RCS, Restricted cubic spline)分析:基于logistic回归模型、南非心脏病数据集(South African Heart Disease) 目录 R语言限制性立方样条(RCS, Restricted cubic spline)分析:基于...
  • 因此,一个更好的解决方法是拟合自变量与因变量之间的非线性关系,限制性立方(Restricted cubic spline,RCS)就是分析非线性关系的最常见的方法之一。 既往文章《stata两种方法制作限制立方条图》中,我们已经介绍了...
  • R语言使用cph函数和rcs函数构建限制性立方样条cox回归模型、使用ggcoxzph函数可视化进行Schoenfeld残差图检验模型是否满足等比例风险
  • 请问R中怎么看限制立方样条曲线与HR=1的截点?请问应该输什么代码看交点数据呢?
  • R语言使用cph函数和rcs函数构建限制性立方样条cox回归模型、使用rms包的Predict函数计算指定连续变量和风险比HR值的关系、可视化连续变量和风险值HR的关系
  • R语言使用cph函数和rcs函数构建限制性立方样条cox回归模型、使用rms包的Predict函数计算指定连续变量在不同分组变量下和风险比HR值的关系、使用ggplot2可视化连续变量在不同分组变量下和风险值HR的关系、美化可视化...
  • R语言使用cph函数和rcs函数构建限制性立方样条cox回归模型、使用cox.zph函数执行PH检验、检验模型是否满足等比例风险
  • R语言使用cph函数和rcs函数构建限制性立方样条cox回归模型、检验模型是否满足等比例风险、是否存在非线性关系、使用rms包的Predict函数计算指定连续变量和风险比HR值的关系并可视化
  • 实例演示R语言制作限制性立方条图

    千次阅读 热门讨论 2021-03-31 11:58:46
    限制性立方样条函数(RCS)在比较非线性关系中很常用。既往我们已经讲过R语言制作限制性立方条图,但是讲得比较简单,中间有些环节没写出来,我也不是很满意,今天重新来说一下。主要是要用到rms包的rcs函数来绘制,...
  • R语言样条回归并绘制限制立方条图

    千次阅读 热门讨论 2021-02-18 10:18:33
    因此,一个更好的解决方法是拟合自变量与因变量之间的非线性关系,限制性立方样条(Restricted cubic spline,RCS)就是分析非线性关系的最常见的方法之一。 最近有粉丝问有序多分类变量可以应用在限制立方条图吗?是...
  • 想问一下立方样条曲线的置信区间的问题,我应该怎么调整置信区间的图像成喇叭形,代码奉上,第一张是自己画的图,第二张是文献的图。 library(ggplot2) library(survival) library(rms) setwd("D:/work/work/zhang 9...
  • 当自变量和因变量之间为非线性关系时,可以将连续型变量转化为分类变量,但是分类变量的类别数目以及节点位置的选择一般会带有主观并且分类变量会损失部分信息;也可以直接拟合自变量和因变量之间的非线性关系,...
  • 大佬们,我用r做限制性立方样条看anova非线性p值的时候结果后面很有error,这个是什么原因呀

空空如也

空空如也

1 2 3 4 5 ... 18
收藏数 355
精华内容 142
关键字:

限制性立方样条

友情链接: displ.zip