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

[Rune Haubo]

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