lsmeans (R): Passen Sie mehrere Vergleiche mit Interaktionstermen an

7

Ich habe ein lsmeans-Problem in R. Ich möchte eine Post-hoc-Analyse einer Interaktion durchführen, ähnlich den Beispielen in der lsmeans-Dokumentation.

Ich bin verwirrt darüber, dass die p-Werte gleich sind, ob ich sie benutze

warp.lm <- lm(breaks ~ wool * tension, data = warpbreaks)

(EIN):

lsmeans(warp.lm, list(pairwise ~ wool|tension, pairwise ~ tension|wool))

oder (B):

lsmeans(warp.lm, list(pairwise ~ wool|tension))

Sollten die p-Werte bei (A) nicht höher sein als bei (B)? (A) führt 9 Post-hoc-Vergleiche durch, (B) nur 3.

Kann ich sicher sein, dass (A) für alle 9 Vergleiche angepasst ist?

Vielen Dank!

Sam
quelle
Was meinst du richtig? Ich schlage vor, die lsmeans-Funktion zuerst ohne die linke Seite zu verwenden und dann zu verwenden pairs, um die gewünschten Kontraste zu erhalten. Lesen Sie die Vignette "using-lsmeans" für einige Beispiele und Diskussionen.
Russ Lenth
@rvl Sorry das war nicht klar. Für beide Befehle erhalten Sie die gleichen p-Werte. Spannung = L: Kontrastschätzung SE df t.Verhältnis p.Wert A - B 16.333333 5.157299 48 3.167 0.0027 Spannung = M: A - B -4.777778 5.157299 48 -0.926 0.3589 Spannung = H: A - B 5.777778 5.157299 48 1.120 0.2682 Sollte die p-Werte nach dem ersten Befehl nicht höher sein als nach dem zweiten? Nach dem ersten Befehl müssen weitere Vergleiche durchgeführt werden. Ich habe die "Paare" -Route ausprobiert, aber die Ergebnisse sind die gleichen wie "paarweise".
Sam
Die p-Wert-Anpassungen werden für jede Gruppierung separat angewendet. Wenn Sie eine Anpassung für den gesamten Satz von Vergleichen wünschen, müssen Sie benutzerdefinierte Kontraste vornehmen.
Russ Lenth
Ich habe einen Antwortvorschlag zum Posten, kann ihn aber nicht posten, da die Frage zurückgestellt wurde. Wenn es aus der Warteschleife genommen werden kann, werde ich es posten.
Russ Lenth
@Scortchi Könnten Sie bitte die Frage aus der Warteschleife nehmen, damit rvl sie beantworten kann? Ich habe es so bearbeitet, dass Sie es hoffentlich klarer finden. Vielen Dank.
Sam

Antworten:

7

Beide lsmeansAnweisungen, die Sie anzeigen, generieren Listen von lsmobjs, und jedes Element dieser Listen wird separat behandelt. Wenn Sie eine Gesamtanpassung für zwei oder mehr Listen zusammen vornehmen möchten, ist dies technisch und erfordert ein wenig Arbeit.

Speichern Sie zunächst die Liste:

lsmlist = lsmeans(warp.lm, list(pairwise ~ wool|tension, 
                  pairwise ~ tension|wool))

Dies erstellt eine Liste von 4 lsmobjs (ursprünglich zwei Listen von zwei)

> names(lsmlist)
[1] "lsmeans of wool | tension"                          
[2] "pairwise differences of contrast, tension | tension"
[3] "lsmeans of tension | wool"                          
[4] "pairwise differences of contrast, wool | wool"

Der Benutzer möchte die drei Vergleiche lsmlist[[2]]mit den sechs in kombinieren lsmlist[[4]]und eine Gesamtmultiplizitätsanpassung für diese 9 Vergleiche vornehmen.

Erstellen Sie zunächst lsmobjaus einem der Ergebnisse ein neues und korrigieren Sie es.

mydiffs = lsmlist[[4]]

Verbinden Sie zunächst die linearen Funktionen für die beiden Vergleichssätze:

mydiffs@linfct = rbind(lsmlist[[2]]@linfct, lsmlist[[4]]@linfct)

Wir müssen auch den gridSchlitz definieren, der die mit jeder linearen Funktion verbundenen Faktoren definiert. Um es einfach zu machen, definiere ich einfach einen Faktor contrastmit 9 Ebenen für die 9 Kontraste (hier könnte etwas ausgefalleneres getan werden).

mydiffs@grid = data.frame(contrast = 1:9)

