0

my question about displaying the homogeneous subsets after a Tukey HSD test is already been asked in the past (see https://stats.stackexchange.com/questions/31547/how-to-obtain-the-results-of-a-tukey-hsd-post-hoc-test-in-a-table-showing-groupe, at the very end of the discussion) but never been answered. I'm running in the same problem now and I desperately need a solution as this function seems to be the only function available for R that displays such groups. If otherwise, please tell me. In short, this is what i get:

Groups, Treatments and means a 140095-001 36.79 b 150004-001 32.1 b 136936-021 31.97 bc 137219-004 31.39 bc 136673-017 31.27 bc 136963-009 30.79 bcd 147328-016 30.63 bcd 0147592-01 30.55 bcde 140094-001 30.02 cde 136730-007 29.7 de 136963-066 29.49 ef 136936-004 28.4 efg 147414-004 28.2 efg 137109-036 28.2 efg 136765-001 28.06 efg 140089-001 27.82 efg 137186-020 27.8 fg 136936-006 27.48 fgh 147350-014 27.43 gh 136992-001 27.36 gh 136730-015 27.18 ghi 0147785-01 27.08 ghi 0147691-01 26.98 ghi 136891-010 26.7 ghij 0147792-01 26.49 ghijk 136947-014 26.3 ghijkl 140097-001 25.8

And this is what I want:

    Groups, Treatments and means
a        2.1     51.17547 
ab       4.1     50.7529 
abc      3.1     47.36229 
 bcd     1.1     45.81229 
  cd     5.1     44.55313 
   de    4.0     41.81757 
    ef   2.0     38.79482 
    ef   1.0     36.91257 
     f   3.0     36.34383 
     f   5.0     35.69507

I already checked the script behind the function but since my knowledge is rather low i couldn't find anything, at least nothing straight and clear. Could anybody help me with a solution? I would be very pleased with an answer.

Thanks in advance, Jelle

Community
  • 1
  • 1
Jelle
  • 1
  • 2

1 Answers1

0

First build e.g. anova model by aov function

my_model <- stats::aov(data$value ~ data$treatment) 
summary(my_model)

Post hoc test with Tukey HSD, use anova model as input data, column treatment in original dataset is grouping factor.

test_result <- agricolae::HSD.test(my_model, "my_model$treatment", group = TRUE, console = TRUE)

And check the results in the console. There is, among others, also table with your asked overview. (not significant in my example)

  Groups, Treatments and means
 a   2008    6.755 
 a   2013    6.147 
 a   1975    4.653 
 a   1997    3.622

It is possible to visualize results for treatments. 

agricolae::bar.group(test_result$groups, ylim = c(0,10), density = 10, border = "blue")

Visualization for Tukey HSD

tlask
  • 198
  • 9