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)1 countryUSA:poly(age,
3)1
-17.0042706 -9.4160841
1.5577738
countryNorway:poly(age, 3)2 countrySweden:poly(age, 3)2 countryUSA:poly(age,
3)2
17.3824147 17.3856575
10.1575695
countryNorway:poly(age, 3)3 countrySweden:poly(age, 3)3 countryUSA: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)1 countryUSA:poly(age, 3)1 countryNorway:poly(age,
3)2
9.4712722 -1.5238183
-17.3344400
countrySweden:poly(age, 3)2 countryUSA:poly(age, 3)2 countryNorway:poly(age,
3)3
-17.3422095 -10.1469673
-3.4087977
countrySweden:poly(age, 3)3 countryUSA: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 +0000
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 = 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.
______________________________________________
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.
