Getting p-values for “multinom” in R (nnet package)

How do I get p-values using the multinom function of nnet package in R?

I have a dataset which consists of “Pathology scores” (Absent, Mild, Severe) as outcome variable, and two main effects: Age (two factors: twenty / thirty days) and Treatment Group (four factors: infected without ATB; infected + ATB1; infected + ATB2; infected + ATB3).

First I tried to fit an ordinal regression model, which seems more appropriate given the characteristics of my dependent variable (ordinal). However, the assumption of odds proportionality was severely violated (graphically), which prompted me to use a multinomial model instead, using the nnet package.

First I chose the outcome level that I need to use as baseline category:

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


Then, I needed to set baseline categories for the independent variables:

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


The 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


The output:

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


For a while, I could not find a way to get the $p$-values for the model and estimates when using nnet:multinom. Yesterday I came across a post where the author put forward a similar issue regarding estimation of $p$-values for coefficients (How to set up and estimate a multinomial logit model in R?). There, one blogger suggested that getting $p$-values from the summary result of multinom is pretty easy, by first getting the $t$values as follows:

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


According to Peter Dalgard, "There's at least a factor of 2 missing for a two-tailed $p$-value. It is usually a mistake to use the $t$-distribution for what is really a $z$-statistic; for aggregated data, it can be a very bad mistake." According to Brian Ripley, "it is also a mistake to use Wald tests for multinom fits, since they suffer from the same (potentially severe) problems as binomial fits. Use profile-likelihood confidence intervals (for which the package does provide software), or if you must test, likelihood-ratio tests (ditto)."

I just need to be able to derive reliable $p$-values.

• You can use model comparisons with likelihood ratio tests for a full and reduced model using nnet's anova() function. – caracal Jul 3 '13 at 14:16

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


Basically, this would be based on the estimated coefficients relative to their standard error, and would use a z test to test against a significant difference with zero based on a two-tailed test. The factor of two corrects the problem Peter Dalgaard referred to above (you need it because you want a two tailed test, not a one tailed one), and it uses a z-test, rather than a t-test, to solve the other problem that you mention.

You can also get the same result (Wald z-tests) using

library(AER)
coeftest(test)


Likelihood ratio tests are generally regarded as more accurate though than Wald z tests (the latter use a normal approximation, LR tests do not), and these can be gotten using

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


If you would like to carry out pairwise Tukey posthoc tests, then these can be obtained using the lsmeans package as explained in my other post!

Also you could be interested in Likehood Ratio test p-values, as seen here:

http://thestatsgeek.com/2014/02/08/wald-vs-likelihood-ratio-test/

Wich you could extract like this (sorry, its a custom function :D)

likehoodmultinom_p <- function(model_lmm)
{

i <- 1

variables <-c("No funciona")
valores <- c("No funciona")

for (var in model_lmm$coefnames[-1]) { # Qutiamos el -1 de coefnames para no obener un NA variables[i] =paste(var) valores[i]= lrtest(model_lmm, var)[[5]][2] i=i+1 ## Contributed to stack at: https://stackoverflow.com/questions/23018238/assesing-the-goodness-of-fit-for-the-multinomial-logit-in-r-with-the-nnet-packag/60835647#60835647 } return (data.frame(variables,valores)) } L_iris= likehoodmultinom_p(iris_fit)  In my function you obtain a df with factors, so you maybe gotta change them a bit to extrac them. I have yet to correct my original function: L_iris= likehoodmultinom_p(iris_fit)  L_iris$valores = as.character(L_iris$valores) # Pass them as chr L_iris$valores = as.numeric(L_iris\$valores) # And as numeric.


Then you can acces them easily. I also usually sort them in base of p-values.