草庐IT

R实战 | 限制性立方样条(RCS)

木舟笔记 2023-05-03 原文

RCS

在科学研究中，我们经常构建回归模型来分析自变量和因变量之间的关系。大多数的回归模型有一个重要的假设就是自变量和因变量呈线性关联。当自变量和因变量之间为非线性关系时，可以将连续型变量转化为分类变量，但是分类变量的类别数目以及节点位置的选择一般会带有主观性并且分类变量会损失部分信息；也可以直接拟合自变量和因变量之间的非线性关系，但是直接构建多项式回归可能存在过度拟合、共线性等问题。因此，一个更好的解决方法是拟合自变量与因变量之间的非线性关系，「限制性立方样条」(Restricted cubic spline，RCS)就是分析非线性关系的最常见的方法之一。

样条（spline）原本是指是一种灵活的细木条或金属条，用来绘制平滑曲线。样条曲线本质是一个分段多项式函数，此函数受限于某些控制点，称为 “节点”，节点放置在数据范围内的多个位置，多项式的类型以及节点的数量和位置决定了样条曲线的类型。立方则指的是 函数为3次多项式。限制是在回归样条的基础上附加要求：样条函数在自变量数据范围两端的两个区间 [X1,X2) 和 (Xn-1,Xn] 内是线性函数。

RCS曲线

RCS节点的数量比位置更重要。由于节点个数的选择和自由度有关, 所以当样本量比较大的时候可以取较多的节点。但是节点越多, 自由度越大, 模型越复杂, 越难解在「«Regression Modeling Strategies»」这本书中，Harrell 建议节点数为4时，模型的拟合效果较好，即同时可以兼顾曲线的平滑程度以及避免过拟合造成的精确度降低。当样本量较大时，5个节点是更好的选择。小样本（n<30）可以选择3个节点。当节点的个数为2时，得到的拟合曲线就是一条直线。大多数研究者推荐的节点为3-5个。

一个例子

❝

Lee DH, Keum N, Hu FB, et al. Predicted lean body mass, fat mass, and all cause and cause specific mortality in men: prospective US cohort study. BMJ. 2018;362:k2575. Published 2018 Jul 3. doi:10.1136/bmj.k2575

❞
Association of predicted body composition and  body mass index (BMI)* with all cause mortality in men
❝

如图，为了探究预测的FM、LBM和BMI与男性全因死亡率的关系，作者分别对这三个因素做了RCS分析。

❞

RCS实战

#加载所需要的包
library(ggplot2)
#install.packages('rms')
library(rms)   
# 导入示例数据
data <- read.csv('test.csv')
head(data)
> head(data)
       age    sex      time death
1 60.57519   Male  3.094579     1
2 42.11447   Male  1.574237     0
3 54.86611   Male  3.239313     0
4 55.82207   Male 12.495977     0
5 52.48256 Female  3.252534     0
6 46.12436   Male  2.836695     0
# 对数据进行打包，整理
dd <- datadist(data) #为后续程序设定数据环境
options(datadist='dd') #为后续程序设定数据环境
# 拟合模型
fit<- cph(Surv(time,death) ~ rcs(age,4) + sex,data=data)  # 节点数设为4
# 非线性检验
# P<0.05为存在非线性关系
anova(fit)
> anova(fit)
                Wald Statistics          Response: Surv(time, death) 

 Factor     Chi-Square d.f. P     
 age        57.75      3    <.0001
  Nonlinear  8.17      2    0.0168
 sex        18.75      1    <.0001
 TOTAL      75.63      4    <.0001
# 查看HR预测表
# 看一下预测的HR所对因的age
HR<-Predict(fit, age,fun=exp,ref.zero = TRUE)
head(HR)
> head(HR)
       age  sex      yhat     lower    upper
1 19.71985 Male 0.7087866 0.2403429 2.090257
2 20.00869 Male 0.7052492 0.2429359 2.047356
3 20.29754 Male 0.7017294 0.2455536 2.005363
4 20.58638 Male 0.6982271 0.2481960 1.964258
5 20.87523 Male 0.6947423 0.2508632 1.924024
6 21.16408 Male 0.6912750 0.2535552 1.884643

Response variable (y):  

Adjust to: sex=Male  

Limits are 0.95 confidence limits
# 绘图
ggplot()+
  geom_line(data=HR, aes(age,yhat),
            linetype="solid",size=1,alpha = 0.7,colour="#0070b9")+
  geom_ribbon(data=HR, 
              aes(age,ymin = lower, ymax = upper),
              alpha = 0.1,fill="#0070b9")+
  theme_classic()+
  geom_hline(yintercept=1, linetype=2,size=1)+
  geom_vline(xintercept=48.89330,size=1,color = '#d40e8c')+ #查表HR=1对应的age
  labs(title = "Stroke Risk", x="Age", y="HR (95%CI)")
RCS1
❝

