Testing proportional odds assumption in R for ordinal logistic regression on ordinal variable with many categories

Question

I'm trying to run an ordinal logistic regression to predict an ordinal variable with 8 categories (bins of # of days spent recreating) based on a categorical (BIPOC/not) and continuous (urbanicity score) variable. I'm having no trouble with the model itself, which I ran using the polr() function from the MASS package in R with the following syntax:

polr(days ~ urbanicity*BIPOC, data = dat1, Hess = T)

However, I'm having trouble figuring out how to test whether my model meets the proportional odds assumption.

I've heard that statistical tests for proportional odds, such as the Brant's test, are too stringent, so I was looking into ways to do graphical assumption testing. However, most of the examples I found online (such as the great resource from https://stats.oarc.ucla.edu/r/dae/ordinal-logistic-regression/) are only for response variables with 3 categories, where to test proportional odds you essentially comparing 2 binomial models with 2 different cutpoints against each other.

Does anyone have experience adapting this type of testing for variables with >3 categories, where one would have to make multiple comparisons?

My current work-around, based on https://peopleanalytics-regression-book.org/ord-reg.html, is to run my own series of binomial models and compare the coefficients against each other. But I'm finding it challenging to determine when two coefficients are sufficiently different that the proportional odds assumption should be considered violated, so other solutions would be much appreciated!

Reference code for my work-around:

m.res <- list()

for (i in 2:8){
  # create binary variable bsed on cutoff
  dat$x <- ifelse(as.numeric(dat$days) < i, 0, 1)
  
  # model based on cutoff
  m <- glm(
    x ~ urbanicity*BIPOC, 
    data = dat, 
    family = "binomial"
  ) 
  
  # store model results in a list
  m.res[[i-1]] <- m
  
}

coefficient_comparison <- data.frame(matrix(nrow=0, ncol=7,))
colnames(coefficient_comparison) <- c("y>1", "y>2", "y>3", "y>4",
                                      "y>5", "y>6", "y>7")

coefficient_comparison <- data.frame(
  oneplus = summary(m.res[[1]])$coefficients[ , "Estimate"],
  twoplus = summary(m.res[[2]])$coefficients[ , "Estimate"],
  threeplus = summary(m.res[[3]])$coefficients[ , "Estimate"],
  fourplus = summary(m.res[[4]])$coefficients[ , "Estimate"],
  fiveplus = summary(m.res[[5]])$coefficients[ , "Estimate"],
  sixplus = summary(m.res[[6]])$coefficients[ , "Estimate"],
  sevenplus = summary(m.res[[7]])$coefficients[ , "Estimate"]
)

coefficient_comparison

Which gives me a table of coefficients for each binomial model, that I can then compare across models:

                  oneplus    twoplus  threeplus    fourplus   fiveplus    sixplus  sevenplus
(Intercept)     3.2413652  0.6931362  0.2917739 -0.01815633 -0.2713922 -0.4564793 -0.6188225
urbanicity     -1.2027440 -1.0300582 -1.2566895 -1.34407778 -1.4980431 -1.6129901 -1.7764729
BIPOC1         -0.7400516 -0.7444939 -0.8024623 -0.85788852 -0.8748921 -0.8564442 -0.8467173
urbanic:BIPOC1  1.2562572  0.8580332  0.6188747  0.47969936  0.3371632  0.2803162  0.1873982