How to add a Scheffe Test to ggplot?

Question

I´m new in R and I taught it myself so please be indulgent with me :)

I want to include the Scheffe-Test to my ggplot. I have tried with the T-Test and the Anova that worked for me. Now that I want to change it to a Scheffe-Test doesnt work for me anymore. For me the reason seems to be the missing group1 and group2 columns but I don´t know how to add them. Can anybody help? I will try to give an minimal working example:

My library:

library(tidyverse)
library(magrittr)
library(gtsummary)
library(ggpubr)
library(rstatix)
library(ggplot2)
library(DescTools)

The data:

> data_clean_bef_tbl
# A tibble: 11 × 3
   time   TEER Assay 
      
 1 -7      182 before
 2 -7      194 before
 3 -4      182 before
 4 -4      190 before
 5 -4      219 before
 6 0       241 before
 7 0       271 before
 8 0       246 before
 9 3       247 before
10 3       284 before
11 3       259 before

The ggplot:

ggplot_vor <- ggplot(data_clean_bef_tbl, aes(x = time, y = TEER)) +
  geom_boxplot(aes(fill = time),
               outlier.shape = NA) +
  scale_fill_manual(values = c("#FFF176", "#AED581", "#4DB6AC", "#7986CB")) +
  xlab("Seeding time of HT29-MTX-E12 cells compared to Caco-2 cells [days]") +
  ylab("TEER [Ω*qcm]") +
  theme_minimal() +
  geom_point(mapping = aes(shape=time),
             position = position_jitterdodge()) +
  stat_summary(fun="mean", fun.min = "min", fun.max = "max", 
               color = "#D32F2F") +
  stat_compare_means(method = "anova", label.y = 410)

The Anova- and Scheffe-Test:

anova_bef <- aov(TEER ~ time, data_clean_bef_tbl)
summary(anova_bef)
ScheffeTest(anova_bef, group = TRUE)

The T-Test that works fine for me:

stat_ttest_vor <- data_clean_bef_tbl %>%
  t_test(TEER ~ time, var.equal = FALSE) %>% 
  adjust_pvalue(method = "bonferroni") %>%
  add_significance("p.adj") %>%
  add_xy_position(x = "time", dodge = 0.8)

ggplot + t-test + Anova

Boxplot_vor <- ggplot_vor +
  stat_pvalue_manual(stat_ttest_vor,
                     label = "p.adj.signif",
                     step.increase = 0.02) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
  stat_compare_means(method = "anova", label.y = 410)
Boxplot_vor

ggplot with Anova and T-test, but I want the same with the Scheffe-Test

When I try the same with this code:

Boxplot_vor <- ggplot_vor +
  stat_pvalue_manual(Scheffe_bef,
                     label = "pval",
                     step.increase = 0.02) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
  stat_compare_means(method = "anova", label.y = 410)
Boxplot_vor

I get the error message:

Error in asserttat_group_columns_exists(data) : data should contain group1 and group2 columns

I´ve tried to add group1 and group2 columns manually to the Scheffe-test but since it´s not a tibble, this isn´t as easy as I thought.

user_14282714 · Answer

An option could be using the Plot2WayANOVA function from the CGPfunctions, documenation. First of all you should have to classes in your Assay column, so I added some random values for "after" class. According to the documentation this function:

Takes a formula and a dataframe as input, conducts an analysis of variance prints the results (AOV summary table, table of overall model information and table of means) then uses ggplot2 to plot an interaction graph (line or bar) . Also uses Brown-Forsythe test for homogeneity of variance. Users can also choose to save the plot out as a png file.

Where you have the argument potshot.method where you can use the Scheffe test:

A character string, one of "hsd", "bonf", "lsd", "scheffe", "newmankeuls", defining the method for the pairwise comparisons. (Default: "scheffe").

The output of this function returns:

A list with 5 elements which is returned invisibly. These items are always sent to the console for display but for user convenience the function also returns a named list with the following items in case the user desires to save them or further process them - $ANOVATable, $ModelSummary, $MeansTable, $PosthocTable, $BFTest, and $SWTest. The plot is always sent to the default plot device

So you can check in the $Posthoctable your results of the Scheffe test. Here is a reproducible example:

