CRITICAL issue: p values depend on the order of factor levels: Games-Howell test

[GegznaV]

The p values of the Games-Howell test depend on the order of factor levels (find the p.adj and p.adj.signif columns in the examples below). I find this as a critical issue, which must be solved as soon as possible.

library(tidyverse)
library(rstatix)
#> 
#> Attaching package: 'rstatix'
#> The following object is masked from 'package:stats':
#> 
#>     filter

# Incorrect p values ------------------------------------------------
levels(PlantGrowth$group)
#> [1] "ctrl" "trt1" "trt2"

rstatix::games_howell_test(PlantGrowth, weight ~ group)
#> # A tibble: 3 x 8
#>   .y.    group1 group2 estimate conf.low conf.high p.adj p.adj.signif
#> * <chr>  <chr>  <chr>     <dbl>    <dbl>     <dbl> <dbl> <chr>       
#> 1 weight ctrl   trt1     -0.371   -1.17      0.430 1     ns          
#> 2 weight ctrl   trt2      0.494   -0.101     1.09  0.113 ns          
#> 3 weight trt1   trt2      0.865    0.114     1.62  0.024 *

# Incorrect p values ------------------------------------------------
PlantGrowth2 <- 
  PlantGrowth %>% 
  mutate(group = fct_reorder(group, weight, .desc = TRUE))  
levels(PlantGrowth2$group)
#> [1] "trt2" "ctrl" "trt1"

rstatix::games_howell_test(PlantGrowth2, weight ~ group)
#> # A tibble: 3 x 8
#>   .y.    group1 group2 estimate conf.low conf.high p.adj p.adj.signif
#> * <chr>  <chr>  <chr>     <dbl>    <dbl>     <dbl> <dbl> <chr>       
#> 1 weight trt2   ctrl     -0.494    -1.09     0.101     1 ns          
#> 2 weight trt2   trt1     -0.865    -1.62    -0.114     1 ns          
#> 3 weight ctrl   trt1     -0.371    -1.17     0.430     1 ns

# Correct p values --------------------------------------------------
PlantGrowth3 <- 
  PlantGrowth %>% 
  mutate(group = fct_reorder(group, weight))
levels(PlantGrowth3$group)
#> [1] "trt1" "ctrl" "trt2"

rstatix::games_howell_test(PlantGrowth3, weight ~ group)
#> # A tibble: 3 x 8
#>   .y.    group1 group2 estimate conf.low conf.high p.adj p.adj.signif
#> * <chr>  <chr>  <chr>     <dbl>    <dbl>     <dbl> <dbl> <chr>       
#> 1 weight trt1   ctrl      0.371   -0.430      1.17 0.475 ns          
#> 2 weight trt1   trt2      0.865    0.114      1.62 0.024 *           
#> 3 weight ctrl   trt2      0.494   -0.101      1.09 0.113 ns

^{Created on 2020-04-27 by the reprex package (v0.3.0)}

You may compare the results with PMCMRplus::gamesHowellTest() and userfriendlyscience::posthocTGH()

[kassambara]

Negative t-statistic value was not correctly handled when computing p-value using the function ptukey(). Fixed now, thank you for pointing this out!

[genghiskhanofnz]

Hi, i find it odd too.
I am analysing the total cost of providing heart failure for 3 different DRG groups: F62A (with complications), F62B (without complications) and F62C (mild).

I think the function games_howell_test is not working properly on my data.
Or, I am doing wrong thing.
Please advise.

My data violates all classic ANOVA assumptions, so welch_anova_test is done. p is significant

f62heart %>% welch_anova_test(TotalAmount ~ drg_factor)

A tibble: 1 x 7

.y. n statistic DFn DFd p method

• 
  1 TotalAmount 429 52.3 2 98.9 3.20e-16 Welch ANOVA

So move on to games_howell_test >> all p report ns.

f62heart %>% games_howell_test(TotalAmount ~ drg_factor, conf.level = 0.95, detailed = TRUE)

A tibble: 3 x 14

.y. group1 group2 n1 n2 estimate conf.low conf.high se statistic df p.adj p.adj.signif method

• 
  1 TotalAmount F62A F62B 171 230 -8522. -10939. -6105. 724. -8.33 202. 1 ns Games-Howell
  2 TotalAmount F62A F62C 171 28 -11291. -13897. -8685. 780. -10.2 188. 1 ns Games-Howell
  3 TotalAmount F62B F62C 230 28 -2769. -4191. -1347. 416. -4.71 48.1 1 ns Games-Howell

Then i decide to test again with other package = userfriendlyscience package >> now p all significant
The data intrinsically = their means are different

posthocTGH(y = f62heart$TotalAmount, x = f62heart$drg_factor)
n means variances
F62A 171 12842 164032960
F62B 230 4320 20372235
F62C 28 1551 7202876

        diff  ci.lo ci.hi    t  df    p

F62B-F62A -8522 -10939 -6105 8.3 202 <.01
F62C-F62A -11291 -13897 -8685 10.2 188 <.01
F62C-F62B -2769 -4191 -1347 4.7 48 <.01

posthocTGH(y = f62heart$TotalAmount, x = f62heart$drg_factor, method = 'tukey')
n means variances
F62A 171 12842 164032960
F62B 230 4320 20372235
F62C 28 1551 7202876

        diff  ci.lo ci.hi   t  df    p

F62B-F62A -8522 -10604 -6440 9.6 426 <.01
F62C-F62A -11291 -15495 -7087 6.3 426 <.01
F62C-F62B -2769 -6896 1358 1.6 426 .26

f62heart.txt

[kassambara]

This issue is fixed. Please install the latest dev version and try.

devtools::install_github("kassambara/rstatix")

[genghiskhanofnz]

Excellent. Downloaded from CRAN version 0.5.0 just now.
Returns nicely = all pairwise p value significant.
Thank you very much. You are the STAR.
Health data need this non-parametric since they are not following any normal distribution.
See here.

we do post hoc games howell now

to see any improvement

f62heart %>% games_howell_test(TotalAmount ~ drg_factor, conf.level = 0.95, detailed = TRUE)

A tibble: 3 x 14

.y. group1 group2 n1 n2 estimate conf.low conf.high se statistic df p.adj p.adj.signif method

• 
  1 TotalAm~ F62A F62B 171 230 -8522. -10939. -6105. 724. 8.33 202. 9.33e-14 **** Games-Ho~
  2 TotalAm~ F62A F62C 171 28 -11291. -13897. -8685. 780. 10.2 188. 0. **** Games-Ho~
  3 TotalAm~ F62B F62C 230 28 -2769. -4191. -1347. 416. 4.71 48.1 6.30e- 5 **** Games-Ho~