Finding confidence interval for response variables in ordinal logistic model

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

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