Getting fitted values using the emmeans and predict functions

Question

I want to get the difference between the "average" scores on a five-point scale using the emmeans library. However, I found that this is only possible for the models of the ordinal library. For the mgcv library, we can only get an approximate result (I'm not sure if this is correct).

I'll give you an example.

# load libraries
library(mgcv)
library(brms)
library(ordinal)
library(emmeans)

# create toy example
y <- c(
  1,2,3,2,1,4,5,
  3,4,5,4,3,1,2
)
ord <- ordered(y)
x <- factor(rep(c('a', 'b'), each = 7))
df <- data.frame(x, y, ord)
item <- 1:5
ndf <- data.frame(x = c('a', 'b')) # for newdata

# create models
mClm <- clm(ord~x, data = df)
mGam <- gam(y~x, family = ocat(R = 5), data = df)
mBrms <- brm(ord~x, family = cumulative("logit"), data = df)

# very close results
predict(mClm, type = 'prob', newdata = ndf)[['fit']] %*% item
predict(mGam, type = 'response', newdata = ndf) %*% item
predict(mBrms, newdata = ndf) %*% item

# latent variable vs "mean.class"
emmeans(mClm, ~x, mode = 'mean.class')
emmeans(mClm, ~x)
emmeans(mGam, ~x)
emmeans(mBrms, ~x)

# is a solution?
predict(mGam, newdata = ndf)
theta <- mGam$family$getTheta(T)
qnorm(plogis(predict(mGam, newdata = ndf), location = mean(theta)), mean = mean(item))

We see that through the predict function we can get a correctly averaged score using three different packages. However, such features are not available in the emmeans functionality.

Are there ways to get the source data for the mgcv (ocat family) and brm (cumulative) models?