How can I make my brms partial proportional odds model more "flexible?"

In response to @harrelfe : Here is a comparison with a maximum likelihood estimate from the ordinal package. A partial proportional odds estimate is implemented in the clmm2() function. It doesn’t have the ability to fit random slopes or multiple levels of random intercept but otherwise is close to what I am trying to fit with brms. As you can see in the plot, the maximum likelihood estimates (triangles) are much better than the brms estimates (circles). Code is below.

code to produce plot

(follows example code in original post)

library(ordinal)

fit_ppo_clmm2 <- clmm2(Rating ~ Variety,
                       nominal = ~ Variety,
                       random = Grower,
                       weights = N_Plants,
                       link = 'logistic',
                       threshold = 'flexible',
                       data = exdata,
                       Hess = TRUE)

pred_grid <- with(exdata, expand.grid(Variety = unique(Variety), Rating = unique(Rating)))

# Predict means for each variety x rating
pred_clmm2 <- predict(fit_ppo_clmm2, newdata = pred_grid)

pred_brm <- predict(fit_ppo, newdata = data.frame(Variety = 1:3), re_formula = ~ 0)

pred_values <- cbind(pred_grid, clmm2 = pred_clmm2, brm = c(pred_brm[, 3:1])) %>%
  tidyr::pivot_longer(c(clmm2, brm), names_to = 'model', values_to = 'pred')

ggplot(observed_p, aes(x = Rating, color = Rating)) +
  geom_point(aes(y = p), alpha = 0.5) +
  geom_point(aes(y = pred, shape = model, group = model), data = pred_values, position = position_dodge(width = 0.5), size = 3) +
  facet_wrap(~ Variety) +
  theme_bw() + see::scale_color_okabeito()