# 1. Load required R packages
suppressPackageStartupMessages(library(ggpubr))
suppressPackageStartupMessages(library(rstatix))

# 2. Data preparation
df <- ToothGrowth
df$dose <- factor(df$dose)

# 3. Statistical tests
res.stats <- df %>%
  group_by(dose) %>%
  t_test(len ~ supp) %>%
  adjust_pvalue() %>%
  add_significance()
res.stats

## # A tibble: 3 x 11
##   dose  .y.   group1 group2    n1    n2 statistic    df       p   p.adj p.adj.signif
##   <fct> <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl>   <dbl>   <dbl> <chr>       
## 1 0.5   len   OJ     VC        10    10    3.17    15.0 0.00636 0.0127  *           
## 2 1     len   OJ     VC        10    10    4.03    15.4 0.00104 0.00312 **          
## 3 2     len   OJ     VC        10    10   -0.0461  14.0 0.964   0.964   ns

# 4. Create a stacked bar plot, add "mean_se" error bars
p <- ggbarplot(
  df, x = "dose", y = "len", add = "mean_se",
   color = "supp", palette = "jco"
  )

# 5. Add p-values to the bar plot using ggpubr verbs
p + stat_pvalue_manual(
  res.stats, x = "dose", y.position = c(30, 45, 60),
  label = "p.adj.signif"
)

# or use ggplot2 verbs
p + geom_text(
  aes(x = dose, y = c(25, 42, 60), label = p.adj.signif),
  data = res.stats
)