AUC柱状图
0.需求
学生听讲座看到一个图，想复现一下，找到了我。。。
他的需求描述如下
目的：
希望将timeROC曲线的9根线的信息展现在柱状图上
疑问：
柱子的高度对应AUC的值；
误差棒的值在哪得到？
上方的统计检验如何添加？
1.数据
rm(list = ls())
load("model_genes_survs.Rdata")
library(survival)
library(timeROC)
result = list()
p = list()
for(i in 1:3){
result[[i]] <-with(survs[[i]], timeROC(T=time,
delta=event,
marker=fp,
cause=1,
times=c(365,1095,1825),
iid = TRUE))
dat = data.frame(fpr = as.numeric(result[[i]]$FP),
tpr = as.numeric(result[[i]]$TP),
time = rep(as.factor(c(365,1095,1825)),each = nrow(result[[i]]$TP)))
library(ggplot2)
p[[i]] = ggplot() +
geom_line(data = dat,aes(x = fpr, y = tpr,color = time),size = 1) +
scale_color_manual(name = NULL,values = c("#92C5DE", "#F4A582", "#66C2A5"),
labels = paste0("AUC of ",c(1,3,5),"-y survival: ",
format(round(result[[i]]$AUC,2),nsmall = 2)))+
geom_line(aes(x=c(0,1),y=c(0,1)),color = "grey")+
theme_bw()+
theme(panel.grid = element_blank(),
legend.background = element_rect(linetype = 1, size = 0.2, colour = "black"),
legend.position = c(0.765,0.125))+
scale_x_continuous(expand = c(0.005,0.005))+
scale_y_continuous(expand = c(0.005,0.005))+
labs(x = "1 - Specificity",
y = "Sensitivity")+
coord_fixed()
}
library(patchwork)
wrap_plots(p)
[图片上传失败...(image-6af064-1691828409852)]
这是上课多个数据批量绘制timeROC曲线的代码。
result是3个数据分别timeROC得到的3个结果，p是三个图。
2.误差棒可以表示什么
误差棒（error bars）是一种用于可视化数据不确定性的图形元素。它们通常用于表示统计数据中的变异程度，如标准差、标准误差、置信区间等。
误差棒一般绘制在数据点或柱状图的上方和下方，以展示数据的范围或不确定性。常见的误差棒类型包括：
1.标准差误差棒（Standard Deviation Error Bars）：基于数据的标准差来表示不确定性，显示了数据分布的散布情况。
2.标准误差误差棒（Standard Error Error Bars）：基于样本均值和样本大小来估计总体均值的不确定性。
3.置信区间误差棒（Confidence Interval Error Bars）：基于统计方法计算出的置信区间来表示不确定性，通常使用较为常见的95%或99%置信水平。
误差棒的长度可以根据需要进行调整，以反映不确定性的大小。较长的误差棒表示较大的不确定性，而较短的误差棒表示较小的不确定性。
绘制误差棒有助于观察数据的变异性和比较不同组之间的差异，同时提供了关于数据可靠性和置信度的信息。这种可视化方法有助于更全面地理解和解释数据分析结果。
3.探索timeROC计算结果
result[[1]]$AUC
## t=365 t=1095 t=1825
## 0.9044373 0.9271197 0.8987427
confint(result[[1]], level = 0.95)$CI_AUC
## 2.5% 97.5%
## t=365 87.58 93.31
## t=1095 89.84 95.59
## t=1825 86.13 93.62
画图所需的数据就是AUC的值和置信区间呗。上面的代码就可以得到。
4.搞数据，画图
dats = lapply(1:3,function(i){
dat = data.frame(x = factor(1:3),
y = result[[i]]$AUC,
confint(result[[i]], level = 0.95)$CI_AUC *0.01)
colnames(dat)[3:4] = c("min","max")
dat$category = paste0("dat",i)
return(dat)
})
dats = do.call(rbind,dats)
dats
## x y min max category
## t=365 1 0.9044373 0.8758 0.9331 dat1
## t=1095 2 0.9271197 0.8984 0.9559 dat1
## t=1825 3 0.8987427 0.8613 0.9362 dat1
## t=3651 1 0.6754528 0.6019 0.7491 dat2
## t=10951 2 0.7570410 0.6986 0.8155 dat2
## t=18251 3 0.7568955 0.6960 0.8178 dat2
## t=3652 1 0.7734769 0.7180 0.8290 dat3
## t=10952 2 0.8281020 0.7802 0.8760 dat3
## t=18252 3 0.8514874 0.8048 0.8982 dat3
数据有咯。画图搞起。
ggplot(dats, aes(x = category, y = y, fill = x)) +
geom_bar(stat = "identity", position = "dodge") +
geom_text(aes(label = round(y,2)), vjust = 4, size = 3,
position = position_dodge(width = 0.9)) +
geom_errorbar(aes(ymin = min, ymax = max),
width = 0.2,
position = position_dodge(width = 0.9)) +
#ggpubr::stat_compare_means(aes(group = x, label = after_stat(p.signif)))+
scale_fill_brewer(palette = "Set2")+
theme_bw()
[图片上传失败...(image-d5b9f7-1691828379528)]
组建比较没画，stat_compare_means只能给出每个组内部的总体p值或者显著性，组内两两比较要手动去搞，这个应该也没啥必要吧。为什么希望一年三年五年生存率的预测准确程度是有差别的啊。
