[R] convergence for proportional odds model

David Duffy David.Duffy at qimr.edu.au
Wed Sep 7 04:39:29 CEST 2005


liu abc <liu2074 at yahoo.com> wrote:
>
> I am using proportional odds model for ordinal responses in
> dose-response experiments. For some samll data, SAS can successfully
> provide estimators of the parameters, but the built-in function polr()
> in R fails. Would you like to tell me how to make some change so I
> can use polr() to obtain the estimators? Or anyone can give me a hint
> about the conditions for the existance of MLE in such a simple case?
> By the way, for the variable "resp" which must be ordered factor, how
> can I do it? Thanks a lot.
>
> Guohui

> The following is one example I used both in SAS and R.
>
> in R:
>
> library(MASS)
> dose.resp = matrix( c(1,1,1,1,2,2,2,3,3,3, 2,2,3,3,4,4,5,4,5,5), ncol=2)
> colnames(dose.resp)= c("resp", "dose")
> polr( factor(resp, ordered=T)~dose, data=dose.resp)
> #Error in optim(start, fmin, gmin, method = "BFGS", hessian = Hess, ...) :
> # initial value in 'vmmin' is not finite

It seems to be the starting values.  Using lrm() from the Design package gave

> dose.resp <- as.data.frame(dose.resp)
> dose.resp$resp <- factor(dose.resp$resp)
> library(Design)
> lrm(resp ~ dose, data=dose.resp)

       Obs  Max Deriv Model L.R.       d.f.          P          C        Dxy
        10      6e-06      11.43          1      7e-04      0.909      0.818
     Gamma      Tau-a         R2      Brier
     0.931        0.6      0.768      0.014

     Coef    S.E.  Wald Z P
y>=2 -10.904 5.137 -2.12  0.0338
y>=3 -14.336 6.287 -2.28  0.0226
dose   3.160 1.399  2.26  0.0239

and giving polr starting values:


> print(m1 <- polr(resp ~ dose, data=dose.resp, start=c(-1, -4, 3)))
Call:
polr(formula = resp ~ dose, data = dose.resp, start = c(-1, -4,
    3))

Coefficients:
    dose
3.158911

Intercepts:
     1|2      2|3
10.90172 14.33296

Residual Deviance: 10.34367
AIC: 16.34367

Even then, summary(m1) gives the same problem (as it refits).  There is
separation in the data, of course, but I presume the ordinality gives
some extra information.

David Duffy.




More information about the R-help mailing list