For your first question, conditional_effects
actually returns a list of data frames, one for each effect
. When you plot them, you get a list of ggplots, which you can then modify to add connecting lines for discrete variables. To add the connecting lines, we just need to look at the structure of the data frames returned by conditional_effects
so we can see what variables we need to use for our call to geom_line
.
I wasn’t able to get your code to run without error, so I’ll use the iris
example to demonstrate. In order to make this relatively general, I’ll create a second categorical variable that will become the faceting variable in the conditional_effects
plots.
library(tidyverse)
library(patchwork) # To lay out multiple plots
library(brms)
# Add an additional grouping variable to iris
set.seed(5)
iris = iris %>% mutate(new.group = sample(c("A","B"), nrow(iris), replace=TRUE))
m = brm(Petal.Width ~ Species + Petal.Length + new.group, data=iris,
cores=4, silent=2, refresh=0)
Now we’ll generate conditional effects for Species
and Petal.Length
by new.group
.
ce = conditional_effects(m, effects=c("Species", "Petal.Length"),
conditions=make_conditions(m, "new.group"))
As you can see below, ce
is a list of two data frames, one for each of the two effects in our call to conditional_effects
.
str(ce)
#> List of 2
#> $ Species :'data.frame': 6 obs. of 10 variables:
#> ..$ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 2 3 1 2 3
#> ..$ Petal.Width : num [1:6] 1.2 1.2 1.2 1.2 1.2 ...
#> ..$ Petal.Length: num [1:6] 3.76 3.76 3.76 3.76 3.76 ...
#> ..$ new.group : Factor w/ 2 levels "A","B": 1 1 1 2 2 2
#> ..$ cond__ : Factor w/ 2 levels "new.group = A",..: 1 1 1 2 2 2
#> ..$ effect1__ : Factor w/ 3 levels "setosa","versicolor",..: 1 2 3 1 2 3
#> ..$ estimate__ : num [1:6] 0.758 1.198 1.605 0.783 1.223 ...
#> ..$ se__ : num [1:6] 0.0888 0.035 0.0679 0.0842 0.0334 ...
#> ..$ lower__ : num [1:6] 0.574 1.128 1.467 0.605 1.153 ...
#> ..$ upper__ : num [1:6] 0.928 1.262 1.744 0.951 1.291 ...
#> ..- attr(*, "effects")= chr "Species"
#> ..- attr(*, "response")= chr "Petal.Width"
#> ..- attr(*, "surface")= logi FALSE
#> ..- attr(*, "categorical")= logi FALSE
#> ..- attr(*, "ordinal")= logi FALSE
#> ..- attr(*, "points")='data.frame': 150 obs. of 4 variables:
#> .. ..$ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
#> .. ..$ resp__ : num [1:150] 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.1 0.2 0.2 ...
#> .. ..$ cond__ : Factor w/ 2 levels "new.group = A",..: 1 1 1 1 1 1 1 1 1 1 ...
#> .. ..$ effect1__: Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
#> $ Petal.Length:'data.frame': 200 obs. of 10 variables:
#> ..$ Petal.Length: num [1:200] 1 1.06 1.12 1.18 1.24 ...
#> ..$ Petal.Width : num [1:200] 1.2 1.2 1.2 1.2 1.2 ...
#> ..$ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
#> ..$ new.group : Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
#> ..$ cond__ : Factor w/ 2 levels "new.group = A",..: 1 1 1 1 1 1 1 1 1 1 ...
#> ..$ effect1__ : num [1:200] 1 1.06 1.12 1.18 1.24 ...
#> ..$ estimate__ : num [1:200] 0.129 0.143 0.156 0.169 0.183 ...
#> ..$ se__ : num [1:200] 0.0332 0.0326 0.0316 0.0311 0.0302 ...
#> ..$ lower__ : num [1:200] 0.0635 0.0786 0.0933 0.1076 0.1222 ...
#> ..$ upper__ : num [1:200] 0.193 0.205 0.217 0.229 0.243 ...
#> ..- attr(*, "effects")= chr "Petal.Length"
#> ..- attr(*, "response")= chr "Petal.Width"
#> ..- attr(*, "surface")= logi FALSE
#> ..- attr(*, "categorical")= logi FALSE
#> ..- attr(*, "ordinal")= logi FALSE
#> ..- attr(*, "points")='data.frame': 50 obs. of 4 variables:
#> .. ..$ Petal.Length: num [1:50] 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.5 1.5 1.6 ...
#> .. ..$ resp__ : num [1:50] 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.1 0.2 0.2 ...
#> .. ..$ cond__ : Factor w/ 2 levels "new.group = A",..: 1 1 1 1 1 1 1 1 1 1 ...
#> .. ..$ effect1__ : num [1:50] 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.5 1.5 1.6 ...
#> - attr(*, "class")= chr "brms_conditional_effects"
Add we can also see the plots that are generated from each of these data frames:
wrap_plots(plot(ce, ask=FALSE, plot=FALSE), ncol=1)

Let’s look at the structure of the Species
data frame (see below). Note that there’s a Species
column and an estimate__
column, which we’ll use below for the x and y axes when we add a line to the plot (if we wish, we can also add lines for the upper and lower credible intervals, using the lower__
and upper__
columns). We also need to group by cond__
so that a separate line will be drawn for each facet.
# Look at data frame for Species
ce$Species
#> Species Petal.Width Petal.Length new.group cond__ effect1__
#> 1 setosa 1.199333 3.758 A new.group = A setosa
#> 2 versicolor 1.199333 3.758 A new.group = A versicolor
#> 3 virginica 1.199333 3.758 A new.group = A virginica
#> 4 setosa 1.199333 3.758 B new.group = B setosa
#> 5 versicolor 1.199333 3.758 B new.group = B versicolor
#> 6 virginica 1.199333 3.758 B new.group = B virginica
#> estimate__ se__ lower__ upper__
#> 1 0.7575476 0.08883262 0.5741506 0.9282307
#> 2 1.1980489 0.03501293 1.1282554 1.2622561
#> 3 1.6054027 0.06787588 1.4674149 1.7437657
#> 4 0.7828697 0.08415274 0.6054553 0.9508139
#> 5 1.2229951 0.03340769 1.1534383 1.2910627
#> 6 1.6309308 0.07029002 1.4909029 1.7744132
# Regenerate the plots
p.ce2 = plot(ce, plot=FALSE)
# Add a line to the Species plot
p.ce2[["Species"]] = p.ce2[["Species"]] +
geom_line(aes(Species, estimate__, group=cond__), colour="red")
wrap_plots(p.ce2, ncol=1)

For your second question, as @frank_hezemans noted, it sounds like a marginal effects plot is what you’re looking for. In addition to emmeans
, the marginaleffects
package is another option.