xinyi lin
2008-Jan-05 13:14 UTC
[R] Likelihood ratio test for proportional odds logistic regression
Hi, I want to do a global likelihood ratio test for the proportional odds logistic regression model and am unsure how to go about it. I am using the polr() function in library(MASS). 1. Is the p-value from the likelihood ratio test obtained by anova(fit1,fit2), where fit1 is the polr model with only the intercept and fit2 is the full polr model (refer to example below)? So in the case of the example below, the p-value would be 1. 2. For the model in which there is only one independent variable, I would expect the Wald test and the likelihood ratio test to give similar p-values. However the p-values obtained from anova(fit1,fit3) (refer to example below) are very different (0.0002622986 vs. 1). Why is this so?> library(MASS) > fit1 <- polr(housing$Sat~1) > fit2<- polr(housing$Sat~housing$Infl) > fit3<- polr(housing$Sat~housing$Cont) > summary(fit1)Re-fitting to get Hessian Call: polr(formula = housing$Sat ~ 1) No coefficients Intercepts: Value Std. Error t value Low|Medium -0.6931 0.2500 -2.7726 Medium|High 0.6931 0.2500 2.7726 Residual Deviance: 158.2002 AIC: 162.2002> summary(fit2)Re-fitting to get Hessian Call: polr(formula = housing$Sat ~ housing$Infl) Coefficients: Value Std. Error t value housing$InflMedium 6.347464e-06 0.5303301 1.196889e-05 housing$InflHigh 6.347464e-06 0.5303301 1.196889e-05 Intercepts: Value Std. Error t value Low|Medium -0.6931 0.3953 -1.7535 Medium|High 0.6932 0.3953 1.7536 Residual Deviance: 158.2002 AIC: 166.2002> summary(fit3)Re-fitting to get Hessian Call: polr(formula = housing$Sat ~ housing$Cont) Coefficients: Value Std. Error t value housing$ContHigh 0.0001135777 0.4330091 0.0002622986 Intercepts: Value Std. Error t value Low|Medium -0.6931 0.3307 -2.0956 Medium|High 0.6932 0.3307 2.0960 Residual Deviance: 158.2002 AIC: 164.2002> anova(fit1,fit2)Likelihood ratio tests of ordinal regression models Response: housing$Sat Model Resid. df Resid. Dev Test Df LR stat. Pr(Chi) 1 1 70 158.2002 2 housing$Infl 68 158.2002 1 vs 2 2 -6.375558e-10 1> anova(fit1,fit3)Likelihood ratio tests of ordinal regression models Response: housing$Sat Model Resid. df Resid. Dev Test Df LR stat. Pr(Chi) 1 1 70 158.2002 2 housing$Cont 69 158.2002 1 vs 2 1 -1.224427e-07 1 Thank you, Xinyi
Prof Brian Ripley
2008-Jan-05 13:53 UTC
[R] Likelihood ratio test for proportional odds logistic regression
On Sat, 5 Jan 2008, xinyi lin wrote:> Hi, > > I want to do a global likelihood ratio test for the proportional odds > logistic regression model and am unsure how to go about it. I am using > the polr() function in library(MASS). > > 1. Is the p-value from the likelihood ratio test obtained by > anova(fit1,fit2), where fit1 is the polr model with only the intercept > and fit2 is the full polr model (refer to example below)? So in the > case of the example below, the p-value would be 1.There is no improvement in fit, as the near-zero coefficients show. You are not calling polr correctly on this example: 'why is this so?' given that it *is* the example on the help page and all you had to do was to read the help (or the book for which this is support software).> 2. For the model in which there is only one independent variable, I > would expect the Wald test and the likelihood ratio test to give > similar p-values. However the p-values obtained from anova(fit1,fit3) > (refer to example below) are very different (0.0002622986 vs. 1). Why > is this so?Because you compared a t-value to a p-value, not at all the same thing.> > >> library(MASS) >> fit1 <- polr(housing$Sat~1) >> fit2<- polr(housing$Sat~housing$Infl) >> fit3<- polr(housing$Sat~housing$Cont) >> summary(fit1) > > Re-fitting to get Hessian > > Call: > polr(formula = housing$Sat ~ 1) > > No coefficients > > Intercepts: > Value Std. Error t value > Low|Medium -0.6931 0.2500 -2.7726 > Medium|High 0.6931 0.2500 2.7726 > > Residual Deviance: 158.2002 > AIC: 162.2002 >> summary(fit2) > > Re-fitting to get Hessian > > Call: > polr(formula = housing$Sat ~ housing$Infl) > > Coefficients: > Value Std. Error t value > housing$InflMedium 6.347464e-06 0.5303301 1.196889e-05 > housing$InflHigh 6.347464e-06 0.5303301 1.196889e-05 > > Intercepts: > Value Std. Error t value > Low|Medium -0.6931 0.3953 -1.7535 > Medium|High 0.6932 0.3953 1.7536 > > Residual Deviance: 158.2002 > AIC: 166.2002 >> summary(fit3) > > Re-fitting to get Hessian > > Call: > polr(formula = housing$Sat ~ housing$Cont) > > Coefficients: > Value Std. Error t value > housing$ContHigh 0.0001135777 0.4330091 0.0002622986 > > Intercepts: > Value Std. Error t value > Low|Medium -0.6931 0.3307 -2.0956 > Medium|High 0.6932 0.3307 2.0960 > > Residual Deviance: 158.2002 > AIC: 164.2002 >> anova(fit1,fit2) > Likelihood ratio tests of ordinal regression models > > Response: housing$Sat > Model Resid. df Resid. Dev Test Df LR stat. Pr(Chi) > 1 1 70 158.2002 > 2 housing$Infl 68 158.2002 1 vs 2 2 -6.375558e-10 1 >> anova(fit1,fit3) > Likelihood ratio tests of ordinal regression models > > Response: housing$Sat > Model Resid. df Resid. Dev Test Df LR stat. Pr(Chi) > 1 1 70 158.2002 > 2 housing$Cont 69 158.2002 1 vs 2 1 -1.224427e-07 1 > > > Thank you, > Xinyi > > ______________________________________________ > R-help at r-project.org mailing list > 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. >-- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595