[R] valid LRT between MASS::polr and nnet::multinom

Steve Taylor steve.taylor at aut.ac.nz
Wed Jul 8 06:18:58 CEST 2015


Dear R-helpers,

Does anyone know if the likelihoods calculated by these two packages are comparable in this way?  

That is, is this a valid likelihood ratio test?

# Reproducable example:
library(MASS)
library(nnet)
data(housing)
polr1 = MASS::polr(Sat ~ Infl + Type + Cont, weights=Freq, data=housing)
mnom1 = nnet::multinom(Sat ~ Infl + Type + Cont, weights=Freq, data=housing)
pll = logLik(polr1)
mll = logLik(mnom1)
res = data.frame(
  model = c('Proportional odds','Multinomial'),
  Function = c('MASS::polr','nnet::multinom'),
  nobs = c(attr(pll, 'nobs'), attr(mll, 'nobs')),
  df = c(attr(pll, 'df'), attr(mll, 'df')),
  logLik = c(pll,mll),
  deviance = c(deviance(polr1), deviance(mnom1)),
  AIC = c(AIC(polr1), AIC(mnom1)),
  stringsAsFactors = FALSE
)
res[3,1:2] = c("Difference","")
res[3,3:7] = apply(res[,3:7],2,diff)[1,]
print(res)
mytest = structure(
  list(
    statistic = setNames(res$logLik[3], "X-squared"),
    parameter = setNames(res$df[3],"df"),
    p.value = pchisq(res$logLik[3], res$df[3], lower.tail = FALSE),
    method = "Likelihood ratio test",
    data.name = "housing"
  ),
  class='htest'
)
print(mytest)

# If you want to see the fitted results:
library(effects)
plot(allEffects(polr1), layout=c(3,1), ylim=0:1)
plot(allEffects(mnom1), layout=c(3,1), ylim=0:1)

many thanks,
    Steve



More information about the R-help mailing list