[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