医学知识

订阅

MedSci主页 / 所有栏目 / 医学知识

如何使用科研神器「R语言」，让真实世界临床研究离SCI更进一步？一篇就够！（附代码）

2022-09-23 烤鸭小笼包 MedSci原创

R语言真实世界研究
关注

这篇文章我们重点介绍在R中的实现。

前两篇文章中，小编向大家介绍了真实世界研究中的混杂因素及其常用的校正方法，特别是倾向性评分匹配和逆概率加权，以及SAS中的实现。这篇文章我们重点介绍在R中的实现。

 

还是先回顾一下数据背景。我们用的数据来源是荷兰格罗宁根市肾脏和血管终末期疾病预防（PREVEND）的队列研究（Joosten, 2014）。

 

原数据含有较多变量，小编对数据进行了简化。其中，RFFT是结局指标，是一个连续型变量；我们想要研究的自变量是是否使用他汀类药物（Statin）。此外，还包括了年龄（Age），性别（Gender），种族（Ethnicity），教育（Education），心血管疾病（CVD），糖尿病（DM），吸烟情况（Smoking），高血压（Hypertension）这些协变量。其中，除了Age是连续型变量外，其他均为分类变量。
 
数据共有4,095名患者，其中904名在治疗组（statin =1），3,190名在对照组 （statin=0）。以下是数据的部分截图：
图片

图1. 原数据部分截图

倾向性评分匹配（PSM）
倾向性评分匹配在R中的实现是在“MatchIt” package里实现的，所以记得使用下列code之前先安装MatchIt包。我们用同样的数据statin进行示范。
 
psm_out <- matchit(Statin ~ Age + Gender + Ethnicity + Education + CVD + DM + Smoking + Hypertension, data = as.data.frame(statin), method = "nearest", distance = "logit", replace = FALSE, caliper = 0.25)
psm_out
plot(psm_out,type = "histogram", interactive = F)
plot(psm_out,type = "QQ", interactive = F, which.xs = "Age")
summary(psm_out, standardize = T)
有了之前SAS操作的基础相信对于R代码的理解会容易一些，因为很多东西大同小异，只是个别语句的syntax略有不同。我们还是逐条来看：
 
/* psm_out <- matchit(Statin ~ Age + Gender + Ethnicity + Education + CVD + DM + Smoking + Hypertension, data = as.data.frame(statin), method = "nearest", exact = “Gender”, distance = "logit", replace = FALSE, caliper = 0.25) */
 
👆这一步把计算PS和匹配同时完成了，用的是matchit() 函数，前半部分Y~X就是建模。
 
data = 明确了我们进行匹配的数据集是statin。
 
Method = 明确了匹配方式是最近邻匹配，在R中用nearest表示，此外还有cem, exact, full, optimal等选项。
 
Exact = 依然表示精确匹配，这里我们依然选择Gender为精确匹配项。
 
Distance =也和之前SAS中的相同，不过SAS中用lps表示，这里用logit。
 
Replace = 表示选择是否允许重复使用对照组的患者与研究组进行匹配。
 
Caliper同样设定为0.25。
/* psm_out */
 
图片
 
👆这一步就是生成一个总结，告诉我们一共4095个患者我们匹配了其中1418个共709对。在上述SAS的操作中我们匹配了720对，还是比较相近的。
 
/* plot(psm_out,type = "histogram", interactive = F) */
👆这一步是对PS生成直方图，观察匹配前后PS的分布。可以看到PS在匹配后的分布非常相似。
 
图片
/* plot(psm_out,type = "QQ", interactive = F, which.xs = "Age") */
 
👆这一步是用QQ图来看变量在匹配前后的分布。因为QQ图只能看连续型变量，而我们这个例子中协变量里的连续型变量只有Age，所以我们看一下Age的分布。这里横坐标表示control group，纵坐标表示treated group。
 

图片

 

/* summary(psm_out, standardize = T) */
 
👆Summary()函数可以生成一个量化的表格，其中我们比较关注的一列依然是SMD，可以看出PS和所有协变量的SMD都在±0.1以内，说明匹配效果较好。
 
图片
逆概率加权（IPW）
 
R中有多种办法可以实现逆概率加权，因为IPW本质上讲就是用logistic回归计算出PS，再根据1/PStreatment和1/(1- PStreatment)两个简单公式对观察值赋予权重。
 
大部分网上教程，是用最基础的包手动进行计算和整理。所以，这里想为大家介绍一个专门为逆概率加权设计的包，”ipw”。
 
IPW有几个主要的function：
 
  • ipwpoint: 用于”point treatment situation”，即混杂因素和我们想探究的结果（outcome）是静态的；
  • ipwtm: 用于随时间变化的自变量和混杂因素；
  • tstartfun: 用于计算基线到随访各个阶段的时间，通常和Cox风险回归结合起来用在生存性分析研究中。
 
