ordinal::clmm supports fitting models to ordinal data with (more that one) random effects. The package's predict function also seems to incorporate estimated values of random effects (random slopes as well as random intercepts) to generate unit-level predictions; if the function's argument newdata contains the data used to fit the model, the predictions are generated for an average unit i.e., random effects are not used and only the fixed-effects part of the model is used for predictions.

When one wishes to use emmeans::emmeans then to get predictions over a reference grid generated using clmm model, does bias-adjustment need to be performed on model predictions obtained for an "average" participant? For example, please consider the following scenario:

The target $Y$ is an ordinal variable with 5 categories. The predictors are $X_{1}$ (categorical variable with two categories), $X_{2}$ (ordinal variable with 3 categories) and $X_3$ (continuous, bounded predictor). For each participant (ID), I make several observations of $Y$ with $X_1 \in$ $\{A, B\}$, $X_{2} \in \{low, med, high\}$; $X_3 = K$ across all observations for the participant. I fit the following cumulative link mixed model to check if the predictors are associated with the target using the logit link function:

fm.model <- Y ~ X1 * X2 * X3 + (1 + X1 | ID)

Assuming that the model fits the data well, I want to estimate predictions, for an average ID, averaged over $X_2$ levels for the following combinations of the predictors:

model.emm <- emmeans(fm.model, ~ X3 | X1, at = list(X1 = c("A", "B"), X3 = c(0.1, 0.7)))
model.emm@grid

| X3  | X1 |
| --- |----|
| 0.1 | A  |
| 0.7 | A  |
| 0.1 | B  |
| 0.7 | B  |

Following an example from the emmeans vignette, I create bias-adjusted weights as follows:

V <- VarCorr(fm.model)$id

sig <- sapply(c(0, 1), function(x) {
    a <- c(1, x)
    sqrt(sum(a * V %*% a))
})

# a vector of four values (a, a, b, b), a > b
SIG <- rep(sig, each = 2)  

I then use SIG vector to obtain bias-adjusted predictions for each of the $Y$ categories as follows:

emmeans(fm.model, ~ Y | X1 + X3, mode = "prob", 
    at = list(X1 = c("A", "B"), X3 = c(0.1, 0.7)), sigma = SIG, bias.adj = T)