I have been struggling with how to compare the effect of conditions on my DV. I have used

‘’‘stan
model ← brm(formula, data = data, family = cumulative(“cloglog”), prior = prior, iter = 5000, warmup = 2000, chains = 3, cores = 3, save_pars = save_pars(all = TRUE))
‘’’
with the formula being Playdate ~ConditionCombined

The model seems to fit, I get a summary, and I also managed to look at conditional effects using condition_effects(model).

Ultimately, what I want is to see if participants respond differently to this item depending on the condition they are in. I suppose I can’t just use conditional_effects, because then the data is treated as continuous. How can I see which condition is different from each of the other conditions (i.e., the ordinal equivalent of simple contrasts)?

But then I am still not sure how to quantify how the effect differs across conditions. I thought there was some way to get the ordinal regression coefficient + 95% HDI’s?

Summary(model) gives me:

But then I only get an estimate of the difference between the Intercept (the Alone condition) and the Music and No Music conditions respectively, but not the difference between the Music and No Music condition relative to one another.

I tried an approach coined by @matti in a previous post using the hypotheses:

And thus gives slightly different conclusions for the 2 contrasts that were in the summary. Did I misunderstand how the summary table works? Or did I incorrectly specify the contrasts?

I’m not familiar with what the hypothesis() method in brms does, but you should simply be able to extract the samples for the coefficients and take the appropriate difference to get any contrast that you want on the scale of the linear predictor.
In summary(model) you have the Music vs Alone and No Music vs Alone. If you want Music vs No Music then you can extract the samples for each coefficient and take the difference and then summarize the resulting vector of samples. It will be approximately 1.23 - -0.68=1.91 for Music vs No Music, when No Music is the reference level. You would do something like:

You can check by changing the reference level in R in your data to “No Music” instead of “Alone” and then re-running the model on the data with the updated reference level.

However, when I ran the analysis again with No Music as a reference factor (I did not know how to change the reference group quickly, so I just put a Z in front of it, because it seems to order alphabetically) I got:

So the bottom line should be Music compared to the reference group No music, but these numbers are different. So something different seems to be happening when calculating the contrast compared to rerunning the analysis with a different reference factor?

Given that I am using a 1-7 Likert scale, would these numbers refer to actual scores on the latent variable? And because this variable is continuous, is that why some of the estimation can exceed 7 even though 7 in the highest value in the Likert scale?

I guess the only question now is, where do the values higher than 7 come from? Are they direct estimates of the difference between the 2 conditions? Because my data is on a 1-7 scale, so should 6 then not be the max value? Or is the model indeed assuming a latent continuous distribution that can surpass 7 from time to time?

Based on your first post, it looks like you are using the complementary-log-log link function for a cumulative distribution. All of the results presented here in terms of the coefficients are on the scale of the linear predictor, so they are on the scale of complementary-log-log. The cumulative distribution assumes the outcome variable arises by categorizing a latent continuous variable. Is there some particular reason you are using the complementary-log-log link function? Something like cumulative("probit") would be easier to interpret. See https://journals.sagepub.com/doi/10.1177/2515245918823199 for a nice intro to ordinal models in brms.

Ah, I see.
Since the beginning of this post, I realized that the gauchit link function seems to work best for my model (i.e., model converges with the least amount of divergent transitions, and the distribution of the different integers in the posterior predictive checks looks most like the data).

But I agree that things are perhaps a bit difficult to interpret. Is there perhaps a way to estimate a more conventional effect size based on these contrasts for readers who are less familiar with these kinds of methods?

And thanks for the paper reference!
It really helped, not just with this analysis, but I also realized that another model I am running should actually be a sequential (instead of a cumulative) model.