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 

The output:

         (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.

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

What about using

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

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


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

set_sum_contrasts() # use sum coding, necessary to make type III LR tests valid

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:


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]
   ## 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.