Re: [R] difference between MASS::polr() and Design::lrm()

From: Rune Haubo <rune.haubo_at_gmail.com>
Date: Mon, 30 Jun 2008 13:58:52 +0200

Dear Vito

No, you are not wrong, but you should center score prior to model estimation:

summary(fm1 <- polr(factor(grade)~I(score - mean(score))))

which gives the same standard errors as do lrm. Now the intercepts refer the median score rather than some potential unrealistic score of 0.

You can always check if standard errors are appropriate by plotting the profile likelihood and the normal approximation (based on the standard errors):

x <- seq(-.043 - .04, -.043 + .04, len = 20) dev <- double(20)
for(i in 1:20)
dev[i] <- polr(factor(grade) ~ 1 + offset(x[i] * (score - median(score))))$deviance

yy <- spline(x, dev, n = 100, method = "natural") plot(yy[[1]], exp(-(yy[[2]] - min(yy[[2]]))/2), type = "l",

        xlab = "Score", ylab = "Normalized profile likelihood") ## Add Normal approximation:
lines(yy[[1]], exp(-(yy[[1]] - coef(fm1))^2 /

                 (2 * vcov(fm1)[1, 1])), lty = 2,
      col = "red")

## Add confidence limits:
level <- c(0.95, 0.99)
lim <- sapply(level, function(x)

              exp(-qchisq(x, df=1)/2) )
abline(h = lim, col = "grey")

## Add confidence interval points:
points(confint(fm1), rep(lim[1], 2), pch = 3)

So these standard errors are the most appropriate you can get, but because the profile likelihood is not symmetric, you could consider reporting the asymmetric profile likelihood confidence interval instead.

Hope this helps
Rune

2008/6/30 vito muggeo <vmuggeo_at_dssm.unipa.it>:
> Dear all,
> It appears that MASS::polr() and Design::lrm() return the same point
> estimates but different st.errs when fitting proportional odds models,
>
> grade<-c(4,4,2,4,3,2,3,1,3,3,2,2,3,3,2,4,2,4,5,2,1,4,1,2,5,3,4,2,2,1)
> score<-c(525,533,545,582,581,576,572,609,559,543,576,525,574,582,574,471,595,
> 557,557,584,599,517,649,584,463,591,488,563,553,549)
>
> library(MASS)
> library(Design)
>
> summary(polr(factor(grade)~score))
> lrm(factor(grade)~score)
>
> It seems to me that results from lrm() are right..
>
> Am I wrong?
> Thanks,
> vito
>
> --
> ====================================
> Vito M.R. Muggeo
> Dip.to Sc Statist e Matem `Vianelli'
> UniversitÓ di Palermo
> viale delle Scienze, edificio 13
> 90128 Palermo - ITALY
> tel: 091 6626240
> fax: 091 485726/485612
> http://dssm.unipa.it/vmuggeo
>
> ______________________________________________
> 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.
>



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. Received on Mon 30 Jun 2008 - 12:05:12 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Mon 30 Jun 2008 - 12:30:55 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive