How to make the plot of both levels of variable with conditional_effects

Hi,
I am running an ordinal model with brms, the model is here:

com.df$GAP <- factor(com.df$GAP, levels = c("g200", "g350", "g500",
                                            "g650", "g800", "g950"))
contrasts(com.df$GAP) <- contr.sdif
com.df$Experiment <- factor(com.df$Experiment)
com.df$Rating <- factor(com.df$Rating, ordered = TRUE)
com.Main.model <- brm(Rating ~ 1 + GAP + Experiment + (1 | sub) + (1 | item),
                      data = com.df,
                      init = 0,
                      family = cumulative("probit", threshold = "flexible"),
                      chains = 4, 
                      warmup = 1000,
                      iter = 5000,
                      thin = 1,
                      cores = 4)
com.Main <- summary(com.Main.model)

Now, i want to plot the probability of participants’ rating on each scores(i.e., 1, 2 and 3 in this model). i run the following code:

lhce <- conditional_effects(com.Main.model, 
                            "GAP", 
                            categorical = TRUE)

but i found that the probability value is only about “Experiment1(实验1)”, so i wonder to get the value of second level of Experiment.

Thanks.

If possible, add also code to simulate data or attach a (subset of) the dataset you work with.

Please also provide the following information in addition to your question:

  • Operating System: Windows
  • brms Version:2.20.4

Welcome wenbo!

conditional_effects uses the reference level for factor variables (and the mean for numeric variables) when generating conditional effects for other variables. You can have conditional_effects produce results for all levels of a factor by using the conditions argument. Here’s an example:

library(brms)

m = brm(Petal.Width ~ Species + Petal.Length, data=iris)

ce = conditional_effects(m, "Petal.Length",
                         conditions=make_conditions(m, "Species"))

ce

So, in your case, I think the code would be:

lhce <- conditional_effects(com.Main.model, 
                            "GAP", 
                            conditions=makes_conditions(com.Main.model, "Experiment"),
                            categorical = TRUE)

Thanks for your code.

i think i nearly finished my figures.

but i still to make it more better, i want to combine these points together like this, is it possible?

lhce <- conditional_effects(com.Main.model, "GAP", 
                            conditions = make_conditions(com.Main.model,
                                                         "Experiment"),
                            categorical = TRUE)

fig1F <- plot(lhce, plot = FALSE)[[1]]
com.df$Rating <- as.numeric(com.df$Rating)
fig1F1 <- fig1F + 
  scale_color_lancet(labels = c("高分", "中分", "低分"))+
  scale_fill_lancet(labels = c("高分", "中分", "低分"))+
  theme_test()+
  theme(legend.position = "top",
        legend.title = element_blank(),
        legend.text = element_text(size = 12))+
  labs(x = NULL, y = "概率")+
  scale_y_continuous(limits = c(0.1, 0.7))+
  theme(axis.title.y = element_text(margin = margin(r = 8), 
                                    size = 12, 
                                    family = "songti",
                                    color = "black"),
        axis.title.x = element_text(margin = margin(t = 8), 
                                    size = 12, 
                                    family = "songti",
                                    color = "black"),
        axis.text.x = element_text(size = 12, 
                                   color = "black", 
                                   family = "TNM"),
        axis.text.y = element_text(size = 12, 
                                   color = "black", 
                                   family = "TNM"))
fig1F1

Then, i still want to make a figure which cross two levels of Experiment, like we usually do in ANOVA test: Main effect, is it possible?

Great thanks for your help.

Wenbo Yu

I’m not sure I understand what you’re asking. First, in the plot image, it looks like you’ve added the red curves by hand. Are those what you want to add to the plot? GAP is coded as a factor in your data. Is it really a continuous variable rather than a set of categories? If so, then maybe it would make sense to model it that way and then conditional_effects would automatically plot it as a curve, rather than as discrete categories. But if GAP really is discrete, I’m not sure it would make sense to connect the points with a line. Also, if GAP really is categorical, it looks like the categories are ordered, so maybe it would make sense to recode GAP as an ordered factor, rather than the current coding which is unordered.

Second, regarding the ANOVA figure you asked about, can you draw a picture to show how you’d like the plot to look?

Thanks for your help.
For the first question,
GAP in my study is indeed a category factor.
I agree with your point of view: linking categorical variables together is not reasonable and the reason why I want to connect the dots is to better showcase the trend of change.

If it is hard to achieve, i may plot the line with PS.
Thanks.

here is my cleaned data and code
data.csv (122.3 KB)

PRIOR:
prior.2 <- c(prior(normal(0, 0.5), class = "b", coef = "GAPg350Mg200"),
             prior(normal(-0.2, 0.5), class = "b", coef = "GAPg500Mg350"),
             prior(normal(-0.2, 0.5), class = "b", coef = "GAPg650Mg500"),
             prior(normal(-0.2, 0.3), class = "b", coef = "GAPg800Mg650"),
             prior(normal(0, 0.5), class = "b", coef = "GAPg950Mg800"),
             prior(normal(-0.3, 0.5), class = "b", coef = "Experiment实验2"))

MODEL:
com.df$GAP <- factor(com.df$GAP, levels = c("g200", "g350", "g500",
                                            "g650", "g800", "g950"))
contrasts(com.df$GAP) <- contr.sdif
com.df$Experiment <- factor(com.df$Experiment)
com.df$Rating <- factor(com.df$Rating, ordered = TRUE)
com.Main.model <- brm(Rating ~ 1 + GAP + Experiment + (1 | sub) + (1 | item),
                      data = com.df,
                      prior = prior.2,
                      sample_prior = "yes",
                      init = 0,
                      family = cumulative("probit", threshold = "flexible"),
                      chains = 4, 
                      warmup = 1000,
                      iter = 5000,
                      thin = 1,
                      cores = 4)
com.Main <- summary(com.Main.model)

FIGURE:
lhce <- conditional_effects(com.Main.model, "GAP", 
                            conditions = make_conditions(com.Main.model,
                                                         "Experiment"),
                            categorical = TRUE)

fig1F <- plot(lhce, plot = FALSE)[[1]]
com.df$Rating <- as.numeric(com.df$Rating)
fig1F1 <- fig1F + 
  scale_color_lancet(labels = c("高分", "中分", "低分"))+
  scale_fill_lancet(labels = c("高分", "中分", "低分"))+
  theme_test()+
  theme(legend.position = "top",
        legend.title = element_blank(),
        legend.text = element_text(size = 12))+
  labs(x = NULL, y = "概率")+
  scale_y_continuous(limits = c(0.1, 0.7))+
  theme(axis.title.y = element_text(margin = margin(r = 8), 
                                    size = 12, 
                                    family = "songti",
                                    color = "black"),
        axis.title.x = element_text(margin = margin(t = 8), 
                                    size = 12, 
                                    family = "songti",
                                    color = "black"),
        axis.text.x = element_text(size = 12, 
                                   color = "black", 
                                   family = "TNM"),
        axis.text.y = element_text(size = 12, 
                                   color = "black", 
                                   family = "TNM"),
        strip.text = element_text(size = 12, family = "songti", color = "black"),
        strip.background = element_rect(fill = "transparent"))
fig1F1

For the second question,
i just want to plot a main effect, which average the “概率” in experimnt 1 and experiment 2. just like this:

i try to obtain the estimates of each rating value across each GAP levels, but conditional_effect only offer results of experiment 1, just the reference level.

Thanks

You might be interested in using the emmeans package to obtain marginal mean effects for one or more variables of interest, averaging across the remaining variables in your model. brms models are directly supported in the emmeans package. You can then extract the posterior draws corresponding to the emmeans object (for example using the function gather_emmeans_draws from the tidybayes package) and then visualise those draws however you like (for example using ggplot2 with the help of ggdist).

Unfortunately I don’t have an opportunity right now to look at your specific data, but this blog post by Andrew Heiss contains lots of useful examples and code.

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.

Completely agree with the advice. OP probably wants something like

... +
  geom_line(
    aes(y = estimate__, group = cond__), 
    position = position_dodge(width = 0.4)
  )
1 Like

Thank you @Ax3man. I forgot about the dodging for multiple lines.

1 Like

thanks for you help.

i followed your advice, but there is still a bug, that is, what i want to line is not the independent variable. Look at the figure i showed, there are three “概率”(probability in English) in each GAP condition. It means participants chose three ratings (高分、中分、低分) with different probility. So it could not work when i set group = Experiment or group = GAP.

I try to set group = Rating, but the figure is not what i want.

thanks for you suggestion

The general idea is that the group aesthetic to geom_line needs to create a “group” of points that we want to connect with a line. In this case, each group will be a facet and an outcome category (I forgot about the outcome categories in my previous code). The data frame produced by conditional_effects uses the column name cats__ for the outcome category column. So, the code in your case would be:

lhce <- conditional_effects(com.Main.model, "GAP", 
                            conditions = make_conditions(com.Main.model,
                                                         "Experiment"),
                            categorical = TRUE)

# Look at the data frame produced by conditional effects
lhce$`GAP:cats__`

fig1F <- plot(lhce, plot = FALSE)[[1]]

# Add the connecting lines
fig1F + 
  geom_line(aes(group=interaction(cond__, cats__)),
            position=position_dodge(0.4))

Note how we grouped by both cond__ (the facet variable) and cats__ (the outcome categories), giving us the six separate set of connecting lines (2 facets * 3 outcome categories) we wanted. Also, we don’t need to include x and y variables inside aes, because fig1F (the ggplot produced by conditional_effects) already has the correct global x and y variables included in the plot mappings.