I was hoping to have a surface plot using the output from marginal_effects (brms). My response is ordinal and fitted using cumulative family. My plan is to present the predicted probability for each response category as a surface plot. However, when I set ordinal=T option, marginal_effects only have two elements and the covriate1:covaraite2 object is not available. Is there any way to get this?

m<-marginal_effects(fit,surface=T)-> return three data frames (covaraite1:, covaraite2: and covaraite1:covaraite2:)

m<-marginal_effects(fit,ordinal=T, surface=T)-> return only two data frames

Due to the nature of the ordinal = TRUE plot, we cannot show a second covariate in the same way as we do otherwise for surface plots since one of the two axes of the plot is already taken by the ordinal categories.

I recommend using the conditions argument to generate a facet plot with different values of the second covariate. See ?marginal_effects for examples.

I am not successful using â€śconditionsâ€ť for now. I see the difficulty, but it would have been great to see how P(Y=k) for a specific response category k change as a function of two covaraiates in a model fitted using te(covariate1,covaraite 2) and then we can compare this with the result for P(Y=kâ€™) for a different response category kâ€™ presented in a different surface plot.

This should result in a 3-D array of (summarized) predictions, which you can convert to a long data.frame containing the necessary variables. This data.frame can then be used to create the desired plot with ggplot2 using geom_raster and facet_wrap(<response category>)