禁止商业或二改转载,仅供自学使用,侵权必究,如需截取部分内容请后台联系作者!
介绍
pROC::roc
函数能够使用一个指标(predictor)去区分两个或多个分组(response),并计算95%置信区间的原理基于以下几个关键点:
- ROC曲线:ROC曲线是一种图形表示,用于展示分类模型在所有可能的分类阈值下的性能。它通过绘制真正类率(True Positive Rate, TPR)和假正类率(False Positive Rate, FPR)来评估模型的分类能力。
- AUC:曲线下面积(Area Under the Curve, AUC)是一个单一的数字,用于描述ROC曲线下的面积。AUC的值介于0和1之间,其中1表示完美的分类器,0.5表示没有区分能力的分类器(相当于随机猜测)。
-
指标转换:
pROC::roc
函数首先将分组变量(response)中的类别标签转换为二进制形式(例如,"healthy"和"cancer"转换为0和1)。这样,可以使用逻辑回归或其他分类方法来估计预测指标(predictor)的概率。 -
排序和阈值:
pROC::roc
函数根据预测指标的概率对样本进行排序,并计算在每个可能的阈值下模型的TPR和FPR。 -
置信区间:
pROC::roc
函数计算AUC的95%置信区间,这是通过使用非参数方法(如自助法)或正态近似方法来实现的。ci = TRUE
参数指示函数计算这个置信区间。 -
模型拟合:在内部,
pROC::roc
可能使用逻辑回归模型来拟合数据,将预测指标作为预测变量,将分组变量作为响应变量。 -
水平设置:
levels
参数指定了响应变量的类别顺序。这很重要,因为ROC曲线是基于类别的正负性来绘制的。在逻辑回归中,通常将较高级别的类别设置为“成功”或“事件”。 -
统计测试:
pROC::roc
函数还包括对AUC是否统计显著不同于0.5(即随机猜测)的测试,这可以通过pROC::summary.roc
函数获得。
通过这些步骤,pROC::roc
函数提供了一种评估和比较不同预测指标或模型在区分两个或多个分组方面性能的方法。这种方法在医学研究、生物统计学和其他领域中非常常用,尤其是在诊断测试评估和风险预测模型的开发中。
加载R包
knitr::opts_chunk$set(warning = F, message = F)
library(tidyverse)
library(pROC)
library(ggrepel)
导入数据
数据可从以下链接下载(画图所需要的所有数据):数据分析:多诊断指标ROC分析
提取码: 请关注WX公zhong号 *生信学习者* 后台发送 *多指标ROC* 获取提取码
data <- read.csv("./inputdata/54-All_Scores.csv", row.names = 1)
head(data)
sample | type | score.meth | score.delfi | score.meth_reg | ens_preds |
---|---|---|---|---|---|
PGDX16571P | healthy | 0.1147934 | 0.2261810 | 0.08935560 | 0.06348261 |
PGDX16601P | cancer | 0.5178563 | 0.9582597 | 0.41919691 | 0.84425808 |
PGDX16605P | cancer | 0.9422967 | 0.9996974 | 0.91752178 | 0.98553102 |
PGDX17725P | cancer | 0.9174640 | 0.9819394 | 0.69399546 | 0.95644096 |
PGDX17735P | cancer | 0.6345916 | 0.9403198 | 0.06706129 | 0.48634758 |
PGDX17740P | cancer | 0.9886795 | 0.9987350 | 0.99318466 | 0.98982497 |
数据预处理
根据score等指标分别计算它们区分healthy和cancer的能力。这段R代码定义了一个名为get_ROC_CI
的函数,用于计算并汇总不同数据集的ROC曲线分析结果,并最终将结果整合到同一个图形上展示。下面是代码的详细解释:数据分析:多诊断指标ROC分析
1-10. get_ROC_CI
函数接受五个参数:
-
inputdata
:输入的数据框,包含用于计算ROC曲线的数据。 -
index
:用于预测的指标列的名称。 -
group
:包含响应变量(如“健康”或“癌症”)的分组列的名称。 -
group_names
:一个向量,包含group
列中的所有可能的组名。 -
tag
:一个字符串,用于标记结果的类型(如DELFI、Methylation或Ensemble)。
12-13. 将inputdata
中相应的列名替换为"Idx"和"Cmp",以便与pROC::roc
函数的要求一致。
15-21. 使用pROC::roc
函数计算ROC曲线。response
参数设置为分组变量,predictor
设置为预测得分,ci = TRUE
表示需要计算95%置信区间,levels
参数指定了分组变量的顺序。
23-26. 使用pROC::coords
函数和Youden指数确定最佳阈值,这将用于最大化真正例和真负例的总和。
28-33. 再次使用pROC::coords
函数,根据最佳阈值获取最佳性能指标,如敏感性、特异性等。
35-39. 将AUC和95%置信区间格式化为一个字符串,包含标签、AUC值和CI的上下限。
41-47. 创建一个新的数据框(tibble),包含ROC曲线的类型(带有标签的AUC和CI)、敏感性(sensitivities)和特异性(specificities)。返回一个列表,包含ROC对象、最佳阈值、最佳性能指标和ROC数据框。
51-55. 分别对三个不同的数据集(Methylation、DELFI、Ensemble)调用get_ROC_CI
函数,并将结果存储在相应的变量中。
57-65. 将三个结果的数据框合并,并使用dplyr::mutate
和factor
函数调整type
列,以确保所有的类型按照相同的顺序排列。这有助于后续在同一图形上统一展示。
get_ROC_CI <- function(
inputdata,
index,
group,
group_names,
tag) {
# inputdata = data
# index = "score.delfi"
# group = "type"
# group_names = c("healthy", "cancer")
# tag = "DELFI"
colnames(inputdata)[which(colnames(inputdata) == index)] <- "Idx"
colnames(inputdata)[which(colnames(inputdata) == group)] <- "Cmp"
# ROC
roc_res <- pROC::roc(
response = inputdata$Cmp,
predictor = inputdata$Idx,
ci = TRUE,
levels = c("healthy", "cancer"))
# 95% CI
roc_CI <- paste0(tag, ": ", round(roc_res$auc[1], 2),
" (", round(roc_res$ci[1], 2), "-",
round(roc_res$ci[3], 2), ")")
# ROC data.frame
roc_df <- tibble(
type = roc_CI,
sens = roc_res$sensitivities,
spec = roc_res$specificities)
res <- list(roc = roc_res,
cutoff = best_threshold,
best = best_performance,
df = roc_df)
return(res)
}
Methylation_res <- get_ROC_CI(
inputdata = data,
index = "score.meth_reg",
group = "type",
group_names = c("healthy", "cancer"),
tag = "Methylation")
DELFI_res <- get_ROC_CI(
inputdata = data,
index = "score.delfi",
group = "type",
group_names = c("healthy", "cancer"),
tag = "DELFI")
Ensemble_res <- get_ROC_CI(
inputdata = data,
index = "ens_preds",
group = "type",
group_names = c("healthy", "cancer"),
tag = "Ensemble")
merge_res <- rbind(Methylation_res$df, DELFI_res$df, Ensemble_res$df) %>%
dplyr::mutate(type = factor(type, levels = c(unique(Methylation_res$df$type),
unique(DELFI_res$df$type),
unique(Ensemble_res$df$type))))
head(merge_res)
type | sens | spec |
---|---|---|
DELFI: 0.89 (0.84-0.94) | 1 | 0.000000000 |
DELFI: 0.89 (0.84-0.94) | 1 | 0.004761905 |
DELFI: 0.89 (0.84-0.94) | 1 | 0.009523810 |
DELFI: 0.89 (0.84-0.94) | 1 | 0.014285714 |
DELFI: 0.89 (0.84-0.94) | 1 | 0.019047619 |
DELFI: 0.89 (0.84-0.94) | 1 | 0.023809524 |
ROC 曲线图
结果:三种指标对分组Healthy和Cancer的区分ROC曲线。合并后的指标能更有效区分Healthy和Cancer(AUC = 0.93)。
上述三个指标对应的区分Healthy和Cancer的阈值
print(paste("Cutoff of Methylation is", round(as.numeric(Methylation_res$cutoff), 4), " (>:Cancer; <=: Healthy)"))
print(paste("Cutoff of DELFI is", round(as.numeric(DELFI_res$cutoff), 4), " (>:Cancer; <=: Healthy)"))
print(paste("Cutoff of Ensemble is", round(as.numeric(Ensemble_res$cutoff), 4), " (>:Cancer; <=: Healthy)"))
[1] "Cutoff of Methylation is 0.2982 (>:Cancer; <=: Healthy)"
[1] "Cutoff of DELFI is 0.3298 (>:Cancer; <=: Healthy)"
[1] "Cutoff of Ensemble is 0.2289 (>:Cancer; <=: Healthy)"
总结
通过应用pROC::roc
函数,我们对预测指标(predictor)进行了效能分析,旨在评估其区分两个不同分组(response)的能力。通过计算得到的AUC值,我们量化了模型的整体分类性能。进一步地,利用Youden指数,我们确定了最优的区分阈值,以实现在灵敏度和特异性之间最佳平衡。
最终,为了综合比较不同指标的分类效能,我们将它们的ROC曲线汇总在单一图形上进行了展示,直观地呈现了每个指标的AUC值和最优阈值。