Uzyskiwanie wartości p dla „multinom” w R (pakiet nnet)

Jak uzyskać wartości p za pomocą multinomfunkcji nnetpakiet w R?

Mam zestaw danych, który składa się z „wyników patologii” (nieobecny, łagodny, ciężki) jako zmiennej wynikowej oraz dwóch głównych efektów: wieku (dwa czynniki: dwadzieścia / trzydzieści dni) i grupy leczenia (cztery czynniki: zainfekowany bez ATB; zainfekowany + ATB1; zainfekowany + ATB2; zainfekowany + ATB3).

Najpierw próbowałem dopasować model regresji porządkowej, który wydaje się bardziej odpowiedni, biorąc pod uwagę cechy mojej zmiennej zależnej (porządkowej). Jednak założenie o proporcjonalności szans zostało poważnie naruszone (graficznie), co skłoniło mnie do użycia modelu wielomianowego zamiast nnetpakietu.

Najpierw wybrałem poziom wyniku, który muszę wykorzystać jako kategorię wyjściową:

Data$Path <- relevel(Data$Path, ref = "Absent")

Następnie musiałem ustawić podstawowe kategorie dla zmiennych niezależnych:

Data$Age <- relevel(Data$Age, ref = "Twenty")
Data$Treat <- relevel(Data$Treat, ref="infected without ATB") 

Model:

test <- multinom(Path ~ Treat + Age, data = Data) 
# weights:  18 (10 variable) 
initial value 128.537638 
iter 10 value 80.623608 
final  value 80.619911 
converged

Wyjście:

Coefficients:
         (Intercept)   infected+ATB1   infected+ATB2   infected+ATB3    AgeThirty
Moderate   -2.238106   -1.1738540      -1.709608       -1.599301        2.684677
Severe     -1.544361   -0.8696531      -2.991307       -1.506709        1.810771

Std. Errors:
         (Intercept)    infected+ATB1   infected+ATB2   infected+ATB3    AgeThirty
Moderate   0.7880046    0.8430368       0.7731359       0.7718480        0.8150993
Severe     0.6110903    0.7574311       1.1486203       0.7504781        0.6607360

Residual Deviance: 161.2398
AIC: 181.2398

Przez jakiś czas nie mogłem znaleźć sposobu na uzyskanie wartości dla modelu i oszacowań przy użyciu . Wczoraj natknąłem się na post, w którym autor przedstawił podobny problem dotyczący szacowania wartości dla współczynników ( Jak skonfigurować i oszacować wielomianowy model logit w R? ). Tam jeden bloger zasugerował, że uzyskanie wartości z wyniku jest dość łatwe, najpierw uzyskując wartości w następujący sposób:p p $p$nnet:multinom$p$$p$summarymultinom$t$

pt(abs(summary1$coefficients / summary1$standard.errors), df=nrow(Data)-10, lower=FALSE) 

         (Intercept)   infected+ATB1   infected+ATB2   infected+ATB3    AgeThirty
Moderate 0.002670340   0.08325396      0.014506395     0.02025858       0.0006587898
Severe   0.006433581   0.12665278      0.005216581     0.02352202       0.0035612114

Według Petera Dalgard „Jest co najmniej współczynnik 2 brakuje na dwustronny -value Zwykle jest błędem użycie. -Dystrybucja za to, co jest naprawdę -statistic; dla danych zagregowanych, może być bardzo zły błąd ”. Według Briana Ripleya „błędem jest również stosowanie testów Walda w przypadku napadów, ponieważ cierpią one na te same (potencjalnie poważne) problemy, co napady dwumianowe. Należy stosować przedziały ufności prawdopodobieństwa profilu (dla których pakiet zapewnia oprogramowanie), lub jeśli musisz przetestować, testy współczynnika wiarygodności (to samo). ”t $p$$t$$z$multinom

Muszę tylko móc uzyskać wiarygodne wartości .$p$

Co z używaniem

z <- summary(test)$coefficients/summary(test)$standard.errors
# 2-tailed Wald z tests to test significance of coefficients
p <- (1 - pnorm(abs(z), 0, 1)) * 2
p

Zasadniczo byłoby to oparte na oszacowanych współczynnikach w stosunku do ich błędu standardowego i wykorzystałoby test az do przetestowania znaczącej różnicy z zerą w oparciu o test dwustronny. Współczynnik dwóch koryguje problem, o którym wspominał Peter Dalgaard (potrzebujesz go, ponieważ chcesz testu dwustronnego, a nie jednostronnego), i używa testu Z, a nie testu T, aby rozwiązać drugi problem, o którym wspominasz.

Możesz również uzyskać ten sam wynik (testy Wald z) za pomocą

library(AER)
coeftest(test)

Testy ilorazu wiarygodności są jednak ogólnie uważane za dokładniejsze niż testy Walda z (te ostatnie wykorzystują normalne przybliżenie, testy LR nie), i można je uzyskać za pomocą

library(afex)
set_sum_contrasts() # use sum coding, necessary to make type III LR tests valid
library(car)
Anova(test,type="III")

Jeśli chcesz przeprowadzić testy posthoc Tukey'a w parach, możesz je uzyskać za pomocą lsmeanspakietu, jak wyjaśniono w moim innym poście !