2015-10-24 35 views
14

Ich passte ein Mutinomialmodell multinom Funktion mit (in diesem Fall auf Daten geben die Diät Präferenz von männlichen und weiblichen und unterschiedlicher Größe Klassen von Alligatoren in verschiedenen Seen):R: Tukey postthoc Tests für nnet multinom Multinomial fit für die Prüfung auf Gesamtunterschiede in Multinomialverteilung

data=read.csv("https://www.dropbox.com/s/y9elunsbv74p2h6/alligator.csv?dl=1") 
head(data) 
    id size sex lake food 
1 1 <2.3 male hancock fish 
2 2 <2.3 male hancock fish 
3 3 <2.3 male hancock fish 
4 4 <2.3 male hancock fish 
5 5 <2.3 male hancock fish 
6 6 <2.3 male hancock fish 

library(nnet) 
fit=multinom(food~lake+sex+size, data = data, Hess = TRUE) 

Die allgemeine Bedeutung meiner Faktoren, die ich bekommen kann

library(car) 
Anova(fit, type="III") # type III tests 
Analysis of Deviance Table (Type III tests) 

Response: food 
    LR Chisq Df Pr(>Chisq)  
lake 50.318 12 1.228e-06 *** 
sex  2.215 4 0.696321  
size 17.600 4 0.001477 ** 

Und Wirkungskurven unter Verwendung ich habe zB für den Faktor „See“

library(effects) 
plot(effect(fit,term="lake"),ylab="Food",type="probability",style="stacked",colors=rainbow(5)) 

enter image description here

Zusätzlich zu den Gesamt Anova-Tests unter Verwendung von mag ich auch auch paarweise Tukey posthoc Tests durchzuführen, obwohl für die allgemeinen Unterschiede in der Multinomialverteilung davon testen Beute Gegenstände werden gegessen, z über verschiedene Seenpaare.

Ich dachte zuerst an die Funktion glht in Paket multcomp, aber dies scheint nicht zu funktionieren, z.B. für Faktor lake:

