Finding confidence interval for response variables in ordinal logistic model

Question

While working with private data, I noticed that the ordinal logistic model fitted using the polr function from the MASS package, along with the confidence intervals provided by broom::tidy, does not display confidence intervals for the transitions between the response variable's categories. Instead, it returns an NA value for these entries in the table. How can I obtain confidence intervals for these transitions between categories? Any solution or did I misunderstand something about the ordinal logistic model? I noticed that in the example of the package function there is also no confidence interval for the transition of categories in the response variable.

I wanted to understand if confidence intervals for Low|Medium and Medium|High are possible to calculate on this exemple here:

#load libraries for models and data
library(MASS)

# fit model
fit <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)

# summarize model fit with tidiers
broom::tidy(fit, exponentiate = TRUE, conf.int = TRUE)

reference of exemple

Answer 1

It might be possible to fix the tidy.polr method (you could post an issue on the broom issues list), but in the meantime you can also do with this with ordinal::clm2() (which is actually more flexible than MASS::polr() too).

library(ordinal)
fit2 <- clm2(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
library(dplyr)  ## mutate, tibble, across
tidy.clm2 <- function(x,
                      conf.int = FALSE,
                      conf.level = 0.95,
                      exponentiate = FALSE,
                      p.values = FALSE,
                      ...) {
    cc <- coef(x)
    vv <- vcov(x)
    res <- tibble(term = names(cc),
                  estimate = cc,
                  std.error = sqrt(diag(vv)),
                  statistic = estimate/std.error)
    if (conf.int) {
        qq <- qnorm((conf.level + 1)/2)
        res <- mutate(res,
                      conf.low = estimate - qq*std.error,
                      conf.high = estimate + qq*std.error)
    }
    if (exponentiate) {
        res <- (res
            |> mutate(across(any_of(c("estimate", "conf.low", "conf.high")), exp))
            |> mutate(across(std.error, ~ . * estimate))
        )
    }
    return(res)
}

broom::tidy(fit2, exponentiate = TRUE, conf.int = TRUE)

The results look commensurate with what polr was giving you ...