Residual deviance difference between multinom() and vglm() function I have used two functions multinom() in package nnet and vglm() in package (VGAM) to make a multinomial logistic regression. 
# packages
require(nnet)
require(VGAM)
require(foreign)

# load data
gator <- read.table("https://onedrive.live.com/redir?resid=EF059954D82A6592!1855&authkey=!AISIAhhB66eHDZ0&ithint=file%2ctxt",
  header = TRUE)

# fit model with multinom () function
fit.multinom <- multinom(cbind(Fish, Bird, Invertebrate, Reptile, Other) ~ Lake + Size,
  data = gator)

# fit model with vglm() function
fit.vglm <- vglm(cbind(Bird, Invertebrate, Reptile, Other, Fish) ~ Lake + Size,
  data = gator, family = multinomial)

# model summary
summary(fit.multinom)
summary(fit.vglm)

The coefficients obtained from these two models are similar. However, the residual deviance are by far different (540 and 52, respectively). 
But, these above differences only observed for response variable in matrix form. When response variable is presented in single column, both approach gives the same residual deviance (359.96) and AIC (375.96).
# load data
ml <- read.dta("http://www.ats.ucla.edu/stat/data/hsbdemo.dta")

# model construction
ml$prog2 <- relevel(ml$prog, ref = "academic")
test.multinom <- multinom(prog2 ~ ses + write, data = ml)
test.vglm <- vglm(prog2 ~ ses + write, data = ml, family="multinomial")

# model summary
summary(test.multinom)
summary(test.vglm)

Could you help me explain why they are different?
 A: The difference is just by a constant. The VGAM package takes the multinomial coefficients into account while the nnet package does not. Note that this does not affect the analysis of deviance or information criteria etc. The multinomial coefficients stay the same (for all models on this data set) - so they just cancel out when looking at differences in deviance or log-likelihood.
I illustrate this with the log-likelihood below. The difference is:
logLik(fit.vglm)
## [1] -74.42948
logLik(fit.multinom)
## 'log Lik.' -270.0401 (df=20)
logLik(fit.vglm) - logLik(fit.multinom)
## 'log Lik.' 195.6107 (df=20)

The multinom version can be easily computed by hand, based on the predicted log-probabilities p and observed frequencies y for each category and observation:
p <- predict(fit.multinom, type = "prob")
y <- model.response(model.frame(fit.multinom))
sum(log(p) * y)
## [1] -270.0401

The vglm version adds the logarithm of the multinomial coefficient $n_i! / (y_{i1}! \dots y_{i5}!)$. The log-factorial $\log(n!)$ can be easily computed in R by lgamma(n + 1). Hence we can easily comput the log-coefficient for each observation and then sum up:
lcf <- apply(y, 1, function(x) lgamma(sum(x) + 1) - sum(lgamma(x + 1)))
sum(lcf)
## [1] 195.6107

This is exactly the difference between the log-likelihoods as computed by the two packages. And this shows up again in the deviances.