library(multcomp) 
summary(glht(fit, mcp(lake = "Tukey"))) 
Error in summary(glht(fit, mcp(lake = "Tukey"))) : 
    error in evaluating the argument 'object' in selecting a method for function 'summary': Error in glht.matrix(model = list(n = c(6, 0, 5), nunits = 12L, nconn = c(0, : 
    ‘ncol(linfct)’ is not equal to ‘length(coef(model))’ 

Alternative war Paket lsmeans dafür zu verwenden, für die ich

versucht, obwohl
lsmeans(fit, pairwise ~ lake | food, adjust="tukey", mode = "prob") 
$contrasts 
food = bird: 
contrast    estimate   SE df t.ratio p.value 
george - hancock -0.04397388 0.05451515 24 -0.807 0.8507 
george - oklawaha 0.03680712 0.03849268 24 0.956 0.7751 
george - trafford -0.02123255 0.05159049 24 -0.412 0.9760 
hancock - oklawaha 0.08078100 0.04983303 24 1.621 0.3863 
hancock - trafford 0.02274133 0.06242724 24 0.364 0.9831 
oklawaha - trafford -0.05803967 0.04503128 24 -1.289 0.5786 

food = fish: 
contrast    estimate   SE df t.ratio p.value 
george - hancock -0.02311955 0.09310322 24 -0.248 0.9945 
george - oklawaha 0.19874095 0.09273047 24 2.143 0.1683 
george - trafford 0.32066789 0.08342262 24 3.844 0.0041 
hancock - oklawaha 0.22186050 0.09879102 24 2.246 0.1396 
hancock - trafford 0.34378744 0.09088119 24 3.783 0.0047 
oklawaha - trafford 0.12192695 0.08577365 24 1.421 0.4987 

food = invert: 
contrast    estimate   SE df t.ratio p.value 
george - hancock  0.23202865 0.06111726 24 3.796 0.0046 
george - oklawaha -0.13967425 0.08808698 24 -1.586 0.4053 
george - trafford -0.07193252 0.08346283 24 -0.862 0.8242 
hancock - oklawaha -0.37170290 0.07492749 24 -4.961 0.0003 
hancock - trafford -0.30396117 0.07129577 24 -4.263 0.0014 
oklawaha - trafford 0.06774173 0.09384594 24 0.722 0.8874 

food = other: 
contrast    estimate   SE df t.ratio p.value 
george - hancock -0.12522495 0.06811177 24 -1.839 0.2806 
george - oklawaha 0.03499241 0.05141930 24 0.681 0.9035 
george - trafford -0.08643898 0.06612383 24 -1.307 0.5674 
hancock - oklawaha 0.16021736 0.06759887 24 2.370 0.1103 
hancock - trafford 0.03878598 0.08135810 24 0.477 0.9635 
oklawaha - trafford -0.12143138 0.06402725 24 -1.897 0.2560 

food = rep: 
contrast    estimate   SE df t.ratio p.value 
george - hancock -0.03971026 0.03810819 24 -1.042 0.7269 
george - oklawaha -0.13086622 0.05735022 24 -2.282 0.1305 
george - trafford -0.14106384 0.06037257 24 -2.337 0.1177 
hancock - oklawaha -0.09115595 0.06462624 24 -1.411 0.5052 
hancock - trafford -0.10135358 0.06752424 24 -1.501 0.4525 
oklawaha - trafford -0.01019762 0.07161794 24 -0.142 0.9989 

Results are averaged over the levels of: sex, size 
P value adjustment: tukey method for comparing a family of 4 estimates 

Dies führt Tests für Unterschiede in dem relativen Anteil jeder spezifischen Art von Nahrungsmitteln. Ich frage mich, ob es auf die eine oder andere Art auch möglich wäre, Tukey-Posthoc-Tests zu erhalten, bei denen die gesamten multinomialen Verteilungen über die verschiedenen Seen hinweg verglichen werden, dh wo Unterschiede im Verhältnis zu den Beutetieren getestet werden Gegenstände gegessen? Ich habe versucht, mit

lsmeans(fit, pairwise ~ lake, adjust="tukey", mode = "prob") 

aber das scheint nicht zu funktionieren:

$contrasts 
contrast     estimate   SE df t.ratio p.value 
george - hancock  3.252607e-19 1.879395e-10 24  0 1.0000 
george - oklawaha -8.131516e-19 1.861245e-10 24  0 1.0000 
george - trafford -1.843144e-18 2.504062e-10 24  0 1.0000 
hancock - oklawaha -1.138412e-18   NaN 24  NaN  NaN 
hancock - trafford -2.168404e-18   NaN 24  NaN  NaN 
oklawaha - trafford -1.029992e-18   NaN 24  NaN  NaN 

Irgendwelche Gedanken?

Oder würde jemand wissen, wie glht für multinom Modelle gemacht werden könnte?

+1

Neben der Tatsache, dass Sie scheinbar die gleiche Frage zweimal gestellt haben, scheint es zu viele Informationen zu diesen Fragen zu geben. Können Sie nicht auf das spezifische Problem eingrenzen, anstatt alles, was Sie bisher getan haben, zu posten? – Molx

+0

Ich habe nicht zweimal die gleiche Frage gestellt - diese Frage bezieht sich auf ein multinomiales Modell, das nnets Multinom verwendet, und die andere Frage nach einem kumulativen Logit-Modell mit proportionalen Quoten unter Verwendung von MASS's polr - das sind zwei völlig verschiedene Dinge! Ich habe meine Frage ein bisschen gekürzt, obwohl ....(Ich dachte nur, dass die anderen Informationen nützlich gewesen sein könnten, um das Problem mit lsmeans zu diagnostizieren) –

+0

In Bezug auf den letzten Teil habe ich bereits versucht, Ihnen in einer persönlichen Mitteilung zu erklären, dass lsmeans von 'lake' Mittelung über' food' beinhalten. was ist Ihre multinomiale Antwort. Diese multinomialen Wahrscheinlichkeiten summieren sich bei allen anderen Faktoren notwendigerweise auf 1 über der Nahrung, wenn Sie also über die Nahrung gehen, erhalten Sie immer $ 0,2 $. Die flockig aussehenden Vergleiche sind also nur Rundungsfehler. Das ist keine Möglichkeit, die Seen zu vergleichen. Wenn Sie paarweise Vergleiche wollen, dann formulieren Sie eine Art lineare Kombination der Wahrscheinlichkeiten für jeden See und vergleichen Sie diese Werte. – rvl

Antwort

2

gerade eine Art Nachricht von Russ Lenth, und er dachte, die Syntax diese paarweise Vergleiche zwischen Seen zu tun, um Unterschiede in der Multinomialverteilung von zu testen, was Lebensmittel die Alligatoren

lsm = lsmeans(fit, ~ lake|food, mode = "latent") 
cmp = contrast(lsm, method="pairwise", ref=1) 
test = test(cmp, joint=TRUE, by="contrast") 
There are linearly dependent rows - df are reduced accordingly 
test 
contrast   df1 df2  F p.value 
george - hancock  4 24 3.430 0.0236 
george - oklawaha  4 24 2.128 0.1084 
george - trafford  4 24 3.319 0.0268 
hancock - oklawaha 4 24 5.820 0.0020 
hancock - trafford 4 24 5.084 0.0041 
oklawaha - trafford 4 24 1.484 0.2383 

Dank wäre essen Russ!

+0

sind diese Tukey korrigierten p-Werte? – WCMC

+0

Guter Punkt - ich denke nicht wirklich ... Aber Sie könnten immer noch eine Korrektur für mehrere Tests mit p.adjust in der multicomp-Paket anwenden ... –

+0

Versuchen Sie dies jetzt, mit 'emmeans' anstelle von' lsmeans'. Haben Sie schon mal eine Interaktion in einem 'nnet :: multinom()' Modell mit 'emmeans' getestet? Es wirft mir ein paar Fehler auf. –