I am writing to follow up on my previous request for assistance with modeling ordinal data and plotting interactions using BRMS (Bayesian Regression Models using Stan). I have made progress in some areas but still don’t have a perfect solution.

On the positive side, I was able to create a model in the form

`df$catDist <- interaction(df$distortion, df$category)`

which allowed me to visualize all interactions effectively using:

`pp_check(brm3, type = "bars_grouped", group = "catDist ", ndraws = 500, prob = 0.95)`

Nonetheless, this model is different from the desired distortion* category model.

Another way, I have tried has been to create predictions for the different interactions myself. However, I have not had too much luck with this compared to the above.

```
# Step 1: Create the new data frame with all combinations of `imAge`, `category`, and `distortion`
new_data <- expand.grid(
imAge = levels(df$imAge),
category = levels(df$category),
distortion = levels(df$distortion)
)
# Step 2: Get the expected predicted values and uncertainty intervals for each combination
predictions <- posterior_epred(brm2, newdata = new_data, ndraws = 500)
# Step 3: Normalize the probabilities for each combination of `imAge`, `category`, and `distortion`
normalized_probs <- apply(predictions, c(1, 2), function(x) x / sum(x))
# Step 4: Calculate the `predicted_rating`, `Q2.5`, and `Q97.5` within each `imAge`, `category`, and `distortion` combination
new_data <- new_data %>%
group_by(imAge, category, distortion) %>%
mutate(
rating_probs = list(colMeans(normalized_probs[imAge, category, distortion, , drop = FALSE])),
predicted_rating = sum(seq_len(5) * rating_probs[[1]]), # Calculate the weighted average
Q2.5 = sum(ifelse(seq_len(5) <= quantile(predictions[imAge, category, distortion, , drop = FALSE], prob = 0.025), rating_probs[[1]], 0)),
Q97.5 = sum(ifelse(seq_len(5) >= quantile(predictions[imAge, category, distortion, , drop = FALSE], prob = 0.975), rating_probs[[1]], 0))
) %>%
ungroup() # Remove the grouping for the resulting data frame
# Reshape the data for plotting
plot_data <- new_data %>%
mutate(rating = 1:5) %>%
pivot_longer(cols = c("Q2.5", "predicted_rating", "Q97.5"), names_to = "interval", values_to = "value") %>%
mutate(interval = factor(interval, levels = c("Q2.5", "predicted_rating", "Q97.5")),
interval_label = ifelse(interval == "predicted_rating", "Predicted Rating", "Uncertainty Interval"))
# Step 5: Create the plot with error bars
ggplot(plot_data, aes(x = factor(rating), y = value, fill = imAge, group = interaction(imAge, category, distortion))) +
geom_bar(stat = "identity", position = position_dodge(width = 0.9), alpha = 0.7, width = 0.6) +
geom_errorbar(aes(ymin = Q2.5, ymax = Q97.5), width = 0.2, position = position_dodge(width = 0.9), alpha = 0.5) +
facet_wrap(~ category + distortion, ncol = 1) +
theme_classic() +
labs(x = "Rating", y = "Probability", title = "Predicted Probabilities with Uncertainty Intervals") +
scale_y_continuous(expand = c(0, 0), limits = c(0, max(plot_data$value) + 0.05)) +
theme(axis.ticks.x = element_blank(),
legend.position = "right",
axis.text = element_text(size = 12),
axis.title = element_text(size = 14))
```

I am thus still facing challenges in visualising the interactions for the desired distortion* category model. Any insights or alternative approaches would be highly appreciated. Please feel free to share any ideas or suggestions to help me move forward.

Best regards,

Simon