Korrigieren Sie abschließend die Zusatzinformationen für die Buchhaltung. Wir verwenden jetzt unseren neuen contrastFaktor als einzige Variable ohne "by" -Variablen. Für die kombinierte Familie von 9 Kontrasten macht die Tukey-Anpassung keinen Sinn, daher verwenden wir die multivariatet( "mvt") Methode:

mydiffs = update(mydiffs, pri.vars = "contrast", by.vars = NULL,
                 adjust="mvt")

Nun können wir uns die resultierende Zusammenfassung ansehen:

> mydiffs

 contrast   estimate       SE df t.ratio p.value
        1 16.3333333 5.157299 48   3.167  0.0205
        2 -4.7777778 5.157299 48  -0.926  0.9119
        3  5.7777778 5.157299 48   1.120  0.8258
        4 20.5555556 5.157299 48   3.986  0.0019
        5 20.0000000 5.157299 48   3.878  0.0026
        6 -0.5555556 5.157299 48  -0.108  1.0000
        7 -0.5555556 5.157299 48  -0.108  1.0000
        8  9.4444444 5.157299 48   1.831  0.3796
        9 10.0000000 5.157299 48   1.939  0.3190

P value adjustment: mvt method for 9 tests

Ich könnte in Betracht ziehen, eine Funktion zum Kombinieren von lsm.listObjektteilen hinzuzufügen , aber dies ist nicht einfach - hauptsächlich aufgrund von Komplikationen beim Erhalten aussagekräftiger Beschriftungen für die Ergebnisse (das gridTeil). Es ist auch ein Problem, dass verschiedene Benutzer unterschiedliche Standardeinstellungen erwarten.

Russ Lenth
quelle
fantastische Antworten wie immer.
Adam Robinsson
Vielen Dank für Ihre Antwort, das ist genial! Wenn es eine Warnmeldung geben könnte, die nützlich wäre, denke ich. Diese Standardeinstellung führt zu vielen Fehlalarmen.
Sam
Bitte beachten Sie - ich habe eine Frage zu diesem Thema gestellt
Russ Lenth
8

Bereitstellung in der kommenden Version von lsmeans

Das nächste Update von lsmeans (2.20 oder höher) enthält eine rbindMethode für ref.gridund lsmobjObjekte. Es macht es einfach, zwei oder mehr Objekte zu einer Familie zu kombinieren, und standardmäßig wird die "mvt"Anpassungsmethode verwendet. Hier ist das vorliegende Beispiel:

> w.t <- pairs(lsmeans(warp.lm, ~ wool | tension))
> t.w <- pairs(lsmeans(warp.lm, ~ tension | wool))

> rbind(w.t, t.w)
 tension wool contrast   estimate       SE df t.ratio p.value
 L       .    A - B    16.3333333 5.157299 48   3.167  0.0203
 M       .    A - B    -4.7777778 5.157299 48  -0.926  0.9119
 H       .    A - B     5.7777778 5.157299 48   1.120  0.8258
 .       A    L - M    20.5555556 5.157299 48   3.986  0.0018
 .       A    L - H    20.0000000 5.157299 48   3.878  0.0027
 .       A    M - H    -0.5555556 5.157299 48  -0.108  1.0000
 .       B    L - M    -0.5555556 5.157299 48  -0.108  1.0000
 .       B    L - H     9.4444444 5.157299 48   1.831  0.3793
 .       B    M - H    10.0000000 5.157299 48   1.939  0.3191

P value adjustment: mvt method for 9 tests 

Wenn Sie Tukey anstelle von mvt angeben, kommt dies dem gleichen Ergebnis ziemlich nahe:

> test(rbind(w.t, t.w), adjust = "tukey")
 tension wool contrast   estimate       SE df t.ratio p.value
 L       .    A - B    16.3333333 5.157299 48   3.167  0.0196
 M       .    A - B    -4.7777778 5.157299 48  -0.926  0.8683
 H       .    A - B     5.7777778 5.157299 48   1.120  0.7727
 .       A    L - M    20.5555556 5.157299 48   3.986  0.0018
 .       A    L - H    20.0000000 5.157299 48   3.878  0.0026
 .       A    M - H    -0.5555556 5.157299 48  -0.108  0.9999
 .       B    L - M    -0.5555556 5.157299 48  -0.108  0.9999
 .       B    L - H     9.4444444 5.157299 48   1.831  0.3467
 .       B    M - H    10.0000000 5.157299 48   1.939  0.2922

P value adjustment: tukey method for comparing a family of 4.772 estimates

Die merkwürdig gebrochene Familiengröße ergibt sich aus .(4.7722)=9

Russ Lenth
quelle