Error that profile() on a model generated with MASS::polr() cannot find data when used in a function

Question

I am looking to use profile() to isolate the summary data from a model made with MASS::polr() inside a function. It works outside of a function, and the function works if using glm() instead of polr(). Can someone explain why and point me in the right direction to get this to work?

Reproducible Data

# for polr() example
df2 <- data.frame(
  a = ordered(c(1,2,3,4,5,1,2,3,4,5,1,2,3,4,5,1,2,3,4,5), levels = 1:5),
  b = ordered(c(1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4), levels = 1:4))

# for glm() example
df3 <- data.frame(
  a = c(1,2,3,4,5,1,2,3,4,5,1,2,3,4,5,1,2,3,4,5),
  b = c(1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4))

This is what I want to put in a function to return the summary data.

model <- polr(a ~ b, data = df2)
summary <- attr(profile(model), "summary")
summary

What have I tried?

fun_polr <- function(x){
  model <- polr(a ~ b, data = x)
  summary <- attr(profile(model), "summary")
  return(summary)
}

When I run fun_polr(df2), I get the error "Error in eval(expr, p) : object 'x' not found."

If I run a similar function with glm() on the df3 data set, I don't get an error.

fun_glm <- function(x){
  model <- glm(a ~ b, data = x)
  summary <- attr(profile(model), "summary")
  return(summary)
}

If I create the function without profile() to simply return the model, it runs fine. So I conclude it is how profile() interacts with polr() that is causing me the problem. fun_polr2(df2) runs without any problem.

fun_polr2 <- function(x){
  model <- polr(a ~ b, data = x)
  return(model)
}

Further insight I also found that profile(polr.model) did not work when the data fed into the model was generated by a function when creating the data frame.

set.seed(123)
df <- data.frame(
  a = ordered(sample(1:5, 100, replace = TRUE), levels = 1:5),
  b = ordered(sample(1:4, 100, replace = TRUE), levels = 1:4))

model <- polr(a ~ b, data = df)
summary <- attr(profile(model), "summary")

The model is generated but the profile() step results in an error, "'data' must be a data.frame, environment, or list." I tried the same approach with generated data for a glm() model and didn't get the same problem.

df4 <- data.frame(
  a = sample(1:5, 100, replace = TRUE),
  b = sample(1:4, 100, replace = TRUE))

model <- glm(a ~ b, data = df4)
summary <- attr(profile(model), "summary")

I expect these observations are related, but am unsure how to resolve. I would appreciate any explanations or help. Might there be any alternatives to using profile() to isolate the summary data?

Answer 1

The workaround is to add Hess = TRUE to the function call.

fun_polr <- function(x){
  model <- polr(a ~ b, data = x, Hess = TRUE)
  summary <- attr(profile(model), "summary")
  return(summary)
}

The default is Hess = FALSE. The issue is this code that is part of the vcov.polr() called by summary.polr().

if (is.null(object$Hessian)) {
    message("\nRe-fitting to get Hessian\n")
    utils::flush.console()
    object <- update(object, Hess = TRUE, start = c(object$coefficients, 
                                                    object$zeta))

If Hess is FALSE, then update() is used to re-run the model. The "object" only has the value of "x" for the data, and is thus kicking it out. Setting the value to TRUE from the beginning will avoid the call to update().