[R] bugs and misfeatures in polr(MASS).... fixed!

[tjb]

The start value generation code in polr is also known to fail quite
frequently. For example, against the Iris data as recently posted to this
list by blackscorpio (	Sep 6, 2010).

> polr(Species~Sepal_Length+Sepal_Width+Petal_Length+Petal_Width,data=iris)
Error in polr(Species ~ Sepal_Length + Sepal_Width + Petal_Length +
Petal_Width,  : 
  attempt to find suitable starting values failed
In addition: Warning messages:
1: In glm.fit(X, y1, wt, family = binomial(), offset = offset) :
  algorithm did not converge
2: In glm.fit(X, y1, wt, family = binomial(), offset = offset) :
  fitted probabilities numerically 0 or 1 occurred

I suggest that simply setting the coefficients beta to zero and the
cutpoints zeta to sensible values will always produce a feasible starting
point for non-pathological data. Here is my code that does this:

    if(missing(start)) {
      # try something that should always work -tjb
      u <- as.integer(table(y))
      u <- (cumsum(u)/sum(u))[1:q]
      zetas <-
         switch(method,
                "logistic"= qlogis(u),
                "probit"=   qnorm(u),
                "cauchit"=  qcauchy(u),
                "cloglog"=  -log(-log(u)) )
      s0 <- c(rep(0,pc),zetas[1],log(diff(zetas)))

Using this start code the problem is not manifested. 

> source('fixed-polr.R')
> polr(Species~Sepal_Length+Sepal_Width+Petal_Length+Petal_Width,data=iris)
Call:
polr(formula = Species ~ Sepal_Length + Sepal_Width + Petal_Length + 
    Petal_Width, data = iris)

Coefficients:
Sepal_Length  Sepal_Width Petal_Length  Petal_Width 
   -2.466724    -6.671515     9.431689    18.270058 

Intercepts:
   setosa|versicolor versicolor|virginica 
            4.080189            42.639320 

Residual Deviance: 11.89857 
AIC: 23.89857 

My change would also likely fix the problem reported by Kevin Coombes on May
6, 2010.
-- 
View this message in context: http://r.789695.n4.nabble.com/bugs-and-misfeatures-in-polr-MASS-fixed-tp3024677p3030405.html
Sent from the R help mailing list archive at Nabble.com.