Ordinal Regression - conditional effects discrepancy

I have fit a cumulative ordinal model using brms. The code for model formulation is below. The issue I have is with conditional_effects().

My understanding is that conditional_effects() by default will plot the mean of the posterior samples of the expected value/mean of the posterior predictive distribution as calculated by posterior_epred(). However if I calculate directly the means using posterior_epred() and summarising they differ from the estimates__ provided in conditional_effects().

I believe these estimates are the probabilities of having an outcome at that particular category and should sum to 1 across all categories for given conditions, they do from posterior_epred() but not from conditional_effects().

This is not the case for a gaussian family mixed effects model for example, where the estimates were the same subject to sampling error. Is there an issue with the estimates being calculated for ordinal models in conditional_effects() or am I misunderstanding what is being calculated here.

I am not able to provide the data to replicate the example as it arises from a clinical trial, I apologise. The snippets of code hopefully provide some insight into how one would recreate this, newdata would need to match the conditions within conditional_effects()

Any comments/help would be greatly appreciated

Thanks,
Dan

model_fit <- brm(WHOScore ~ (1 + time | TNO) + CareStatusID_txt + Age +
                           time + TrtID_txt + TrtID_txt:time, data = model_data, 
                           family=cumulative("logit"), iter = 5000, file = 'modelfit')

conditions <- make_conditions(model_fit , c('CareStatusID_txt', 'time'))

conditional_effects(model_fit, effects = TrtID_txt, conditions, categorical = TRUE)

posterior_epred(model_fit, newdata) %>% data.frame() %>% summarise_all(mean)

@paul.buerkner is there a bug with conditional_effects estimates for ordinal regression models?, the estimates don’t match those calculated directly with posterior_epred

So I have resolved the discrepancy, no bug!, I have realised that by default conditional_effects return the median of the samples rather than the mean, this can be changed by changing robust to FALSE. It also be default does not include group-level effects, again resolved by setting re_formula to NULL. This will result in matching estimates from posterior_epred used above.

My question now is around the appropriateness of median or mean. The means seem to have the properties that they partition 0-1 for the category probabilities as one would expect for estimates for given conditions. They also seem to be much more inline with the raw category proportions when looked at for the conditions.

The medians however do not sum to 1 and can be drastically different to the means by up to 0.2 in my model above.

For ordinal models the mean seems to be preferable. Would be interested to hear any others thoughts on this.

Mean sounds sensible in this case. I choose the median by default in conditional_effects because sometimes the mean behaves weirdly relative to the credible interval boundaries and then I get lots of questions of people suspecting there is a problem. So I go with median which behaves consistently with the credible intervals as all are based on quantiles.

Thanks Paul, appreciate you taking the time to respond. In this instance I have chosen to push ahead with the Mean.