R中的下划线图

debugcn 发表于 Dev

NM_

介绍和当前工作

[注意：对于那些感兴趣的人，我在最后提供了用于重现示例的代码。]

我有一些数据，并且进行了ANOVA分析，并获得了Tukey的成对比较：

model1 = aov(trt ~ grp, data = df)
anova(model1)

> TukeyHSD(model1)
          diff         lwr       upr     p adj
B-A 0.03481504 -0.40533118 0.4749613 0.9968007
C-A 0.36140489 -0.07874134 0.8015511 0.1448379
D-A 1.53825179  1.09810556 1.9783980 0.0000000
C-B 0.32658985 -0.11355638 0.7667361 0.2166301
D-B 1.50343674  1.06329052 1.9435830 0.0000000
D-C 1.17684690  0.73670067 1.6169931 0.0000000

我还可以绘制Tukey的成对比较

> plot(TukeyHSD(model1))

从Tukey的置信区间和图可以看出A-B，B-C并且A-C没有显着差异。

问题

我被要求创建一个称为“下划线”的东西，其描述如下：

我们在实线上绘制组均值，并在组均值之间绘制一条线段，以表明这两个特定组之间没有显着差异。

获得手段并不困难：

> aggregate(df$trt ~ df$grp, FUN = mean)
  df$grp   df$trt
1      A 2.032086
2      B 2.066901
3      C 2.393491
4      D 3.570338

期望的输出

使用本示例中的数据，所需的图应如下图所示：

组之间存在没有明显差异的线段（即A-B，B-C和之间的线段，A-C如Tukey所示）。

注意：请注意，上面的图未按比例绘制，它是在主题演讲中创建的，仅用于说明目的。

有没有办法使用R（使用基本R或诸如的库ggplot2）来获得上述的“下划线图” ？

编辑

这是我用来创建以上示例的代码：

library(data.table)

set.seed(3)
A = runif(20, 1,3)
A = data.frame(A, rep("A", length(A)))
B = runif(20, 1.25,3.25)
B = data.frame(B, rep("B", length(B)))
C = runif(20, 1.5,3.5)
C = data.frame(C, rep("C", length(C)))
D = runif(20, 2.75,4.25)
D = data.frame(D, rep("D", length(D)))

df = list(A, B, C, D)
df = rbindlist(df)

colnames(df) = c("trt", "grp")

eipi10

这是下划线图的ggplot版本。我们将加载tidyverse包，它会加载ggplot2，dplyr并从tidyverse其他几个包。我们创建一个系数数据框以绘制组名，系数值和垂直段，并创建一个非有效对的数据框以生成水平下划线。

library(tidyverse)

model1 = aov(trt ~ grp, data=df)

# Get coefficients and label coefficients with names of levels
coefs = coef(model1)
coefs[2:4] = coefs[2:4] + coefs[1]
names(coefs) = levels(model1$model$grp)

# Get non-significant pairs
pairs = TukeyHSD(model1)$grp %>% 
  as.data.frame() %>% 
  rownames_to_column(var="pair") %>% 
  # Keep only non-significant pairs
  filter(`p adj` > 0.05) %>% 
  # Add coefficients to TukeyHSD results
  separate(pair, c("pair1","pair2"), sep="-", remove=FALSE) %>% 
  mutate(start = coefs[match(pair1, names(coefs))],
         end = coefs[match(pair2, names(coefs))]) %>% 
  # Stagger vertical positions of segments
  mutate(ypos = seq(-0.03, -0.04, length=3))

# Turn coefs into a data frame
coefs = enframe(coefs, name="grp", value="coef")

ggplot(coefs, aes(x=coef)) +
  geom_hline(yintercept=0) +
  geom_segment(aes(x=coef, xend=coef), y=0.008, yend=-0.008, colour="blue") +
  geom_text(aes(label=grp, y=0.011), size=4, vjust=0) +
  geom_text(aes(label=sprintf("%1.2f", coef)), y=-0.01, size=3, angle=-90, hjust=0) +
  geom_segment(data=pairs, aes(group=pair, x=start, xend=end, y=ypos, yend=ypos),
               colour="red", size=1) +
  scale_y_continuous(limits=c(-0.05,0.04)) +
  theme_void()