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

John Fox jfox at mcmaster.ca
Wed Jul 8 15:25:10 CEST 2015


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.taylor at 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 at 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.



More information about the R-help mailing list