0

I have my data below and I would like to plot a multiple box plot with the ANOVA test result shown on the plot.

> combined
   SampleID     chao1 Samples Sgroup Bgroup Duration
1     BSS21 1275.1071   BSS20    BSS      S      20d
2     BSS22 1575.4972   BSS20    BSS      S      20d
3     BSS23 1381.2963   BSS20    BSS      S      20d
4     BSS41 1090.0254   BSS40    BSS      S      40d
5     BSS42 1103.6522   BSS40    BSS      S      40d
6     BSS43 1065.7177   BSS40    BSS      S      40d
7     BSS61 1077.8776   BSS60    BSS      S      60d
8     BSS62 1123.5759   BSS60    BSS      S      60d
9     BSS63 1201.3571   BSS60    BSS      S      60d
10    BSW21  937.0231   BSW20    BSW      W      20d
11    BSW22  970.0462   BSW20    BSW      W      20d
12    BSW23 1070.1560   BSW20    BSW      W      20d
13    BSW41 1894.8606   BSW40    BSW      W      40d
14    BSW42 1825.0271   BSW40    BSW      W      40d
15    BSW43 1869.3494   BSW40    BSW      W      40d
16    BSW61 1332.4078   BSW60    BSW      W      60d
17    BSW62 1354.4593   BSW60    BSW      W      60d
18    BSW63 1365.2961   BSW60    BSW      W      60d
19     BW21 1533.9137    BW20     BW      W      20d
20     BW22 1643.1564    BW20     BW      W      20d
21     BW23 1572.8900    BW20     BW      W      20d
22     BW41 1678.0270    BW40     BW      W      40d
23     BW42 1596.9105    BW40     BW      W      40d
24     BW43 1684.8433    BW40     BW      W      40d
25     BW61 1060.2059    BW60     BW      W      60d
26     BW62 1127.0738    BW60     BW      W      60d
27     BW63 1097.7083    BW60     BW      W      60d
28     SS21 1751.0145    SS20     SS      S      20d
29     SS22 1662.5932    SS20     SS      S      20d
30     SS23 1806.3628    SS20     SS      S      20d
31     SS41 1302.9245    SS40     SS      S      40d
32     SS42 1126.5082    SS40     SS      S      40d
33     SS43 1122.6136    SS40     SS      S      40d
34     SS61 1429.4972    SS60     SS      S      60d
35     SS62 1402.5714    SS60     SS      S      60d
36     SS63 1493.1477    SS60     SS      S      60d
37     SW21 1559.5000    SW20     SW      W      20d
38     SW22 1387.1173    SW20     SW      W      20d
39     SW23 1563.9524    SW20     SW      W      20d
40     SW41 1439.0355    SW40     SW      W      40d
41     SW42 1508.0054    SW40     SW      W      40d
42     SW43 1425.1602    SW40     SW      W      40d
43     SW61 1488.0000    SW60     SW      W      60d
44     SW62 1398.9880    SW60     SW      W      60d
45     SW63 1497.8553    SW60     SW      W      60d
46     W011 1377.8092    W010    WCW      W      10d
47     W012 1304.3725    W010    WCW      W      10d
48     W013 1413.2292    W010    WCW      W      10d
49     W021 1377.8092    W010     BW      W      10d
50     W022 1304.3725    W010     BW      W      10d
51     W023 1413.2292    W010     BW      W      10d
52     W031 1377.8092    W010     SW      W      10d
53     W032 1304.3725    W010     SW      W      10d
54     W033 1413.2292    W010     SW      W      10d
55     W041 1377.8092    W010    BSW      W      10d
56     W042 1304.3725    W010    BSW      W      10d
57     W043 1413.2292    W010    BSW      W      10d
58     W051 1377.8092    W010     SS      W      10d
59     W052 1304.3725    W010     SS      W      10d
60     W053 1413.2292    W010     SS      W      10d
61     W061 1377.8092    W010    BSS      W      10d
62     W062 1304.3725    W010    BSS      W      10d
63     W063 1413.2292    W010    BSS      W      10d
64    WCW21 1246.5794   WCW20    WCW      W      20d
65    WCW22 1249.2180   WCW20    WCW      W      20d
66    WCW23 1134.3462   WCW20    WCW      W      20d
67    WCW41 1074.9192   WCW40    WCW      W      40d
68    WCW42  887.7191   WCW40    WCW      W      40d
69    WCW43  990.3733   WCW40    WCW      W      40d
70    WCW61  864.2727   WCW60    WCW      W      60d
71    WCW62  934.5111   WCW60    WCW      W      60d
72    WCW63  801.5696   WCW60    WCW      W      60d

I tried the ggpubr package in r and apply the compare_means function. For the anova test, I got:

> compare_means(chao1~Duration,data=combined,method="anova",group.by = "Sgroup")

# A tibble: 6 x 7
  Sgroup .y.             p      p.adj p.format p.signif method
  <fct>  <chr>       <dbl>      <dbl> <chr>    <chr>    <chr> 
1 BSS    chao1 0.00394     0.0079     0.00394  **       Anova 
2 BSW    chao1 0.000000164 0.00000098 1.6e-07  ****     Anova 
3 BW     chao1 0.00000283  0.000014   2.8e-06  ****     Anova 
4 SS     chao1 0.0000996   0.0004     1.0e-04  ****     Anova 
5 SW     chao1 0.160       0.16       0.16005  ns       Anova 
6 WCW    chao1 0.000118    0.0004     0.00012  ***      Anova

Which is exactly what I wanted. Then I make the plot:

ggboxplot(combined,x="Samples", y="chao1",palette = "jco", add = "jitter",short.panel.labs = FALSE)+facet_wrap(~Sgroup,scales = "free")

and got:

enter image description here

I have tried so hard by the stat_compare_means function but I still cannot add my anova test result on the plot

Any suggestion will be helped!

Lennon Lee
  • 194
  • 1
  • 14

1 Answers1

1

I've had problems with facet_wrap in the past, it maybe that its not supported. Hopefully I'm proven wrong! If you want an alternative, you can use the cowplot package to build the grids yourself. It's more work, but allows more flexibility perhaps. Here is an example with the data used when you run ?stat_compare_means

# packages
library(cowplot)
library(ggpubr)

# data from ggpubr
data("ToothGrowth")
head(ToothGrowth)

# define groups, build plots in a list and subset
groups <- unique(ToothGrowth$supp)

pl <- lapply(groups, function(g){
  p <- ggboxplot(ToothGrowth[ToothGrowth$supp==g,], x = "dose", y = "len",
                 color = "dose", palette = "npg") + 
    stat_compare_means()
})

# use cowplot to bring together
cowplot::plot_grid(plotlist = pl, ncol = 2)
Jonny Phelps
  • 2,687
  • 1
  • 11
  • 20