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

Question

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.

Answer 1

You can use that custom function that I`ve created, for example you can just get true or false for Hipotesis test.

The TRUE ones should represent P-Value > 0,5.

zWald_test <- function(x){
  zWald_modelo<- (summary(x)$coefficients / 
                                summary(x)$standard.errors)
  a <- t(apply(zWald_modelo, 1, function(x) {x < qnorm(0.025, lower.tail = FALSE)} ))
  b <- t(apply(zWald_modelo, 1, function(x) {x > -qnorm(0.025, lower.tail = FALSE)} ))
  ifelse(a==TRUE & b==TRUE, TRUE, FALSE)
}

or if you want to confirm the p-values you can use the following function:

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

Answer 2

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
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!

Answer 3

As already said by OP (by his quote of B Ripley), wald tests is not really good for multinomial models, we really should use likelihoodratio tests. Below I show an easy way of getting that via functions from the MASS package, using an example from the help page of nnet::multinom. The workhorse function used is MASS::dropterm:

> library(nnet)
> example(birthwt)
(bwt.mn <- multinom(low  ~ . , bwt) )


brthwt> bwt <- with(birthwt, {
brthwt+ race <- factor(race, labels = c("white", "black", "other"))
brthwt+ ptd <- factor(ptl > 0)
.
.
.
Call:
multinom(formula = low ~ ., data = bwt)

Coefficients:
(Intercept)         age         lwt   raceblack   raceother   smokeTRUE 
 0.82320102 -0.03723828 -0.01565359  1.19240391  0.74065606  0.75550487 
    ptdTRUE      htTRUE      uiTRUE        ftv1       ftv2+ 
 1.34375901  1.91320116  0.68020207 -0.43638470  0.17900392 

Residual Deviance: 195.4755 
AIC: 217.4755 
> confint(bwt.mn)   
                  2.5 %       97.5 %
(Intercept) -1.61649875  3.262900795
age         -0.11309745  0.038620883
lwt         -0.02953168 -0.001775495
raceblack    0.14190092  2.242906899
raceother   -0.16438896  1.645701076
smokeTRUE   -0.07755089  1.588560638
ptdTRUE      0.40173272  2.285785295
htTRUE       0.50053490  3.325867418
uiTRUE      -0.22990670  1.590310835
ftv1        -1.37601313  0.503243725
ftv2+       -0.71550657  1.073514417
> MASS::dropterm(bwt.mn, trace=FALSE, test="Chisq") 
# weights:  11 (10 variable)
initial  value 131.004817 
iter  10 value 98.297550
.
.
.
Single term deletions

Model:
low ~ age + lwt + race + smoke + ptd + ht + ui + ftv
       Df    AIC    LRT  Pr(Chi)   
<none>    217.48                   
age     1 216.42 0.9419 0.331796   
lwt     1 220.95 5.4739 0.019302 * 
race    2 219.23 5.7513 0.056380 . 
smoke   1 218.67 3.1982 0.073717 . 
ptd     1 223.58 8.1085 0.004406 **
ht      1 222.93 7.4584 0.006314 **
ui      1 217.59 2.1100 0.146342   
ftv     2 214.83 1.3582 0.507077   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

There are many other ways ... but using functions from MASS together with functions from nnet seems prudent.

Answer 4

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.

Answer 5

One way to think about the p-value is a likelihood test that your fit is better than some simpler fit with fewer terms (or possibly no terms, the constant fit). Below is some code.

# Multinomial fit
fit <- nnet::multinom(cyl ~ mpg + hp, data=datasets::mtcars)

# Multinomial fit with one or more terms dropped
base_fit <- nnet::multinom(cyl ~ 1, data=datasets::mtcars)
base_fit2 <- nnet::multinom(cyl ~ mpg, data=datasets::mtcars)

# p-value that the fit is better than the base_fit
result <- lmtest::lrtest(fit, base_fit)
p_val1 <- result$`Pr(>Chisq)`[[2]]

# p-value that the fit is better than the base_fit2
result <- lmtest::lrtest(fit, base_fit2)
p_val2 <- result$`Pr(>Chisq)`[[2]]

Results:

> p_val1
[1] 6.250054e-14
> p_val2
[1] 0.0003148036

Apparently dropping hp makes the fit worse (p=0.0003), and dropping both mpg and hp makes it much worse (p near zero).