图中是三个大组,内部用渐变色画亚组。
红色和蓝色大组分别画了半透明的红色和蓝色误差线。
这个图，可以比较组内的情况，然后还可以合并组内的数据，再比组间。还有这种可以直接实现出图的代码。
文献中数据产生步骤
即复现图i和j
多组柱状图(mean ± se; se = sd/sqrt(n))，某些组同时存在亚组，并进行组间比较和亚组内比较。有些对数据的描述不允许使用se（视觉上感觉离散度更小），要使用sd。该绘图代码在出图前需要不断调整坐标或文字位置，需要在理解的基础上加以使用。
##计算se
se <- function(x) sd(x)/sqrt(length(x))
###multicenter_validation_*.txt，每个亚组的观测值对应图中的一个bar，放到一个单独的文件里，只需提供diff这一列。
prefix <- "multicenter_validation_" # 输入文件的文件前缀
# 加载图中绿色大组对应的数据
full1 <- read.table(paste0(prefix,"full_NJDTH_UI.txt"),sep = "\t",check.names = F,stringsAsFactors = F,header = T,row.names = NULL)
# 加载图中红色大组对应的数据
semi1 <- read.table(paste0(prefix,"semi_NJDTH_UI.txt"),sep = "\t",check.names = F,stringsAsFactors = F,header = T,row.names = NULL)
semi2 <- read.table(paste0(prefix,"semi_PZTH_UI.txt"),sep = "\t",check.names = F,stringsAsFactors = F,header = T,row.names = NULL)
semi3 <- read.table(paste0(prefix,"semi_SYCCYHC_UI.txt"),sep = "\t",check.names = F,stringsAsFactors = F,header = T,row.names = NULL)
# 加载图中蓝色大组对应的数据
manu1 <- read.table(paste0(prefix,"manual_NJDTH_UI.txt"),sep = "\t",check.names = F,stringsAsFactors = F,header = T,row.names = NULL)
manu2 <- read.table(paste0(prefix,"manual_NJDTH_GE.txt"),sep = "\t",check.names = F,stringsAsFactors = F,header = T,row.names = NULL)
manu3 <- read.table(paste0(prefix,"manual_NJDTH_Philips.txt"),sep = "\t",check.names = F,stringsAsFactors = F,header = T,row.names = NULL)
manu4 <- read.table(paste0(prefix,"manual_JSPPH_Siemens.txt"),sep = "\t",check.names = F,stringsAsFactors = F,header = T,row.names = NULL)
manu5 <- read.table(paste0(prefix,"manual_GCPH_Toshiba.txt"),sep = "\t",check.names = F,stringsAsFactors = F,header = T,row.names = NULL)
# 1. 计算各亚组均值（注意均值的顺序是后期调整的，根据均值大小[柱子高低]）
avg <- c(mean(full1$diff),0, # 这里的0是为了柱状图（full和semi）之间的空隙
mean(semi1$diff),
mean(semi2$diff),
mean(semi3$diff),0, # 这里的0是为了柱状图（semi和manu）之间的空隙
mean(manu3$diff),
mean(manu2$diff),
mean(manu1$diff),
mean(manu4$diff),
mean(manu5$diff))
# 2. 计算个亚组标准误（注意柱子的顺序要和均值一致，这里的manu组为32145）
var <- c(se(full1$diff),0,
se(semi1$diff),
se(semi2$diff),
se(semi3$diff),0,
se(manu3$diff),
se(manu2$diff),
se(manu1$diff),
se(manu4$diff),
se(manu5$diff))
# 3. 计算有多个亚组的组均值
avg_semi <- mean(c(semi1$diff,semi2$diff,semi3$diff))
avg_manu <- mean(c(manu1$diff,manu2$diff,manu3$diff,manu4$diff,manu5$diff))
# 4. 计算有多个亚组的组se
var_semi <- se(c(semi1$diff,semi2$diff,semi3$diff))
var_manu <- se(c(manu1$diff,manu2$diff,manu3$diff,manu4$diff,manu5$diff))
# 5. semi亚组的方差分析（可换成kruskal）
tmp <- data.frame(diff = c(semi1$diff,semi2$diff,semi3$diff),
class = as.factor(rep(c("A","B","C"), c(nrow(semi1),nrow(semi2),nrow(semi3)))))
saov <- summary(aov(diff~class,tmp)) # intersemi p=6.04e-09
# 6. manual亚组的方差分析（可换成kruskal）
tmp <- data.frame(diff = c(manu1$diff,manu2$diff,manu3$diff,manu4$diff,manu5$diff),
class = as.factor(rep(c("A","B","C","D","E"), c(nrow(manu1),nrow(manu2),nrow(manu3),nrow(manu4),nrow(manu5)))))
maov <- summary(aov(diff~class,tmp)) # intermanul p<2e-16
# 7. 组间t检验（可换成wilcox非参检验，未校正）
fvss.t <- t.test(full1$diff,c(semi1$diff,semi2$diff,semi3$diff))$p.value # full vs semi p=2.194e-06
fvsm.t <- t.test(full1$diff,c(manu1$diff,manu2$diff,manu3$diff,manu4$diff,manu5$diff))$p.value # full vs manul p=1.557547e-18
svsm.t <- t.test(c(semi1$diff,semi2$diff,semi3$diff),c(manu1$diff,manu2$diff,manu3$diff,manu4$diff,manu5$diff))$p.value # semi vs manul p=0.1779296
# 设置颜色 #
darkblue <- "blue"
seagreen <- "green"
sun <- "red"
red <- "darkred"
# 编写label，注意和均值的柱子一致
lab <- c("Full-NJDTH-UI","",
"Semi-NJDTH-UI","Semi-PZTH-UI","Semi-SYCCYHC-UI","",
"Manual-NJDTH-Philips","Manual-NJDTH-GE","Manual-NJDTH-UI","Manual-JSPPH-Siemens","Manual-GCPH-Toshiba")
# 绘制基本柱状图
pdf("fancybar.pdf",width = 6,height = 5)
par(bty="o", mgp = c(2,0.5,0), mar = c(7.1,4.1,2.1,4.1),tcl = -.25,las = 1) # 基本设置
par(xpd = T) # 允许图像超出（一般是为了画图例用）
bar <- barplot(avg,
border = F, # 柱子不显示轮廓
ylab = "Scanning length error (mm)", # y轴标签
ylim = c(0,80), # 高度要比最大的均值要大一些，给添加P值预留空间
yaxt = "n", # 不显示y轴
xaxt = "n", # 不显示x轴
col = c(seagreen,NA, # 注意颜色的NA是均值为0的部分，为空白柱
sun, ggplot2::alpha(sun,0.6), ggplot2::alpha(sun,0.3), NA, # 注意颜色的渐变透明
darkblue, ggplot2::alpha(darkblue,0.8), ggplot2::alpha(darkblue,0.6), ggplot2::alpha(darkblue,0.4), ggplot2::alpha(darkblue,0.2)))
基本柱状图
axis(side = 2,at = seq(0,60,10)) # 手动添加y轴
添加y轴
# 添加垂直的误差线
segments(bar, # 由bar对象确定位置，也就是每个柱子的中心点
avg - var, # 下侧最低点，为mean- se
bar,
avg + var, # 顶部最高点，为mean + se
lwd = 1.5)
垂直误差线
# 添加水平的误差线
arrows(bar, # 同上位置添加水平线
avg - var,
bar,
avg + var,
lwd = 1.5, # 线宽
angle = 90,
code = 3, length = 0.05)
水平误差线
# 根据bar对象添加x轴标签
text(x = bar,
y = par("usr")[3] - 1, # 在图像垂直最低点再往下1单位添加x轴标签
srt = 45,
adj = 1, labels =lab, xpd = TRUE, cex = 0.6)
x轴标签
# 添加图例
par(xpd = T)
legend(par("usr")[2]-0.5, # 水平最右侧向内0.5
par("usr")[4], # 垂直最顶上
legend = c("Full","Semi","Manual"), # 原文中带有渐变的图例是后期AI的，因为考虑到柱子本身存在渐变
fill = c(seagreen,sun,darkblue),
cex=0.8,
border=NA,
y.intersp=1,
x.intersp=0.2,
bty = "n") # 不要图例边框
图例
# 添加semi亚组的背景透明区块（这里要好好体会坐标）
par(new = T,xpd = F) # 新添加图层，允许元素超过边界
rect(xleft = bar[3]-0.5,# 红色semi组的左侧外缘在bar对象x轴第三个点-0.5的位置（每个柱子宽1，左右各0.5）
ybottom = 0, # y轴底为0
xright = bar[5]+0.5, # 同理x轴右侧在第五个点向右再+0.5
ytop = avg_semi, # y顶就是组均值
col = ggplot2::alpha(sun,0.1), # 透明色
border = NA)
添加semi亚组的背景透明区块
# 同理添加误差线（线稍微粗一些，颜色透明）
segments(bar[4], # 在bar的第四个点，也就是红色亚组最中间的点上添加
avg_semi - var_semi,
bar[4],
avg_semi + var_semi,
lwd = 4,
col = ggplot2::alpha(red,0.4))
arrows(bar[4],
avg_semi - var_semi,
bar[4],
avg_semi + var_semi,
lwd = 4, angle = 90,
code = 3, length = 0.05,col = ggplot2::alpha(sun,0.4),lty = 2)
添加误差线
网友评论