Paarweise Vergleiche aller Level einer Variable in GLM?

Varianzanalyse, Diskriminanzanalyse, Kontingenzanalyse, Faktorenanalyse, Clusteranalyse, MDS, ....

Moderator: EDi

Antworten
Romunkel

Paarweise Vergleiche aller Level einer Variable in GLM?

Beitrag von Romunkel »

Hallo liebes Forum,

vorab: ich bin kein Experte, nur einfacher R Anwender. Bitte seht mir nach wenn ich das vermeintlich Logische nicht immer gleich sehe.

Ich habe bei einem Journal ein Manuskript eingereicht, indem ich logistische (glm()) und negativ-binomiale Regressionsmodelle (glm.nb(), MASS package) verwende. Darin habe ich zwei unabhängige, kategoriale Variablen A und B mit 6 bzw. 4 Leveln, sowie eine binominale (logit) bzw. ganzzahlig diskrete (Zähl-) Variable (negbin) als abhängige Variable. Meine Stichprobengröße ist n=256 und es kommt erschwerend dazu, dass nicht für jede Level-Kombination aus den beiden Faktoren Daten existieren.
Die Modelle spucken mir ja immer die Estimates, SE und p-values für die übrigen 5 bzw. 3 Level im Vergleich zum Referenzlevel aus, wobei ich mir die 95% Wald Konfidenzintervalle berechne, indem ich die Estimates +/- die SE*1.96 rechne (dann entsprechen diese den p-Werten im Modell im Gegensatz zu confint()). Schließlich nutze ich exp(), um die Odds Ratios bzw. Rate ratios und 95%KI der verglichenen Level zu erhalten. Soweit so gut.

Der Reviewer des Journals möchte nun, dass ich die Ergebnisse für alle Level-Kombinationen innerhalb der beiden Variablen zeige, natürlich korrigiert für mehrfaches Testen, was also 15 paarweise Vergleiche für Variable A bzw. 6 für Variable B bedeutet. Mir ist klar, dass damit der alpha-Fehler massiv steigt aber ich versuchs dennoch. Er empfiehlt dazu das lsmeans package. Wenn ich das aber mit Holm-Korrektur rechne

Code: Alles auswählen

lsmeans(model,pairwise~A,adjust="holm")
, kriege ich teilweise um den Faktor 14 erhöhte p-Werte! Das liegt aber nicht an der Holm-Korrektur, denn wenn ich es ohne Korrektur rechne

Code: Alles auswählen

lsmeans(model,pairwise~A)
sind sie imer noch viel höher als im glm output.

Kann ich nicht die beiden GLM einfach mehrmals mit 5 (Variable A) bzw. 3 (Variable B) verschiedenen Referenzlevel rechnen und dann die p-Werte manuell korrigieren? Oder habe ich die lsmeans Funktion falsch verwendet? Oder etwas völlig falsch verstanden?

Viele Grüße und Danke,
Romunkel
Zuletzt geändert von Romunkel am Di Okt 31, 2017 12:41 pm, insgesamt 1-mal geändert.
Benutzeravatar
EDi
Beiträge: 1599
Registriert: Sa Okt 08, 2016 3:39 pm

Re: Paarweise Vergleiche aller Level einer Variable in GLM?

Beitrag von EDi »

Hast du ein reproduzierbares Beispiel für uns? z.b. simulierte Daten...

Ich vermute der Fehler liegt in der Interpretation von lsmeans, bzw. falscher annahme was gerechnet wird (default ist nicht "keine korrektur").

Auch nutze ich immer contrast() und summary() von lsmeans.
Bitte immer ein reproduzierbares Minimalbeispiel angeben. Meinungen gehören mir und geben nicht die meines Brötchengebers wieder.

Dieser Beitrag ist lizensiert unter einer CC BY 4.0 Lizenz
Bild.
Romunkel

Re: Paarweise Vergleiche aller Level einer Variable in GLM?

Beitrag von Romunkel »

Ok, ich hab hier mal meinen Datensatz:

Code: Alles auswählen