可以看到，年龄<49岁，死亡风险随年龄变化不是很明显；年龄>49岁后，死亡风险随年龄的增加而增加。

❞

分组

### 性别分组
HR1 <- Predict(fit, age, sex=c('Male','Female'),
               fun=exp,type="predictions",
               ref.zero=TRUE,conf.int = 0.95,digits=2)
HR1
ggplot()+
  geom_line(data=HR1, aes(age,yhat, color = sex),
            linetype="solid",size=1,alpha = 0.7)+
  geom_ribbon(data=HR1, 
              aes(age,ymin = lower, ymax = upper,fill = sex),
              alpha = 0.1)+
  scale_color_manual(values = c('#0070b9','#d40e8c'))+
  scale_fill_manual(values = c('#0070b9','#d40e8c'))+
  theme_classic()+
  geom_hline(yintercept=1, linetype=2,size=1)+
  labs(title = "Stroke Risk", x="Age", y="HR (95%CI)")
RCS2

示例数据和代码领取

点赞、在看 本文，分享至朋友圈集赞25个并保留30分钟，截图发至微信mzbj0002领取。

「木舟笔记2022年度VIP可免费领取」。

木舟笔记2022年度VIP企划

「权益：」

  1. 「2022」年度木舟笔记所有推文示例数据及代码(「在VIP群里实时更新」)。

    资源合集

  2. 木舟笔记「科研交流群」。

  3. 「半价」购买跟着Cell学作图系列合集(免费教程+代码领取)|跟着Cell学作图系列合集。

「收费：」

「99￥/人」。可添加微信：mzbj0002 转账，或直接在文末打赏。

参考

  1. https://blog.csdn.net/weixin_43645790/article/details/125285467

  2. 限制性立方样条图，一种美的不行的趋势性分析方法（附R语言详细教程

  3. Restricted cubic splines. A spline is a drafting tool for drawing… | by Peter Flom | Towards Data Science

  4. R语言绘制限制性立方样条（Restricted cubic spline，RCS） - 简书 (jianshu.com)

往期内容

  1. CNS图表复现|生信分析|R绘图 资源分享&讨论群！

  2. 组学生信| Front Immunol |基于血清蛋白质组早期诊断标志筛选的简单套路

  3. (免费教程+代码领取)|跟着Cell学作图系列合集

  4. Q&A | 如何在论文中画出漂亮的插图？

  5. 跟着 Cell 学作图 | 桑葚图(ggalluvial)

  6. R实战 | Lasso回归模型建立及变量筛选

  7. 跟着 NC 学作图 | 互作网络图进阶(蛋白+富集通路)(Cytoscape)

  8. R实战 | 给聚类加个圈圈(ggunchull)

  9. R实战 | NGS数据时间序列分析(maSigPro)

  10. 跟着 Cell 学作图 | 韦恩图(ggVennDiagram)

木舟笔记矩阵
RCS实战codexffxff0c机器学习python人工智能大数据数据挖掘

