tomka - 1 year ago 47

R Question

I may have found some sort of bug

`polr`

`MASS`

`R`

`coef()`

`summary`

The problem occurs in a function of type:

`pol_me <- function(d){`

some_x <- d[,1]

mod <- polr(some_x ~ d[,2])

pol_sum <- summary(mod)

return(pol_sum)

}

To illustrate, I simulate some data for an ordinal regression model.

`library(MASS)`

set.seed(2016)

n=1000

x1 <- rnorm(n)

x2 <- 2*x1 +rnorm(n)

make_ord <- function(y){

y_ord <- y

y_ord[y < (-1)] <- 1

y_ord[y >= (-1) & y <= (1)] <- 2

y_ord[y >= 1] <- 3

y_ord <- as.factor(y_ord)

}

x1 <- make_ord(x1)

dat <- data.frame(x1,x2)

When we now call the function:

`library(MASS)`

pol_me(d = dat)

We get error

`Error in eval(expr, envir, enclos) : object 'some_x' not found`

I do not think this should logically happen at this point. In fact when we define alternative function in which the model command is replaced by a linear model

`lm`

`mod <- lm(as.numeric(some_x) ~ d[,2])`

The resulting function works fine.

Is this really a bug or a programming problem in my code and how can I get

`pol_me`

`polr`

`summary`

`coef`

`pol_me`

Answer Source

`summary(polr(dat[,1] ~ dat[,2]))`

returns semi-error message `Re-fitting to get Hessian`

and it's the cause of the error. `polr`

's argument `Hess = T`

will solve your problem. (`?polr`

says *Hess: logical for whether the Hessian (the observed information matrix) should be returned. Use this if you intend to call summary or vcov on the fit.*)

```
pol_me <- function(d){
some_x <- d[,1]
mod <- polr(some_x ~ d[,2], Hess = T) # modify
pol_sum <- summary(mod)
return(pol_sum)
}
```