A<-c("A1","A1","A5","A1","A1","A1","A1","A5","A1","A5","A5","A5","A5","A6","A5","A5","A5","A1","A5","A1","A6","A4","A5","A4","A5","A5","A4","A4","A4","A2","A5","A5","A2","A1","A1","A1","A1","A5","A6","A6","A6","A5","A5","A5","A5","A5","A6","A6","A6","A1","A4","A4","A5","A5","A4","A4","A4","A4","A4","A4","A5","A5","A5","A5","A5","A5","A1","A5","A5","A1","A5","A1","A5","A1","A1","A5","A1","A5","A5","A5","A5","A5","A5","A1","A1","A5","A5","A5","A5","A6","A6","A6","A6","A6","A6","A6","A4","A4","A4","A4","A5","A5","A4","A5","A4","A4","A4","A4","A4","A4","A5","A4","A4","A4","A4","A4","A4","A4","A4","A4","A2","A4","A1","A1","A5","A1","A5","A1","A1","A1","A3","A1","A1","A5","A1","A1","A6","A4","A4","A4","A4","A4","A4","A4","A5","A4","A4","A4","A5","A4","A4","A5","A1","A1","A5","A2","A1","A5","A1","A2","A1","A1","A1","A1","A1","A1","A1","A6","A6","A6","A6","A6","A5","A2","A4","A4","A4","A5","A2","A2","A5","A5","A4","A5","A4","A5","A5","A4","A5","A5","A4","A4","A5","A1","A3","A1","A5","A1","A1","A4","A1","A1","A2","A5","A1","A1","A1","A1","A5","A5","A6","A6","A6","A6","A4","A4","A4","A4","A4","A4","A5","A4","A4","A4","A4","A4","A4","A4","A4","A4","A4","A4","A5","A2","A4","A4","A4","A1","A5","A1","A5","A4","A2","A1","A5","A5","A4","A5","A5","A5","A1","A1","A5","A5","A5","A1")
B<-c(„B2","B2","B2","B2","B2","B2","B2","B2","B2","B2","B2","B2","B2","B4","B2","B2","B2","B2","B2","B2","B4","B3","B3","B3","B3","B3","B3","B3","B3","B1","B1","B1","B1","B1","B1","B1","B2","B2","B4","B4","B4","B3","B3","B3","B1","B1","B4","B4","B4","B4","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B2","B2","B2","B2","B2","B2","B2","B2","B2","B2","B2","B4","B4","B4","B4","B4","B4","B4","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B2","B4","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B2","B4","B4","B4","B4","B4","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B2","B2","B4","B4","B4","B4","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B3","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1","B1")
C<-c(0,0,0,0,0,1,0,1,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,4,0,1,1,0,0,0,0,0,0,26,3,58,30,7,0,0,1,0,2,3,14,0,5,10,2,3,8,11,7,2,9,22,31,27,10,4,9,1,4,1,10,1,2,35,1,11,5,3,0,2,6,6,4,1,5,5,2,6,0,0,1,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,7,1,0,0,101,29,15,1,37,25,35,27,21,2,14,40,4,22,27,3,0,2,7,1,4,0,0,0,0,3,9,2,3,3,0,1,13,5,8,52,13,45,10,19,5,10,20,0,4,13,1,0,0,0,0,0,0,0,5,1,0,0,1,4,0,1,1,0,2,1,1,0,1,0,0,14,14,4,7,5,14,1,6,16,2,4,30,5,6,3,4,6,0,8,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,1,0,0,0,0,1,1,0,0,0,0,1,0,1,2,0,1,0,0,0,1,1,0,1,0,0,1,0,0,0)
mydata<-data.frame(A,B,C)
Jetzt das logit (modelL) und negbin Modell (modelN):

Code: Alles auswählen

modelL=glm(I(C>0)~A+B,data=mydata,family="binomial")
modelN=glm.nb(C~A+B, data=mydata[mydata$C>0,])
Im negativ-binomialen Modell nutze ich nur die Zähldaten, die nicht 0 sind. Damit ist n=139. Im logit Modell ist n=256.
Am Beispiel des negbin Modells ist A2 (p=0.001), A5 (p=0.006) und A6 (p=0.007) gegenüber A1 signifikant (Das NA in Variable B4 kommt daher, weil alle A6 nur in B4 beobachtet wurden und umgekehrt bis auf eine Ausnahme auch. Somit sind beide Level kollinear. Das musste ich leider hier dulden und diskutieren)

Wenn ich jetzt mithilfe von lsmeans alle 15 Level-Kombinationen der 6 Level von A Holm-korrigiert haben möchte

Code: Alles auswählen

lsmeans(modelN,pairwise~A,adjust="holm")
dann ist nur noch A2 (p=0.022) signifikant, während A5 (p=0.0819) und A6 (NA) rausfallen. Das NA bei A6 rührt wahrscheinlich daher, dass A1 nur ein einziges mal mit A6 im selben B-Level (B4) vorkamen. Die anderen p-Werte sind jetzt um den Faktor 11 bis 15 erhöht. Katastrophe!

Daher meine Frage, ob ich lsmeans richtig verwende bzw. ob man die Modelle nicht mehrmals mit verschiedenen Referenzleveln rechnen kann und dann den p-Wert separat korrigiert?

Dank dir vielmals!
Benutzeravatar
EDi
Beiträge: 1599
Registriert: Sa Okt 08, 2016 3:39 pm

Re: Paarweise Vergleiche aller Level einer Variable in GLM?

Beitrag von EDi »

dann ist nur noch A2 (p=0.022) signifikant, während A5 (p=0.0819) und A6 (NA) rausfallen. Das NA bei A6 rührt wahrscheinlich daher, dass A1 nur ein einziges mal mit A6 im selben B-Level (B4) vorkamen. Die anderen p-Werte sind jetzt um den Faktor 11 bis 15 erhöht. Katastrophe!
Erscheint mir so korrekt. Das ist der Preis den man für 15 tests berappen muss (bei 20 tests ist ja im Mittel einer falsch).

Hier wie ich es gemacht hätte:

Code: Alles auswählen

lsm <- lsmeans(modelN, ~A)
cc <- contrast(lsm, method = 'pairwise')
summary(cc, adjust = 'none', infer = TRUE) # no adjustment, same as tests from summary(modelN)
summary(cc, adjust = 'holm', infer = TRUE) # using holms procedure
Die NAs kommen daher, dass die Daten sehr unbalanced sind und man nicht alles schätzen kann (vorallem A6 und B4 machen Probleme).
Ohne korrektur kommt das gleiche bei rum wie bei den Wald tests für die Model-koeffizienten.
Daher meine Frage, ob ich lsmeans richtig verwende bzw. ob man die Modelle nicht mehrmals mit verschiedenen Referenzleveln rechnen kann und dann den p-Wert separat korrigiert?
Technisch machbar, aber wozu gibts denn lsmeans? Man muss natürlich alle tests berücksichtigen und würde zum gleichen Ergebniss kommen.

Was du machen könntest:
* Anstatt Holm eine andere (weniger konservative) Methode nutzen.
* Weniger Tests machen (z.b. A6 rauswerfen, dann sinds nur noch 10 Tests - vielleicht sind auch nicht alle vergleiche interessant, sondern nur zu A1; vielleicht gibts es a priori hypothesen)
* Nicht anpassen (kann man auch durchs Review bekommen...)
* Anstatt auf p-Wert zu schauen, Effektgrößen und Konfidenzintervalle betrachten (vieeeel informativer als p-Werte).
Bitte immer ein reproduzierbares Minimalbeispiel angeben. Meinungen gehören mir und geben nicht die meines Brötchengebers wieder.

Dieser Beitrag ist lizensiert unter einer CC BY 4.0 Lizenz
Bild.
Romunkel

Re: Paarweise Vergleiche aller Level einer Variable in GLM?

Beitrag von Romunkel »

Oh danke, das hilft mir sehr, ebenso wie die Tipps!
Mit deinem Befehl seh ich jetzt auch die unkorrigierten p-Werte für alle Levelvergleiche, das hatte ich nicht hingekriegt.

Jetzt wundert mich nur noch, warum das glm die Estimates, SE und p-Werte für A6 gegen A1 berechnen kann und die contrast()-Funktion nicht?
Ist das nur ein Manko der contrast()-Funktion? Oder macht es denn tatsächlich keinen Sinn diese Werte aus dem glm (z.B. für A6 vs. A1) zu verwenden? Sonst könnte ich die Modelle ja mit relevel() mit verschiedenen Referenzleveln berechnen bis ich alle Daten habe. Die 15 p-Werte müsste ich dann aber doch separat korrigieren (Die Holm-Methode hat hier sowieso scheinbar für 15 Tests korrigiert, obwohl nur 10 berechnet wurden (der Rest NA)).
Romunkel

Re: Paarweise Vergleiche aller Level einer Variable in GLM?

Beitrag von Romunkel »

Ich glaub ich kann mir meine Frage selbst beantworten. Wenn man nämlich die Level A6 mit A1 vergleicht kommen für verschiedene Level der Variable B als Referenz verschiedene Estimates, SE und p-Werte angezeigt, während diese bei Vergleichen innerhalb A1-A5 konstant bleibt, wie es ja sein sollte:

Code: Alles auswählen

summary(glm.nb(C~A+relevel(B,ref="B1"),data=mydata[mydata$C>0,]))
summary(glm.nb(C~A+relevel(B,ref="B2"),data=mydata[mydata$C>0,]))
Da stimmt was im glm output nicht und das lsmeans package liegt wohl richtig, wenn es NA anzeigt.
Benutzeravatar
EDi
Beiträge: 1599
Registriert: Sa Okt 08, 2016 3:39 pm

Re: Paarweise Vergleiche aller Level einer Variable in GLM?

Beitrag von EDi »

Da stimmt was im glm output nicht und das lsmeans package liegt wohl richtig, wenn es NA anzeigt.
lsmeans rechnet die Kontraste für A mit dem Mittel von B, das macht die summary von glm nicht - ich vermute (!) daher kommt der unterschied ?
Das grundlegende Problem ist aber die Datenbasis, die das nicht einfach macht.

Deine Methode funktioniert auch, für p-werte anzupassen kannst du dann p.adjust() verwenden.
Bitte immer ein reproduzierbares Minimalbeispiel angeben. Meinungen gehören mir und geben nicht die meines Brötchengebers wieder.

Dieser Beitrag ist lizensiert unter einer CC BY 4.0 Lizenz
Bild.
Antworten