Получение p-значений для «multinom» в R (пакет nnet)

Как я могу получить p-значения, используя multinomфункцию nnetпакета в R?

У меня есть набор данных, который состоит из «Оценка патологии» (Отсутствует, Легкая, Тяжелая) в качестве переменной результата и двух основных эффектов: Возраст (два фактора: двадцать / тридцать дней) и Группа лечения (четыре фактора: инфицированные без АТБ; инфицированные + ATB1; зараженный + ATB2; зараженный + ATB3).

Сначала я попытался приспособить модель порядковой регрессии, которая кажется более подходящей, учитывая характеристики моей зависимой переменной (порядковый номер). Однако предположение о пропорциональности шансов было серьезно нарушено (графически), что побудило меня вместо этого использовать полиномиальную модель с использованием nnetпакета.

Сначала я выбрал уровень результата, который мне нужно использовать в качестве базовой категории:

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

Затем мне нужно было установить базовые категории для независимых переменных:

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

Модель:

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

Выход:

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

Некоторое время я не мог найти способ получить для модели и оценки при использовании . Вчера я натолкнулся на пост, где автор выдвинул аналогичную проблему, касающуюся оценки значений для коэффициентов ( Как установить и оценить модель многочленного логита в R? ). Там один блоггер предположил, что получить из результата довольно просто, сначала получив значения следующим образом:р р $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

По словам Питера Далгарда, «для двухстороннего значения отсутствует, по крайней мере, коэффициент 2. Обычно ошибочно использовать -распределение для действительно статистики; для агрегированных данных это может быть очень плохая ошибка ". По словам Брайана Рипли, «также ошибочно использовать тесты Вальда для подгонок, поскольку они страдают от тех же (потенциально серьезных) проблем, что и при биномиальных подгонках. Используйте доверительные интервалы вероятности профиля (для которых пакет предоставляет программное обеспечение), или если вы должны проверить, тесты отношения правдоподобия (то же самое). "т $p$$t$$z$multinom

Мне просто нужно иметь возможность получить надежные .$p$

Как насчет использования

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

По сути, это будет основано на оценочных коэффициентах относительно их стандартной ошибки и будет использовать z-тест для проверки на значительную разницу с нулем на основе двустороннего теста. Коэффициент два исправляет проблему Питера Далгаарда, о которой говорилось выше (она нужна вам, потому что вам нужен двухсторонний тест, а не односторонний), и для решения другого используется z-тест, а не t-тест. проблема, которую вы упоминаете.

Вы также можете получить тот же результат (z-тесты Вальда), используя

library(AER)
coeftest(test)

Тесты отношения правдоподобия обычно считаются более точными, чем тесты Вальда (последние используют нормальное приближение, тесты LR - нет), и их можно получить, используя

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

Если вы хотите провести парные тесты Tukey posthoc, то их можно получить с помощью lsmeansпакета, как описано в моем другом посте !