Re: [R] valid LRT between MASS::polr and nnet::multinom
Dear Steve, The short answer is no, but the test you propose is in my experience usually a close approximation to a valid test. The proportional-odds and multinomial-logistic regression models differ in two ways: the po model has an equal-coefficients (parallel-regressions) assumptions (except for intercepts) and so has fewer parameters; the po model models cumulative logits, while the mnl model models individual-category logits. As a consequence of the latter difference, the po model is not a specialization of the mnl model, as required by the LR test. A proper test of the po (equal-coefficients) assumption is to fit a cumulative-logit model with this constraint. You can do this with the vglm() function in the VGAM package using the cumulative family. Here is an example: -- snip - library(effects) # for WVS data library(nnet) # for multinom() library(MASS) # for polr() library(VGAM) # for vglm() Loading required package: stats4 Loading required package: splines mod.polr - polr(poverty ~ gender + religion + degree + country*poly(age,3), +data=WVS) coef(mod.polr) gendermale religionyes degreeyes 0.1691953 0.1684846 0.1413380 countryNorway countrySweden countryUSA -0.3217821 -0.5732783 0.6040006 poly(age, 3)1 poly(age, 3)2 poly(age, 3)3 19.9101983 -10.2208416 6.1157062 countryNorway:poly(age, 3)1 countrySweden:poly(age, 3)1countryUSA:poly(age, 3)1 -17.0042706 -9.4160841 1.5577738 countryNorway:poly(age, 3)2 countrySweden:poly(age, 3)2countryUSA:poly(age, 3)2 17.3824147 17.3856575 10.1575695 countryNorway:poly(age, 3)3 countrySweden:poly(age, 3)3countryUSA:poly(age, 3)3 3.5181428 2.3652443 -8.4004861 logLik(mod.polr) 'log Lik.' -5182.602 (df=20) mod.vglm.p - vgam(poverty ~ gender + religion + degree + country*poly(age,3), +data=WVS, family=cumulative(parallel=TRUE)) coef(mod.vglm.p) # same within rounding as polr except for sign (Intercept):1 (Intercept):2 gendermale 0.2139034 2.0267774 -0.1691716 religionyes degreeyes countryNorway -0.1684541 -0.1413316 0.3218158 countrySweden countryUSA poly(age, 3)1 0.5733987 -0.6040348 -19.9340368 poly(age, 3)2 poly(age, 3)3 countryNorway:poly(age, 3)1 10.2120529 -6.1558215 17.0535729 countrySweden:poly(age, 3)1countryUSA:poly(age, 3)1 countryNorway:poly(age, 3)2 9.4712722 -1.5238183 -17.3344400 countrySweden:poly(age, 3)2countryUSA:poly(age, 3)2 countryNorway:poly(age, 3)3 -17.3422095 -10.1469673 -3.4087977 countrySweden:poly(age, 3)3countryUSA:poly(age, 3)3 -2.2649306 8.4423869 logLik(mod.vglm.p) # same (within rounding error) [1] -5182.601 mod.multinom - multinom(poverty ~ gender + religion + degree + country*poly(age,3), + data=WVS) # weights: 60 (38 variable) initial value 5911.632725 iter 10 value 5162.749080 iter 20 value 5007.247684 iter 30 value 4995.375350 iter 40 value 4989.909216 iter 50 value 4987.650806 iter 60 value 4987.190310 iter 70 value 4987.131548 iter 80 value 4987.075422 final value 4987.073531 converged logLik(mod.multinom) 'log Lik.' -4987.074 (df=38) length(coef(mod.multinom)) [1] 38 mod.vglm.np - vgam(poverty ~ gender + religion + degree + country*poly(age,3), +data=WVS, family=cumulative(parallel=FALSE)) logLik(mod.vglm.np) # close but not the same as multinom [1] -4988.865 length(coef(mod.vglm.np)) # same no. of coefs [1] 38 snip -- Best, John John Fox, Professor McMaster University Hamilton, Ontario, Canada http://socserv.mcmaster.ca/jfox/ On Wed, 8 Jul 2015 04:18:58 + Steve Taylor steve.tay...@aut.ac.nz wrote: 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 =
[R] valid LRT between MASS::polr and nnet::multinom
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
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
