Model condition differences for ordinal Likert scale data

Modeling
,
1

Hi everyone,

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)?

Thanks!

2

Howdy! Have you tried conditional_effects(model, categorical=TRUE)?

3

Hi there,

I did and then I get this graph:

image
image717×441 13.7 KB

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:
image

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:

Hypotheses <- c(
"Music-NoMusic" = "ConditionCombinedMusic + ConditionCombinedNoMusic = 0",
"Music-Alone" = "ConditionCombinedMusic + Intercept[1] = 0",
"NoMusic-Alone" = "ConditionCombinedNoMusic + Intercept[1] = 0"
)

hypothesis(model, Hypotheses)

which yields:

image
image859×128 4.3 KB

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?

Thanks!

4

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:

s1 <- as_draws_df(model)
str(s1)
diff <- (s1$b_ConditionCombinedMusic - s1$b_ConditionCombinedNoMusic)
mean(diff)
quantile(diff, c(0.025, 0.975))

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.

5

Thanks for your reply!

I tried your suggestion, which yielded:

Estimate 1.907107
     2.5%     97.5% 
0.5380316 3.5037748

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:

Population-Level Effects: 
                        Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept[1]              -13.60      4.27   -22.96    -6.57 1.00    21826    12604
Intercept[2]              -11.02      3.42   -18.74    -5.41 1.00    31833    28161
Intercept[3]               -7.68      2.48   -13.34    -3.72 1.00    41753    48501
Intercept[4]               -3.28      1.12    -5.89    -1.55 1.00    46803    53743
Intercept[5]               -0.86      0.43    -1.84    -0.14 1.00    68207    59560
Intercept[6]               -0.63      0.40    -1.52     0.05 1.00    69734    58690
ConditionCombinedZalone     0.47      0.65    -0.73     1.84 1.00    58281    45437
ConditionCombinedZmusic     3.10      1.84     0.55     7.64 1.00    40228    27569

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?

Thanks!

6

Howdy,

It sounds like you are after something like the risk difference (the difference in probabilities) between the conditions. You can try use the {marginaleffects} package. An example of how you can use it is here: https://solomonkurz.netlify.app/blog/2023-05-21-causal-inference-with-ordinal-regression/

You’ll need the avg_contrasts() function.

7

Thanks Stephen!

I now used the avg_comparisons of the marginaleffects package:

comparisons_data <-avg_comparisons(model, variables = list(ConditionCombined = "pairwise"), type = "link", comparison = "difference")

print(comparisons_data)

This gets me the output I want, but I was wondering about some of the numbers within the contrasts:

image
image734×113 3.35 KB

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?

Thanks!

8

See reprex below. You did this?

library(brms)
str(kidney)
kidney$y <- rbinom(76, 4, p=0.2)
kidney$y <- factor(kidney$y, ordered=TRUE)

#model with reference level for "disease" = "other"
m1 <- brm(y ~ disease, data=kidney, family=cumulative("cloglog"), cores=4)
m1

#calculate contrast for diseaseGN vs diseaseAN, when the reference level is "other"
s1 <- as_draws_df(m1)
str(s1)
diff <- (s1$b_diseaseGN - s1$b_diseaseAN)
mean(diff)
quantile(diff, c(0.025, 0.975))


#compare to changing the reference level in R and re-running model 
#see coefficient for 'diseaseGN' below in m2, same as 'diff' calculated above for m1
kidney$disease <- relevel(kidney$disease, ref="AN")
m2 <- brm(y ~ disease, data=kidney, family=cumulative("cloglog"), cores=4)
m2
9

Hi everyone,

thanks to you both for your help!
No I get more or less the same outcome for all 3 methods:

Modeldraws <- as_draws_df(model)
str(Modeldraws)
diff <- Modeldraws <- (Modeldraws$b_ConditionCombinedMusic - Modeldraws$b_ConditionCombinedNoMusic)
mean(diff)
quantile(diff, c(0.025, 0.975))

gives

[1] 3.134915
> quantile(diff, c(0.025, 0.975))
     2.5%     97.5% 
0.5700196 7.7296409

Then

data$ConditionCombined <- relevel(factor(data$ConditionCombined), ref = "No Music")

Model2 <- brm(formula, data = data, family = cumulative("cauchit"), prior = prior, iter = 50000, warmup = 20000, chains = 3, cores = 3, save_pars = save_pars(all = TRUE))

gives

                                           estimate  lolim   uplim
ConditionCombinedMusic     3.12     0.55     7.73

and

comparisons_data <-avg_comparisons(model, variables = list(ConditionCombined = "pairwise"), type = "link", comparison = "difference")

print(comparisons_data)

gives

image
image748×125 3.44 KB

So for all 3 the results are now aligned.

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?

Thanks again!

10

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.

11

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.