这里我们主要讲ipwpoint。
 
图片
wt <- ipwpoint(exposure = Statin, family = "binomial", link = "logit", numerator = ~1, denominator = ~ Age + Gender + Ethnicity + Education + CVD + DM + Smoking + Hypertension, data = statin)
summary(wt$ipw.weights)
summary(wt$num.mod)
summary(wt$den.mod)
statin$sw <- wt$ipw.weights
ipw_out <- (svyglm(RFFT ~ Statin, design = svydesign(~ Age + Gender + Ethnicity + Education +
          CVD + DM + Smoking + Hypertension, weights = ~sw, data = statin)))
coef(ipw_out)
confint(ipw_out)
下面是逐条解释：
 
/* wt <- ipwpoint(exposure = Statin, family = "binomial", link = "logit", numerator = ~1, denominator = ~ Age + Gender + Ethnicity + Education + CVD + DM + Smoking + Hypertension, data = statin) */
 
👆这一步生成的是权重，我们把这一步骤产生的结果命名为wt。在ipwpoint()函数中，
 
exposure = 是我们想探究的exposure变量，即是否用statin，
 
family =和link = 定义了我们所用的函数，因为我们仍然是用logistic回归来计算PS所以两个选项分别是binomial和logit。
 
还记得权重的计算公式是1/PStreatment和1/(1- PStreatment)，所以这里numerator = ~1，
 
denominator = 中把我们用于生成PS的协变量加进去。
 
Weights = 这里我们选择sw (standardized weight)，最后的data定义我们的数据来源是statin。
 
/* summary(wt$ipw.weights) */
/* summary(wt$num.mod) */
/* summary(wt$den.mod) */
 
👆这两步是对生成的结果进行总结，ipw.weights总结了所生成的权重的分布，而num.mod和den.mod分别总结了numerator和denominator的模型。
 
图片
图片
/* statin$sw <- wt$ipw.weights */
 
👆将生成的权重转移到原数据statin中去。
 
/* ipw_out <- (svyglm(RFFT ~ Statin, design = svydesign(~ Age + Gender + Ethnicity + Education +
          CVD + DM + Smoking + Hypertension, weights = ~sw, data = statin))) */
/* coef(ipw_out) */
/* confint(ipw_out) */
 
👆Svyglm()函数需要调用”survey”包，这里通过weights = ~sw这一权重选项，可以探究在校正权重后statin对outcome变量RFFT的影响，结果显示statin对RFFT的效应系数为-2.5(95%CI -5.5 – 0.5)，即不显著。
 

图片

通过学习我们可以发现，SAS和R都能比较好地完成倾向性评分匹配和逆概率加权。小编个人认为，相比较而言SAS在这方面的功能更加系统明确，且有清晰的评价手段。而R的实现方式会更加多样化，但也更考验coder灵活调用多种函数的能力。
 
希望这个系列的文章，能在大家开展「真实世界研究」的倾向性评分匹配和逆概率加权时有所帮助~

版权声明：
本网站所有内容来源注明为“梅斯医学”或“MedSci原创”的文字、图片和音视频资料，版权均属于梅斯医学所有。非经授权，任何媒体、网站或个人不得转载，授权转载时须注明来源为“梅斯医学”。其它来源的文章系转载文章，或“梅斯号”自媒体发布的文章，仅系出于传递更多信息之目的，本站仅负责审核内容合规，其内容不代表本站立场，本站不负责内容的准确性和版权。如果存在侵权、或不希望被转载的媒体或个人可与我们联系，我们将立即进行删除处理。
在此留言
评论区 (0)
#插入话题

相关资讯

不平衡数据机器学习的四种处理策略---采用R语言实现

在对不平衡的分类数据集进行建模时，机器学习算法可能并不稳定，其预测结果甚至可能是有偏的，而预测精度此时也变得带有误导性。实际上，经典的统计学建模（如回归），同样也是不稳定的。那么，这种结果是为何发生的呢?到底是什么因素影响了这些算法的表现? 在不平衡的数据中，任一算法都没法从样本量少的类中获取足够的信息来进行精确预测。因此，机器学习算法常常被要求应用在平衡数据集上。那我们该如何处理不平衡数据

做临床预测模型一定要学R语言吗？

如何利用临床预测模型发表高分SCI-问题2

用R语言巧妙处理不平衡数据的方法

在对不平衡的分类数据集进行建模时，机器学**算法可能并不稳定，其预测结果甚至可能是有偏的，而预测精度此时也变得带有误导性那么，这种结果是为何发生的呢？到底是什么因素影响了这些算法的表现？   在不平衡的数据中，任一算法都没法从样本量少的类中获取足够的信息来进行精确预测因此，机器学**算法常常被要求应用在平衡数据集上那我们该如何处理不平衡数据集？本文会介绍一些相关方

R语言：线性回归+样条回归+Nomogram，全搞定

载入我们用到的数据集