0
$\begingroup$

I have a dataset that I have been trying to run Tukey HSD tests on and obtain a compact letter display. I originally used the agricolae package with function HSD.Test(). After reviewing my results, I was curious about comparisons between groups that seemed different from what I expected- not that I refuse to accept these results but because it honestly surprised me after I had reviewed raw data. I went and found the multcomp package to check if it would provide the same results from the first package with the functions glht() and cld(). Below is the code that shows the difference in compact letter displays, which I am ultimately wanting to use in a table. My question is: should I trust one over the other? Or is there a fundamental difference in how these are calculated between packages that would aid me in how to reconcile the different results I am getting.

Here is my setup and run of the agricolae.

#import libraries
library (agricolae)
library(multcomp)

#fit analysis of variance model
SPLW<-aov(sepalLengthToWidth~TAXON, data=garden)

#get compact letter comparisons from agricolae
hsd<-HSD.test(SPLW, "TAXON")
hsd

$statistics
   MSerror Df     Mean       CV
  0.219363 63 3.019053 15.51354

$parameters
   test name.t ntr StudentizedRange alpha
  Tukey  TAXON   6         4.156707  0.05

$means
     sepalLengthToWidth       std  r      Min      Max      Q25      Q50      Q75
AGCP           3.003620 0.4542374 24 2.185729 4.048658 2.781238 3.090727 3.254949
CPTN           2.516216 0.1663721  3 2.369989 2.697232 2.425709 2.481429 2.589330
NOAL           3.517557 0.5924274 13 2.439724 4.167752 3.237401 3.745421 3.983766
OUMT           2.633716 0.3395092  6 2.128205 3.056782 2.449174 2.646832 2.866976
TMAL           3.837945 0.7807588  3 3.290098 4.731935 3.390950 3.491803 4.111869
WGCP           2.781738 0.3981362 20 1.854508 3.302083 2.608329 2.882435 3.011886

$comparison
NULL

$groups
     sepalLengthToWidth groups
TMAL           3.837945      a
NOAL           3.517557      a
AGCP           3.003620     ab
WGCP           2.781738      b
OUMT           2.633716      b
CPTN           2.516216      b

attr(,"class")
[1] "group"

Then I use the multcomp method on the same aov model.

#get compact letter comparisons from multcomp
tuk <- glht(SPLW, linfct = mcp(TAXON = "Tukey"))
tuk.cld <- cld(tuk)   
tuk.cld

AGCP CPTN NOAL OUMT TMAL WGCP 
 "a"  "a"  "b"  "a"  "b"  "a" 

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.