Potential bug in R's `polr` function when run from a function environment?

Question

I may have found some sort of bug polr function (ordinal / polytomous regression) of the MASS library in R. The problem seems to be related to use of coef() on the summary object, but maybe is not.

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.

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 on a numerical dependent variable, i.e.

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 to run?

Answer 1

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)
}