R help - Sep 2005 - convergence for proportional odds model

Hey, everyone,
 
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")
dose.resp
#> dose.resp
#      resp dose
# [1,]    1    2
# [2,]    1    2
# [3,]    1    3
# [4,]    1    3
# [5,]    2    4
# [6,]    2    4
# [7,]    2    5
# [8,]    3    4
# [9,]    3    5
#[10,]    3    5
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
#In addition: Warning message:
#fitted probabilities numerically 0 or 1 occurred in: 
#glm.fit(X, y1, wt, family = binomial(), offset = offset) 
 
in SAS
NOTE: PROC LOGISTIC is fitting the cumulative logit model. The probabilities
      modeled are summed over the responses having the lower Ordered Values in
      the Response Profile table.
NOTE: Convergence criterion (GCONV=1E-8) satisfied.
NOTE: There were 10 observations read from the data set WORK.ONE.


		
---------------------------------
 Click here to donate to the Hurricane Katrina relief effort.
	[[alternative HTML version deleted]]

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.