Search code examples
rplotbayesianstanordinal

Brms: plotting three-way interaction in ordinal regression

(crossposting from mcstan)

Operating system: Ubuntu 18.04.5 LTS

brms version: 2.13.5

I have run a Bayesian ordinal regression using Buerkner's brms package (which provides a user-friendly interface to stan) and now am trying to plot the effect of three categorical predictors (Morphology, Cluster2, CountryExperiment) on the response variable (a Likert scale with 7 points). Following closely the information written by Paul Buerkner on conditional_effects for his package, I used the following code:

#model
fit_sc1 <- brm(
  formula = response ~ 1 + Cluster2*Morphology*CountryExperiment,
  data = datafile.for.model,
  family = cumulative("probit")
)
#plot
conditions <- make_conditions(fit_sc1, "Morphology")
conditional_effects(fit_sc1, "Cluster2:CountryExperiment", conditions = conditions, categorical=TRUE)

But I get the following warning:

Error: No valid effects detected.
In addition: Warning message:
Interactions cannot be plotted directly if 'categorical' is TRUE. Please use argument 'conditions' instead.

Note that the three-way interaction is plotted if I remove categorical=TRUE from the code above: the plot shows the three predictors in facets. But then the response variable is treated as a continuous variable (instead of ordinal) on the plot, and this is likely invalid for ordinal families.

enter image description here

I get a nice plot treating the response variable as ordinal, but only if I remove one of the effects. But then it's no longer a three-way interaction and this is really what I am interested in...

enter image description here

It looks like Buerkner's conditional_effects only works for linear regression three-way interactions... Has anyone managed to get around this issue? I suspect that it is possible to use the regression's posterior samples to do that and create a plot from scratch (that would be the best option), but I am not confident enough with the interpretation of intercepts in ordinal regression to do it all by myself.

Thanks!

Solution

  • I managed to get what I want using the ordinal regression's posterior samples to reconstruct the response's probabilities for all cells in the factorial design. I used Equation (5) p 79 in Buerkner & Vuorre's paper to find the probability of each of my 7 responses (along a 7-point Likert scale) for each of the 12 conditions (all combinations of my three variables). Equation (5) defines the probability P(Response=k) of each response category k as a function of (i) the 6 thresholds that split the underlying continuous scale, (ii) the predictor variables, and (iii) the corresponding regression coefficients. It was tedious but it works... It would be nice to develop a function that does that automatically.