如何在 R 中绘制带有条件 p 值和置信区间的 ANOVA？

Question

我有 df1：

    Rate       Dogs        MHI_2018  Points      Level AGE65_MORE P_Elderly
1  0.10791173  0.00000000    59338 236.4064         C       8653  15.56267
2  0.06880040  0.00000000    57588 229.4343         C      44571  20.44335
3  0.08644537  0.00000000    50412 200.8446         C      10548  18.23651
4  0.29591635  0.00000000    29267 116.6016         A       1661  16.38390
5  0.05081301  0.00000000    37365 148.8645         B       3995  20.29980
6  0.02625200  0.00000000    45400 180.8765         D      20247  17.71748
7  0.80321285  0.02974862    39917 159.0319         D       6562  19.52105
8  0.07682852  0.00000000    42132 167.8566         D       5980  22.97173
9  0.18118814  0.00000000    47547 189.4303         B       7411  16.78482
10 0.07787555  0.00000000    39907 158.9920         B       2953  22.99665
11 0.15065913  0.00000000    39201 156.1793         C       2751  20.72316
12 0.33362247  0.00000000    46495 185.2390         B       2915  19.45019
13 0.03652168  0.00000000    49055 195.4382         B      10914  19.92988
14 0.27998133  0.00000000    42423 169.0159         A       2481  23.15446
15 0.05407451  0.00000000    40203 160.1713         A       7790  21.06202
16 0.07233796  0.00000000    39057 155.6056         A       2629  19.01765
17 0.08389061  0.00000000    45796 182.4542         B      15446  18.51106
18 0.05220569  0.00000000    34035 135.5976         B       6921  18.06578
19 0.05603418  0.00000000    39491 157.3347         B      12322  17.26133
20 0.15875536  0.00000000    60367 240.5060         C      12400  15.14282

有

    AOV <- aov(Rate~Level, data = df)
    TukeyHSD(AOV)
$Level
            diff        lwr       upr     p adj
B-A -0.066558621 -0.3783957 0.2452784 0.9272012
C-A -0.061063140 -0.4026635 0.2805372 0.9551663
D-A  0.126520253 -0.2624089 0.5154494 0.7890519
C-B  0.005495482 -0.2848090 0.2958000 0.9999404
D-B  0.193078874 -0.1516699 0.5378277 0.4049948
D-C  0.187583392 -0.1843040 0.5594708 0.4923479

我现在想用均值和置信区间绘制这些数据的图。如果 p_adj

Answer 1

一种解决方案是使用ggsignif 包，但首先您需要准备TukeyHSD 的输出，以便在ggsignif 中使用：

AOV <- aov(Rate~Level, data = df)
t <-as.data.frame(TukeyHSD(AOV)$Level)

library(tidyverse)

MAX <-df %>% group_by(Level) %>% summarise(Max = max(Rate))

T1 <- t %>% rownames_to_column("Group") %>%
  mutate(Start = sub("^(.).*","\\1",Group),
         End = sub(".*(.)$","\\1",Group)) %>%
  left_join(.,MAX, by = c("Start" = "Level")) %>%
  left_join(.,MAX, by = c("End" = "Level")) %>%
  mutate(End = factor(End)) %>%
  rowwise() %>%mutate(ypos = max(Max.x, Max.y)*(1+0.25*as.numeric(End)))

Source: local data frame [6 x 10]
Groups: <by row>

# A tibble: 6 x 10
  Group     diff    lwr   upr `p adj` Start End   Max.x Max.y  ypos
  <chr>    <dbl>  <dbl> <dbl>   <dbl> <chr> <fct> <dbl> <dbl> <dbl>
1 B-A   -0.0666  -0.378 0.245   0.927 B     A     0.334 0.296 0.417
2 C-A   -0.0611  -0.403 0.281   0.955 C     A     0.159 0.296 0.370
3 D-A    0.127   -0.262 0.515   0.789 D     A     0.803 0.296 1.00 
4 C-B    0.00550 -0.285 0.296   1.00  C     B     0.159 0.334 0.500
5 D-B    0.193   -0.152 0.538   0.405 D     B     0.803 0.334 1.20 
6 D-C    0.188   -0.184 0.559   0.492 D     C     0.803 0.159 1.41 

现在，您可以根据 T1 数据集绘制数据并添加重要性：

library(ggsignif)
library(ggplot2)

ggplot(df, aes(x = Level, y = Rate))+
  geom_jitter(width = 0.2)+
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", width = 0, color = "red") +
  stat_summary(fun = "mean", geom = "errorbar",  aes(ymax = ..y.., ymin = ..y..), col = "red", width = 0.5) +
  geom_signif(data = subset(T1,`p adj` <0.5), manual = TRUE,
              aes(xmax = End, xmin = Start, y_position= ypos, annotations = round(`p adj`,3)))