Different p-value in ggplot2 stat_compare_means and wilcox.test

Question

I try to add p-values to my ggplot using the stat_compare_means function. However, the p-values I get within the ggplot differs from the result of a basic wilcox.test.

I used paired testing in both cases, and also used the wilcoxon test within the ggplot.

I tried to search my question and but couldn't find an exact answer. I updated R (v. 3.5.2), R-Studio (v. 1.1.463), and all packages. In the following I added a few lines of codes with an example. I am new to R and the statistic, so forgive me if I ask in a newbie way.

library("ggplot2")  
library("ggpubr")


c1 <- c( 798.3686, 2560.9974,  688.3051,  669.8265, 2750.6638, 1136.3535,  
         1335.5696, 2347.2777, 1149.1940,  901.6880, 1569.0731 ,3915.6719,  
         3972.0250 ,5517.5016, 4616.6393, 3232.0120, 4020.9727, 2249.4150,  
         2226.4108, 2582.3705, 1653.4801, 3162.2784, 3199.1923, 4792.6118)  
c2 <- c(0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1)  

test <-data.frame(c2,c1)  

test$c2 <- as.factor(test$c2)  

ggplot(test, aes(x=c2, y=c1)) +  
  stat_compare_means(paired = TRUE)  

wilcox.test( test$c1~ test$c2, paired= TRUE)

Result of the stat_compare_means within the ggplot

Result of the Wilcoxon signed rank test:

data: test$c1 by test$c2
V = 0, p-value = 0.0004883
alternative hypothesis: true location shift is not equal to 0

As you can see, the result is p = 0.0025 within the ggplot and p= 0.0004883 with the basic wilcox.test function. Do you know why it's different? And which value is the correct one?

PS: I tried to do the same with ToothGrowths. In that case the result of the stat_compare_means and the wilcox.test show the same results: p = 0.004313. I have no clue why it doesn’t work with my data :/

dipetkov dipetkov · Accepted Answer · 2019-02-23T14:18:37

In one case the p-value is exact and in the other it is a normal approximation.

wilcox.test( test$c1~ test$c2, paired = TRUE, exact = TRUE)
# Wilcoxon signed rank test
# 
# data:  test$c1 by test$c2
# V = 0, p-value = 0.0004883
# alternative hypothesis: true location shift is not equal to 0

wilcox.test( test$c1~ test$c2, paired = TRUE, exact = FALSE)
# Wilcoxon signed rank test with continuity correction
# 
# data:  test$c1 by test$c2
# V = 0, p-value = 0.002526
# alternative hypothesis: true location shift is not equal to 0

According to help(wilcox.test), if the samples contain less than 50 values (as in your case), the exact p-value is computed (unless you specify otherwise).

stat_compare_means has a method.args argument but it doesn't seem to pass the exact = TRUE specification correctly. Instead you can compute the p-value exactly how you want it first and then add it to the plot:

exact_pvalue <-
  wilcox.test( test$c1~ test$c2, paired = TRUE, exact = TRUE) %>%
  # Format the test output as a tibble
  broom::tidy() %>%
  # Format the p-value
  mutate(pval_fmt = format.pval(p.value, digits = 2)) %>%
  # Specify position in (c1, c2) coordinates
  mutate(c1 = 5518, c2 = 0)
exact_pvalue
# A tibble: 1 x 7
#  statistic  p.value method                    alternative pval_fmt    c1    c2
#      <dbl>    <dbl> <chr>                     <chr>       <chr>    <dbl> <dbl>
#1         0 0.000488 Wilcoxon signed rank test two.sided   0.00049   5518     0

ggplot(test, aes(x=c2, y=c1)) +
  geom_text(aes(label = glue::glue("Wilcoxon p = {pval_fmt}")), 
            data = exact_pvalue)

You can generalize this approach to perform multiple tests simulataneously and create a faceted plot at the end. Needs a heavier use of tidyverse magic.

library("tidyverse")

test2 <-
  # Fake data with two subsets to run to test on (in this case the p-value
  # will be the same because the subsets actually contain the same data).
  bind_rows(test, test, .id = "subset") %>%
  # Group by subset and nest the data columns. This creates a "list of
  # tibbles" column called "data".
  group_by(subset) %>%
  nest() %>%
  # Use `purrr::map` to perform the test on each group.
  mutate(wilcox = map(data, ~ wilcox.test(.x$c1 ~ .x$c2,
                                          paired = TRUE, exact = TRUE))) %>%
  # And again `purrr::map` to tidy the test results.
  # Now we have two list columns, one with the data and the other with 
  # the test results
  mutate(wilcox = map(wilcox, broom::tidy))
test2
# A tibble: 2 x 3
# subset data              wilcox
# <chr>  <list>            <list>
#   1 1      <tibble [24 x 2]> <tibble [1 x 4]>
#   2 2      <tibble [24 x 2]> <tibble [1 x 4]>

test2 %>%
  unnest(data) %>%
  ggplot(aes(c1, c2)) +
  # Plot the raw data
  geom_point() +
  # Add the p-value
  geom_text(data = test2 %>% unnest(wilcox),
            # Specify the aestetic mapping so that the p-value is
            # plotted in the top right corner of each plot.
            aes(x = Inf, y = Inf, label = format.pval(p.value, digits = 2)),
            inherit.aes = FALSE, hjust = "inward", vjust = "inward") +
  # Do this for each subset in its own subplot.
  facet_wrap(~ subset)

Different p-value in ggplot2 stat_compare_means and wilcox.test

1 Answers