data_clean_bef_tbl <- read.table(text="   time   TEER Assay 
 -7      182 before
 -7      194 before
 -4      182 before
 -4      190 before
 -4      219 before
 0       241 before
 0       271 before
 0       246 before
 3       247 before
 3       284 before
 3       259 before
 -7      183 after
 -7      195 after
 -4      186 after
 -4      192 after
 -4      216 after
 0       242 after
 0       274 after
 0       246 after
 3       248 after
 3       283 after
 3       252 after", header = TRUE)

library(ggplot2)
library(CGPfunctions)
library(DescTools)
library(pwr)
anova_bef <- aov(TEER ~ time, data_clean_bef_tbl)
summary(anova_bef)
#>             Df Sum Sq Mean Sq F value   Pr(>F)    
#> time         1  20920   20920   72.94 4.18e-08 ***
#> Residuals   20   5736     287                     
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Plot2WayANOVA(TEER ~ time*Assay, 
              data = data_clean_bef_tbl, 
              xlab = "Seeding time of HT29-MTX-E12 cells compared to Caco-2 cells [days]",
              ylab = "TEER [Ω*qcm]",
              plottype = "bar", 
              posthoc.method = "scheffe",
              overlay.type = "box",
              ggplot.component = theme_minimal())
#> 
#> Converting time to a factor --- check your results
#> 
#> Converting Assay to a factor --- check your results
#> 
#>              --- WARNING! ---
#>      You have one or more cells with less than 3 observations.
#> # A tibble: 8 × 3
#>   time  Assay  count
#>      
#> 1 -7    after      2
#> 2 -7    before     2
#> 3 -4    after      3
#> 4 -4    before     3
#> 5 0     after      3
#> 6 0     before     3
#> 7 3     after      3
#> 8 3     before     3
#> 
#>              --- WARNING! ---
#>      You have an unbalanced design. Using Type II sum of 
#>             squares, to calculate factor effect sizes eta and omega.
#>             Your two factors account for 0.851 of the type II sum of 
#>             squares.
#> term       |     sumsq |   meansq | df | statistic | p.value | etasq | partial.etasq | omegasq | partial.omegasq | epsilonsq | cohens.f | power
#> -----------------------------------------------------------------------------------------------------------------------------------------------
#> time       | 22665.697 | 7555.232 |  3 |    26.594 |  < .001 | 0.850 |         0.851 |   0.810 |           0.777 |     0.818 |    2.387 | 1.000
#> Assay      |     0.182 |    0.182 |  1 |     0.001 |   0.980 | 0.000 |         0.000 |  -0.011 |          -0.048 |    -0.011 |    0.007 | 0.050
#> time:Assay |    13.152 |    4.384 |  3 |     0.015 |   0.997 | 0.000 |         0.003 |  -0.031 |          -0.155 |    -0.031 |    0.058 | 0.053
#> Residuals  |  3977.333 |  284.095 | 14 |           |         |       |               |         |                 |           |          |
#> 
#> Table of group means
#> # A tibble: 8 × 15
#> # Groups:   time [4]
#>   time  Assay  TheMean TheSD TheSEM CIMuliplier LowerBoundCI UpperBoundCI
#>                                  
#> 1 -7    after     189   8.49   6          12.7          113.         265.
#> 2 -7    before    188   8.49   6          12.7          112.         264.
#> 3 -4    after     198  15.9    9.17        4.30         159.         237.
#> 4 -4    before    197  19.5   11.2         4.30         149.         245.
#> 5 0     after     254  17.4   10.1         4.30         211.         297.
#> 6 0     before    253. 16.1    9.28        4.30         213.         293.
#> 7 3     after     261  19.2   11.1         4.30         213.         309.
#> 8 3     before    263. 18.9   10.9         4.30         216.         310.
#> # … with 7 more variables: LowerBoundSEM , UpperBoundSEM ,
#> #   LowerBoundSD , UpperBoundSD , N , LowerBound ,
#> #   UpperBound 
#> 
#> Post hoc tests for all effects that were significant
#> 
#>   Posthoc multiple comparisons of means: Scheffe Test 
#>     95% family-wise confidence level
#> 
#> $time
#>            diff    lwr.ci    upr.ci   pval    
#> -4--7  9.000000 -38.85861  56.85861 0.9976    
#> 0--7  64.833333  16.97473 112.69194 0.0048 ** 
#> 3--7  73.666667  25.80806 121.52527 0.0015 ** 
#> 0--4  55.833333  13.02729  98.63937 0.0067 ** 
#> 3--4  64.666667  21.86063 107.47271 0.0018 ** 
#> 3-0    8.833333 -33.97271  51.63937 0.9958    
#> 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Testing Homogeneity of Variance with Brown-Forsythe
#> Brown-Forsythe Test for Homogeneity of Variance using median
#>       Df F value Pr(>F)
#> group  7  0.0859 0.9984
#>       14
#> 
#> Testing Normality Assumption with Shapiro-Wilk
#>    *** Possible violation of the assumption.  You may 
#>             want to plot the residuals to see how they vary from normal ***
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  MyAOV_residuals
#> W = 0.83071, p-value = 0.001577
#> 
#> Bayesian analysis of models in order
#> # A tibble: 4 × 4
#>   model                            bf support                    margin_of_error
#>                                                             
#> 1 time                      82116.    " data support is extreme"     0.000000283
#> 2 time + Assay              32418.    " data support is extreme"     0.0270     
#> 3 time + Assay + time:Assay  9663.    " data support is extreme"     0.0129     
#> 4 Assay                         0.385 " data support is anecdot…     0.0000144
#> 
#> Interaction graph plotted...

^{Created on 2022-07-12 by the reprex package (v2.0.1)}

How to add a Scheffe Test to ggplot?

1 Answers