Wie erhält man die Ergebnisse eines Tukey-HSD-Post-Hoc-Tests in einer Tabelle mit gruppierten Paaren?

13

Ich würde gerne einen TukeyHSD-Post-Hoc-Test nach meiner Zwei-Wege-Anova mit R durchführen und eine Tabelle mit den sortierten Paaren erhalten, die nach signifikanten Unterschieden gruppiert sind. (Entschuldigung für die Formulierung, ich bin noch neu mit Statistiken.)

Ich hätte gerne so etwas:

Bildbeschreibung hier eingeben

Also gruppiert mit Sternen oder Buchstaben.

Irgendeine Idee? Ich habe die Funktion HSD.test()aus dem agricolaePaket getestet , aber es scheint, dass sie keine bidirektionalen Tabellen verarbeitet.

stragu
quelle

Antworten:

22

Die agricolae::HSD.testFunktion macht genau das, aber Sie müssen es wissen lassen, dass Sie an einem Interaktionsbegriff interessiert sind . Hier ist ein Beispiel mit einem Stata-Datensatz:

library(foreign)
yield <- read.dta("http://www.stata-press.com/data/r12/yield.dta")
tx <- with(yield, interaction(fertilizer, irrigation))
amod <- aov(yield ~ tx, data=yield)
library(agricolae)
HSD.test(amod, "tx", group=TRUE)

Dies ergibt die unten gezeigten Ergebnisse:

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 

Sie stimmen mit dem überein, was wir mit den folgenden Befehlen erhalten würden:

. webuse yield
. regress yield fertilizer##irrigation
. pwcompare fertilizer#irrigation, group mcompare(tukey)

-------------------------------------------------------
                      |                           Tukey
                      |     Margin   Std. Err.   Groups
----------------------+--------------------------------
fertilizer#irrigation |
                 1 0  |   36.91257   1.116571    AB    
                 1 1  |   45.81229   1.116571      CDE 
                 2 0  |   38.79482   1.116571    AB    
                 2 1  |   51.17547   1.116571         F
                 3 0  |   36.34383   1.116571    A     
                 3 1  |   47.36229   1.116571       DEF
                 4 0  |   41.81757   1.116571     BC   
                 4 1  |    50.7529   1.116571        EF
                 5 0  |   35.69507   1.116571    A     
                 5 1  |   44.55313   1.116571      CD  
-------------------------------------------------------
Note: Margins sharing a letter in the group label are
      not significantly different at the 5% level.

Das Multcomp- Paket bietet auch eine symbolische Visualisierung ('kompakte Buchstabenanzeigen', siehe Algorithmen für kompakte Buchstabenanzeigen: Vergleich und Auswertung für weitere Einzelheiten) von signifikanten paarweisen Vergleichen, obwohl sie nicht in tabellarischer Form dargestellt werden. Es verfügt jedoch über eine Plotmethode, mit der die Ergebnisse mithilfe von Boxplots bequem angezeigt werden können. Die Präsentationsreihenfolge kann ebenfalls geändert werden (Option decreasing=), und es stehen viel mehr Optionen für mehrere Vergleiche zur Verfügung. Es gibt auch das multcompView- Paket, das diese Funktionen erweitert.

Hier ist das gleiche Beispiel analysiert mit glht:

library(multcomp)
tuk <- glht(amod, linfct = mcp(tx = "Tukey"))
summary(tuk)          # standard display
tuk.cld <- cld(tuk)   # letter-based display
opar <- par(mai=c(1,1,1.5,1))
plot(tuk.cld)
par(opar)

Behandlungen, die denselben Buchstaben haben, unterscheiden sich auf der gewählten Ebene nicht signifikant (Standard: 5%).

Bildbeschreibung hier eingeben

Übrigens gibt es ein neues Projekt, das derzeit auf R-Forge gehostet wird und vielversprechend aussieht: Factorplot . Es enthält zeilen- und buchstabenbasierte Darstellungen sowie eine Matrixübersicht (über ein Ebenendiagramm) aller paarweisen Vergleiche. Ein Arbeitspapier finden Sie hier: Factorplot: Verbesserte Darstellung einfacher Kontraste in GLMs

chl
quelle
Vielen Dank für diese ausführliche Antwort! Ich werde diese verschiedenen Methoden ausprobieren, sobald ich ein paar Minuten Zeit habe. Prost!
Stragu
Ich habe die multcomp package-Funktion ausprobiert. Wenn ich die Funktion 'cld ()' verwende, erhalte ich die Fehlermeldung 'Fehler: sapply (split_names, length) == 2 ist nicht alles WAHR'. Irgendeine Idee warum?
Stragu
1
@chtfn Es scheint ein Problem mit Variablenlabels zu geben. Ein kurzer Blick auf den Quellcode zeigt an, dass diese Fehlermeldung stammt, von insert_absorb()der versucht wird, zwei Behandlungen zu extrahieren. Sie können vielleicht versuchen, das Trennzeichen zu ändern, das Sie für die Codierungsebenen Ihres Interaktionsbegriffs verwendet haben. Ohne ein funktionierendes Beispiel ist es schwer zu sagen, was passiert ist.
chl
Ich habe es herausgefunden: Ich hatte Punkte in den Namen meiner Genotypen und Behandlungen und da qlht () einen Punkt verwendet, um die Paarnamen zu teilen, ist es ausgeflippt. Vielen Dank für all Ihre Hilfe, chl! :)
stragu
3
Ich bemerkte , dass ich heute nun hinzufügen muß console=TRUEin HSD.test()um die Tabellen zu erhalten, falls jemand versucht , dies und sieht kein Ergebnis. Wahrscheinlich ein Update von agricolae.
stragu
2

Es gibt eine Funktion mit dem Namen TukeyHSD, die gemäß der Hilfedatei eine Reihe von Konfidenzintervallen für die Unterschiede zwischen den Mitteln der Niveaus eines Faktors mit der angegebenen familienbezogenen Wahrscheinlichkeit der Abdeckung berechnet. Die Intervalle basieren auf der studentisierten Bereichsstatistik, Tukeys "Honest Significant Difference" -Methode. Macht das was du willst?

http://stat.ethz.ch/R-manual/R-patched/library/stats/html/TukeyHSD.html

smillig
quelle
Danke für Ihre Antwort. Ja, ich habe diese Funktion ausprobiert, aber sie gibt mir unformatierte Vergleichslisten. Ich möchte, dass sie wie im Bild in meiner Frage gruppiert sind, eine klare Sicht darauf haben, welche Gruppe sich von welcher unterscheidet, und schließlich die Gruppennamen in meine Diagramme einfügen (zum Beispiel: a, ab, abc, bc , c)
stragu