How do I get p-values using the
multinom function of
nnet package in
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
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")
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
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.