[R] difference between MASS::polr() and Design::lrm()
rune.haubo at gmail.com
Mon Jun 30 16:13:00 CEST 2008
(question to the authors of MASS below)
2008/6/30 vito muggeo <vmuggeo at dssm.unipa.it>:
> Dear Haubo,
> many thanks for your reply.
> Yes you are right, by scaling the "score", I get the same results.
> However it sounds strange to me. I understand that the SE and/or t-ratio of
> the intercepts depend on the scale of the explanatory variable, but I do not
> understand why the t-ratio of the slope also depends on the scale.. In
> classical regression models (LM, GLM, ...) the scale of the covariate does
> not affect the t ratio! As far as I know, classical references (e.g.
> Agresti, 2002) do not discuss a such point. Could yu suggest me additional
> reference? Namely why should I center the explanatory variables?
The wrong standard errors is a computational issue. Centering is a
statistical and (sometimes) a computational issue.
The t-value is just the ratio of the estimate to its standard error,
and when the program does not give you the correct standard error,
naturally the t-value is also wrong. It is often a good idea to center
covariates for inferential purposes and a very good reference is
Venables' "Exegeses on Linear Models":
I find it natural to expect polr to give the correct standard error
without centering the covariate. I do not know why it does not, but
perhaps is has something to do with the scaling of the parameters
parsed to the BFGS-optimizer in optim. Optimization of the intercepts
seems to be performed on something like
summary(fm1 <- polr(factor(grade)~score))
which seems to distort the whole variance-covariance matrix of the
parameters and therefore provides wrong standard errors for the
score-parameter as well. Optimization is much better scaled if you do
summary(fm2 <- polr(factor(grade)~1 + I(score - median(score))))
This is merely a guess, and perhaps the authors of MASS, could provide
One other thing, that puzzles me, is the fact that the parameter
estimate/standard error ratio is reported as a t-value rather than a
z-value. Likelihood theory (as I know it) seems to suggest a z-value
rather than the t-value. A t-value is also provided by the special
case; the ordinary logistic regression model:
summary(glm(grade > 3 ~ score, family = binomial))
as well as by lrm in package Design.
> Of course profile-Lik based CI may be very usefuls at this aim..but this is
> another story..
I agree, but the topic is closely related to standard errors ;-)
> many thanks again,
> Rune Haubo ha scritto:
>> Dear Vito
>> No, you are not wrong, but you should center score prior to model
>> 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
>> 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 -
>> yy <- spline(x, dev, n = 100, method = "natural")
>> plot(yy[], exp(-(yy[] - min(yy[]))/2), type = "l",
>> xlab = "Score", ylab = "Normalized profile likelihood")
>> ## Add Normal approximation:
>> lines(yy[], exp(-(yy[] - 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, 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
>> Hope this helps
>> 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,
>>> It seems to me that results from lrm() are right..
>>> Am I wrong?
>>> 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
>>> R-help at r-project.org mailing list
>>> PLEASE do read the posting guide
>>> and provide commented, minimal, self-contained, reproducible code.
> 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
Rune Haubo Bojesen Christensen
Master Student, M.Sc.
Phone: (+45) 30 26 45 54
Informatics and Mathematical Modelling
Section for Statistics
Technical University of Denmark, Build.321, DK-2800 Kgs. Lyngby, Denmark
More information about the R-help