24
$\begingroup$

I ran a multinomial logit model in JMP and got back results which included the AIC as well chi-squared p-values for each parameter estimate. The model has one categorical outcome and 7 categorical explanatory vars.

I then fit what I thought would build the same model in R, using the multinom function in the nnet package.

The code was basically:

fit1 <- multinom(y ~ x1 + x2 + ... xn, data=mydata);
summary(fit1);

However, the two give different results. With JMP the AIC is 2923.21, and with nnet::multinom the AIC is 3116.588.

So my first question is: Is one of the models wrong?

The second thing is, JMP gives chi-squared p-values for each parameter estimate, which I need. Running summary on the multinom fit1 does not - it just gives the estimates, AIC and Deviance.

My second question is thus: Is there a way to get the p-values for the model and estimates when using nnet::multinom?

I know mlogit is another R package for this and it looks like its output includes the p-values; however, I have not been able to run mlogit using my data. I think I had the data formatted right, but it said I had an invalid formula. I used the same formula that I used for multinom, but it seems like it requires a different format using a pipe and I don't understand how that works.

Improve this question
$\endgroup$
2
  • 2
    $\begingroup$ You can set the argument Hess=TRUE to get the Hessian back from multinom and then calculate the p-values manually. But I suggest you use the mlogit library (nnet can have convergence issues when the covariates are not scaled correctly). The vignettes for mlogit are pretty good, and should help you get your data set up correctly. The vignettes can be found at the usual place: cran.r-project.org/web/packages/mlogit $\endgroup$
    – Jason Morgan
    Commented Apr 18, 2011 at 23:26
  • $\begingroup$ I wrote up one way to think about the p-values in stats.stackexchange.com/a/502170. $\endgroup$
    – dfrankow
    Commented Dec 25, 2020 at 17:02

3 Answers 3

Reset to default
10
$\begingroup$

Im sure you've already found your solutions as this post is very old, but for those of us who are still looking for solutions - I have found Multinomial Probit and Logit Models in R is a great source for instructions on how to run a multinomial logistic regression model in R using mlogit package. If you go to the econometrics academy website she has all the scripts, data for R and SAS and STATA I think or SPSS one of those.

Which kind of explains how/why and what to do about transforming your data into the format of the "long" format vs "wide". Most likely you have a wide format, which requires transformation.

Multinomial Probit and Logit Models

Improve this answer
$\endgroup$
3
$\begingroup$

In general, differences in AIC values between two different pieces of software are not entirely surprising. Calculating the likelihoods often involves a constant that is the same between different models of the same data. Different developers can make different choices about what to leave in or out of that constant. When you should worry is when the differences in AIC values between two models differ. Actually I just noticed an argument to multinom() allows you to change how rows with identical X values are collapsed, and that this affects the baseline of the deviance, and hence the AIC. You could try different values of the summ argument and see if that makes the deviances agree. We don't know what JMP is doing! :)

If the estimated coefficients and standard errors are the same, then you're good. If the coefficients are not the same, don't forget that JMP might choose a different baseline outcome to calculate the coefficients for. multinom() makes different choices from mlogit(), for example.

Getting p-values from the summary() result of multinom() is pretty easy. I can't reproduce your models, so here's the example from the help page on multinom():

library("nnet")
data("Fishing", package = "mlogit")
fishing.mu <- multinom(mode ~ income, data = Fishing)
sum.fishing <- summary(fishing.mu) # gives a table of outcomes by covariates for coef and SE
str(sum.fishing)
# now get the p values by first getting the t values
pt(abs(sum.fishing$coefficients / sum.fishing$standard.errors),
  df=nrow(Fishing)-6,lower=FALSE)

I agree that figuring out the mlogit package is a bit of a challenge! Read the vignettes, carefully. They do help.

Improve this answer
$\endgroup$
6
  • $\begingroup$ How would I use the other (generic) variables from the Fishing dataset with multinom? $\endgroup$
    – gregmacfarlane
    Commented Apr 23, 2014 at 16:38
  • $\begingroup$ @gmacfarlane Just add the variables you want to the formula in multinom(mode ~ income + price.beach, ... $\endgroup$
    – atiretoo
    Commented May 26, 2014 at 16:25
  • $\begingroup$ @atiretoo I was looking for a way to get my pvals so thanks! but just for clarification - where does the 6 come from with the df? like what should I be counting to get my df? I have a continuous variable and a categorical variable (4 categories) in my model. So would that be a df of 5? Also, fishing is the entire data set right? Which is sample size-degree of freedom. $\endgroup$
    – Kerry
    Commented Jul 16, 2014 at 18:37
  • $\begingroup$ @atiretoo I apologize, I just found that if using package nnet to run the logit regression - it actually calculates the effective df and it is stored within the nnet object. So, I can extract the df from the model and use that. I have not checked if mlogit object will have the same information. $\endgroup$
    – Kerry
    Commented Jul 16, 2014 at 18:44
  • $\begingroup$ @atiretoo Any reason you estimate p-values with "pt" which use the T-distribution? $\endgroup$
    – eXpander
    Commented Oct 8, 2014 at 8:49
1
$\begingroup$

You could also try running a multinomial logit using the glmnet package. I'm not sure how to force it to keep all variables, but I'm sure it's possible.

Improve this answer
$\endgroup$

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