I am trying to use conditional_effects() to visualize a 3-way interaction. My target plot would show how the size of the interaction between two variables (on the y-axis) varies as a function of the third variable (x-axis).

Unfortunately, the default way of plotting 3-way interactions with conditional_effects() does not quite do this. Here’s a reproducible example with the “epilepsy” data. The key effect is the interaction between the age of a patient and the two-way interaction between his number of visits to the hospital and treatment (zAge:visit:Trt).

Run model with “zAge” (continuous), “visit” and “Trt” (both discrete).

As you see, the plot puts the values of “visit” on different facets. Instead, I would like to show how the “visit:Trt” interaction grows/shrinks as a function of patients’ age, because this is the visualization that most straightforwardly maps to my research question.

Has anyone managed to achieve this kind of plot, either with conditional_effects() or with any other package? Any suggestions on how to achieve it?

Am I following correctly that you would like to display the interaction among the two predictors visit and Trt on the y-axis as a function of Age on the x-axis? If so, yes this is feasible.

Thanks! I’ll go through it to see if I can apply your solution to the reproducible example that I used in my post. In the meantime, if anyone else has suggestions please let me know :)

Not sure how to do it with conditional_effects(), but can’t this be done with fitted()?
Specifically, with fitted(…, summary=F) you’d then be able to compute a difference of differences of the posteriors (yielding your 2-way) at different values of the continuous predictor.

Hi Paul, thanks for your reply! But like I mentioned in my post, even after I use the conditions argument I can’t manage to make the 2-way interaction be plotted on the y-axis as a function of zAge. Could you say a bit more on how you think this is feasible?

I am sorry, I overlooked this in your original post. I should have read more thoroughly.

The problem is that brms does not plot effects directly but only the underlying variables or predictions. So it cannot display interaction effects on some axis (I understand this is what you want?).

You will have to create this plot manually unfortunately.

Yes, indeed, that is exactly what I want (and interaction effect on one of the axis). Ok, I think that I should be able to get the interaction by doing manual subtractions… I just wasn’t sure how to get the error/credible interval around it. Thanks Paul!

If you compute the interactions manually on a per-draw basis, you can obtain the CIs after doing all transformations via posterior_summary for example.

Oh, this is cool! Thanks for the suggestion, I didn’t know about emmeans!

I am trying to implement your ggplot() plotting option after having ran the brms model. Two follow-up questions:

ggplot() gives me a warning that lower.CL and upper.CL don’t exist. I think this might be related to running brms instead of lm, so it should be fixable by replacing these variables with lower.HDP and upper.HDP, right? This is a sample output from printing c_df after running the model in brms (sorry about the distorted alignment):

Yup - for mcmc models emmeans returns lower.HDP and upper.HDP.

I think the ribbon is pointy because HDI intervals aren’t parametric the same way CIs are - so they’re not as smooth. But that’s perfectly fine - it’s just your Bayesian showing :)

BTW, depending on the size of the data, using cov.red = unique might take a long time.
You can replace this with cov.red = function(x) {seq(min(x), max(x), length.out = k)} and adjust k (larger k, smoother ribbon (but will never be 100% smooth, due to point 2 above).

Another good option is the tidybayes package. Using the tidybayes package can give you a bit more control for handling different kinds of models – like zero-inflated or ordinal models with or without random effects, among many others – and how the output gets displayed. And for my way of thinking, the tidybayes syntax is a bit more explicit about what’s going on.

To provide quick and dirty usage for your example…
draws_interaction <- expand.grid(Trt=c(0,1), visit = c(1,2), zAge=seq(-1.75, 2.25, 0.05)) %>% add_fitted_draws(fit3way) %>% compare_levels(.value, by=Trt) %>% compare_levels(.value, by=visit) %>% median_hdi(.width=c(.5,.8,.95))

It’s worth pointing out that you have a good deal of flexibility about where and how you summarize the interaction this way. Only want to look at some values of zAge, change the values supplied within expand.grid(…). Prefer a different way of calculating the center or spread of the model predictions, use a different point_interval function. There’s a wormhole of customization to explore!

Well… that’s frustrating. I copy-pasted out of a working script on my end. Could be a package version issue? I am using tidybayes 1.1.0, but I wouldn’t have expected usage to change that much.

Yeah, it’s really puzzling. I was using the latest version of tidybayes (2.0.2) but I installed your version (1.1.0) to double-check. The error still persists with the older version.

I also thought that the problem might be that in your code Trt and visit are numeric, not categorical. But I converted them to factors prior to running compare_levels and the error is still there. So I don’t know why it’s not working…

One more thought… My code runs with warnings (that I’ve always ignored, since I can still get from A to B).

Warning messages:
1: unnest() has a new interface. See ?unnest for details.
Try `df %>% unnest(c(.lower, .upper))`, with `mutate()` if needed

unnest is a tidyr function that maybe called internally within compare_levels or median_hdi and it looks like the usage for that function has changed in the last handful of months (Sept. 2019). Perhaps tidybayes has not been updated to reflect the new usage? Regardless, it doesn’t seem like that provides an explanation for why the code works for me, but not you.