有关R实战 | 限制性立方样条(RCS)的更多相关文章

  1. 微信小程序开发入门与实战（Behaviors使用） - 2

    @作者:SYFStrive @博客首页:HomePage📜：微信小程序📌：个人社区（欢迎大佬们加入）👉：社区链接🔗📌：觉得文章不错可以点点关注👉：专栏连接🔗💃：感谢支持,学累了可以先看小段由小胖给大家带来的街舞👉微信小程序（🔥）目录自定义组件-behaviors    1、什么是behaviors    2、behaviors的工作方式    3、创建behavior    4、导入并使用behavior    5、behavior中所有可用的节点    6、同名字段的覆盖和组合规则总结最后自定义组件-behaviors    1、什么是behaviorsbehaviors是小程序中，用于实现

  2. ruby-on-rails - 限制 will_paginate 中的页数 - 2

    因此，在使用Sphinx时，搜索限制为1000个结果。但是，如果will_paginate生成的结果分页链接超过1000个，请不要考虑这一点，并提供指向超过1000/per_page的页面的链接。设置最大页数或类似内容的明显方法是什么？干杯。 最佳答案 我认为最好将参数:total_entries提交给方法paginate:@posts=Post.paginate(:page=>params[:page],:per_page=>30,:total_entries=>1000)will_paginate将仅为显示1000个结果所需的页

  3. ruby - 评论限制 - 2

    在ruby​​1.9中，放宽了行结束位置的条件，因此我们现在可以用句号开始一行来显示方法调用。当我们混淆了链式和非链式方法，并希望显示下一个非链式方法的开始位置时，这很方便。如果没有这个新功能，我们能做的最好的可能就是使用缩进:method1(args1).method2(args2).method3(args3)method4(args4).method5(args5).method6(args6)或插入一个空行。但这很不方便，因为我们必须注意缩进，同时，不要忘记在每个方法调用之后加上链中最后一个方法调用之后的句点。正因为如此，我制造了很多错误，要么有一个额外的周期，要么有一个缺失的

  4. ruby - 如何在特定队列中推送作业并使用 sidekiq 限制工作人员数量？ - 2

    我知道我们可以做到:sidekiq_optionsqueue:"Foo"但在这种情况下，Worker只分配给一个队列:“Foo”。我需要在特定队列中分配作业(而不是worker)。使用Resque很容易:Resque.enqueue_to(queue_name,my_job)另外，为了并发问题，我需要限制每个队列的Worker数量为1。我该怎么做？ 最佳答案 您可能会使用https://github.com/brainopia/sidekiq-limit_fetch然后:Sidekiq::Client.push({'class'=>

  5. ruby-on-rails - 如何限制模型每天创建一条记录？ - 2

    业务逻辑:用户每天只能为日记创建一个条目。在创建条目之前，它必须查询记录以确定是否已经为今天创建了条目。我正在寻找解决此问题的最佳方法的建议。我对如何在客户端实现它有一些想法，但我真的很想在模型层进行验证。任何帮助将不胜感激。 最佳答案 在日志表上创建唯一索引:add_index:journal_entries,[:user_id,:created_on],unique:true然后只能创建一条具有给定user_id和日期的记录，如果违反，数据库将引发异常。请注意，created_on必须是date列，而不是datetime。这是唯

  6. ruby-on-rails - 使用 PostgreSQL 适配器限制 ActiveRecord 迁移 5.0 中的文本列 - 2

    我的迁移看起来像这样classCreateQuestionings现在，当我运行$rakedb:migrate:reset时，在我的db/schema.rb中看不到限制:create_table"questionings",force::cascadedo|t|t.text"body",null:falseend我做错了吗还是这是一个错误？顺便说一下，我使用的是rails5.0.0.beta3和ruby​​2.3.0p0。 最佳答案 t.text在PostgreSQL和textdoesn'tallowforsizelimits中生成

  7. 你真正了解什么是接口测试么？接口实战一“篇”入魂 - 2

    最近在工作中，看到一些新手测试同学，对接口测试存在很多疑问，甚至包括一些从事软件测试3，5年的同学，在聊到接口时，也是一知半解；今天借着这个机会，对接口测试做个实战教学，顺便总结一下经验，分享给大家。计划拆分成4个模块跟大家做一个分享，（接口测试、接口基础知识、接口自动化、接口进阶）感兴趣的小伙伴记得关注，希望对你的日常工作和求职面试，带来一些帮助。注：文章较长有5000多字，希望小伙伴们认真看完，当然有些内容对小白同学不是太友好，如果你需要详细了解其中的一些概念或者名词，请在文章之后留言，后续我将针对大家的疑问，整理输出一些大家感兴趣的文章。随着开发模式的迭代更新，前后端分离已不是新的概念，

  8. ruby - 使用 gsub 时如何限制替换次数？ - 2

    在Ruby中如何限制String#gsub的替换次数？在PHP中，这可以通过preg_replace轻松完成，它接受一个参数来限制替换，但我不知道如何在Ruby中做到这一点。 最佳答案 您可以在gsub循环中创建一个计数器并递减它。str='aaaaaaaaaa'count=5pstr.gsub(/a/){ifcount.zero?then$&elsecount-=1;'x'end}#=>"xxxxxaaaaa" 关于ruby-使用gsub时如何限制替换次数？，我们在StackOverf

  9. ruby - Ruby 进程如何限制其 CPU 使用率？ - 2

    假设我希望Ruby进程使用的CPU不超过15%。是否可以？怎么办？ 最佳答案 您可以尝试使用Process.setrlimit来自标准核心:Setstheresourcelimitoftheprocess.这看起来只是setrlimit的包装器来自C库，因此它可能仅在Unix-ish平台上可用。setrlimit不支持CPU百分比限制，但它支持以秒为单位限制CPU时间。如果您只是想让您的Ruby进程不占用整个CPU，那么您可以尝试使用Process.setpriority来调整它的优先级。这只是libc的setpriority的包装

  10. ruby - 仅使用必需的提供程序构建 Fog gem 并限制依赖项 - 2

    我正在使用出色的Foggem来访问Rackspace云文件服务。我面临的挑战是，我正在努力使访问CloudFiles的服务保持轻量级，而且Fog似乎通过其灵active具有很多我永远不需要的依赖项和代码。有没有人尝试过构建Fog的精简副本，只包含一部分提供者，从而限制依赖性？例如，专门针对Rackspace云文件API，我希望能够在没有net-ssh、net-scp、nokogirigems以及亚马逊、Rackspace和其他20个未使用的提供商的所有未使用代码的情况下处理所有内容用过的。我希望避免在每次这些未使用的提供程序之一发现错误时升级gem，同时减少我的内存占用。如果任何人在这

随机推荐