How to set up and estimate a multinomial logit model in R? 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.
 A: 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. 
A: 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
A: 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.
