Huge difference in result of vglm() and multinomial() for mlogit

Question

I am doing multinomial logistic regression model for iris data set,

library(VGAM)
mlog1 <- vglm(Species ~ ., data=iris, family=multinomial())
coef(mlog1)

and the coefficients are:

 (Intercept):1  (Intercept):2 Sepal.Length:1 Sepal.Length:2  Sepal.Width:1 
     34.243397      42.637804      10.746723       2.465220      12.815353 
 Sepal.Width:2 Petal.Length:1 Petal.Length:2  Petal.Width:1  Petal.Width:2 
      6.680887     -25.042636      -9.429385     -36.060294     -18.286137 

Then I use multinom() function and do the same thing:

library(nnet)
mlog2 <- multinom(Species ~ ., data=iris)

Coefficients:

Coefficients:
           (Intercept) Sepal.Length Sepal.Width Petal.Length Petal.Width
versicolor    18.69037    -5.458424   -8.707401     14.24477   -3.097684
virginica    -23.83628    -7.923634  -15.370769     23.65978   15.135301

It seems to be a big gap between these two results? Where did I do wrong? How can I fix them and get the similar result?

Forensic

Answer 1

The gap is due to two factors: (1) The multinomial() family in VGAM chooses the reference to be the last level of the response factor by default while multinom() in nnet always uses the first level as the reference. (2) The species categories in the iris data can be separated linearly, thus leading to very large coefficients and huge standard errors. Where exactly the numeric optimization stops when the log-likelihood virtually doesn't change further, can be somewhat different across implementations but is practically irrelevant.

As an example without separation consider a school choice regression model based on data from the German Socio-Economic Panel (1994-2002) in the AER package:

data("GSOEP9402", package = "AER")
library("nnet")
m1 <- multinom(school ~ meducation + memployment + log(income) + log(size),
  data = GSOEP9402)
m2 <- vglm(school ~ meducation + memployment + log(income) + log(size),
  data = GSOEP9402, family = multinomial(refLevel = 1))

Then, both models lead to esssentially the same coefficients:

coef(m1)
##                (Intercept) meducation memploymentparttime memploymentnone
## Realschule   -6.366449  0.3232377           0.4422277       0.7322972
## Gymnasium   -22.476933  0.6664295           0.8964440       1.0581122
##            log(income) log(size)
## Realschule   0.3877988 -1.297537
## Gymnasium    1.5347946 -1.757441

coef(m2, matrix = TRUE)
##                     log(mu[,2]/mu[,1]) log(mu[,3]/mu[,1])
## (Intercept)                 -6.3666257        -22.4778081
## meducation                   0.3232500          0.6664550
## memploymentparttime          0.4422720          0.8964986
## memploymentnone              0.7323156          1.0581625
## log(income)                  0.3877985          1.5348495
## log(size)                   -1.2975203         -1.7574912

Raglan