Animated violin/boxplots plots with p-values for paired Wilcoxon tests

Question

I am trying to add p values for paired wilcox test in R. I am using the following code. The code below creates violins (density distributions) of an outcome reading (bicep) for two diets (treatment). These violins are animated over time 1, time 2, time 3. And the top of the graph prints p-values. I would like these p-values to be paired values such that

For Diet 'a' Bicep reading at time 2 is compared to time 1, and bicep reading at time 3 is compared to time 1.

And the same for Diet 'b'. So, there should be two separate p-values printed on top of the violins at time 2 and time 3. Indicating paired tests (Time 2 vs Time 1 and Time 3 vs Time 1) for both Diet 'a' and Diet 'b'.

What should be the correct code for this test? I have tried something here below based on a suggestion I got yesterday, but I ran into an error. I also think the code below just does paired tests for Time 2 vs Time 1, and Time 3 vs Time 2. Which is not what I want.

Thanks for reading.

 library(ggplot2)
 library(ggpubr)
 library(gganimate)
 library(tidyverse)

Example Data

 structure(list(code = c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 
4L, 4L, 5L, 5L, 5L, 6L, 6L, 6L, 7L, 7L, 7L, 8L, 8L, 8L), diet = c("a", 
"a", "a", "b", "b", "b", "a", "a", "a", "b", "b", "b", "a", "a", 
"a", "b", "b", "b", "a", "a", "a", "b", "b", "b"), time = c(1L, 
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
3L, 1L, 2L, 3L, 1L, 2L, 3L), bicep = c(8L, 7L, 7L, 9L, 9L, 9L, 
11L, 10L, 9L, 11L, 11L, 12L, 12L, 11L, 10L, 9L, 9L, 9L, 12L, 
10L, 8L, 12L, 12L, 12L)), class = "data.frame", row.names = c(NA, 
-24L))

Example Code

example3 %>%
  group_by(time) %>% 
  mutate(p=pairwise.wilcox.test(example3$bicep, interaction(example3$diet, example3$time), p.adjust.method = "none")$p.value,
         max=max(bicep, na.rm = T)) %>% 
  ggplot() + 
  geom_violin(aes(x=diet, y=bicep, fill=diet)) +
  geom_text(data = . %>% distinct(p, max, time), 
            aes(x=1.5, y = max+.5, label=as.character(round(p,2))),
            size=12) +
  transition_time(time) +
  ease_aes('linear')

This is the error I get

 Error in mutate_impl(.data, dots) : 
  Column `p` must be length 8 (the group size) or one, not 25
In addition: There were 15 warnings (use warnings() to see them)

Answer 1

If I understand you problem right then there is very easy solution. You got this error because of wrong syntax inside mutate. There is no need to call for values with $ when you using mutate and pipes %>%:

This code gives a desired animation with minor warnings:

example3 %>%
  group_by(time) %>% 
  mutate(p=pairwise.wilcox.test(bicep, interaction(diet, time), p.adjust.method = "none")$p.value,
         max=max(bicep, na.rm = T)) %>% 
  ggplot() + 
  geom_violin(aes(x=diet, y=bicep, fill=diet)) +
  geom_text(data = . %>% distinct(p, max, time), 
            aes(x=1.5, y = max+.5, label=as.character(round(p,2))),
            size=12) +
  transition_time(time) +
  ease_aes('linear')

UPDATE

In case of independent p-values you just need to add facet_wrap(), for example. It seems to be the easiest:

example3 %>%
  group_by(time) %>% 
  mutate(p = pairwise.wilcox.test(bicep, interaction(diet, time), p.adjust.method = "none")$p.value,
         max = max(bicep, na.rm = T)) %>%
  ggplot() + 
  geom_violin(aes(x = diet, y = bicep, fill = diet)) +
  geom_text(data = . %>% distinct(p, max, time), 
            aes(x = 1, y = max+.5, label = as.character(round(p,2))),
            size = 12) +
  facet_wrap(~diet, scales = "free_x") + # add facets
  theme(legend.position = "none") +
  transition_time(time) +
  ease_aes('linear')