tomka tomka - 3 months ago 13
R Question

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

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.

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

Edit note: Since first submission I edited this post tracing down the problem to the
polr
,
summary
or
coef
function called in the environment of
pol_me
. When called in the global environment, the problem (obviously) does not occur.

Answer

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