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

- Operating System: Mac OS Sierra
- brms Version: 2.4.0 (I think)

Hi!

I ran an experiment where I obtained perceptual data on listeners’s sensitivity for recognizing a token in two different levels of noise (low and high). I have two types of tokens from three speakers ( for example “pat”/“bat”. In total I therefore have 12 datapoints per listener (2 levels noise x two tokens x 3 speakers.

Let’s say the following model best fit my data:

model7B.model <- brm(sensitivity_scores ~ noise_level * token + speaker + speaker:token + (1|listener), data=mydata, iter = 20000, control = list(adapt_delta = 0.99), init = 0)

It would be great if it were possible to display the marginal effects plots for the interaction between conditions, for example noise_level and token. I’ve experimented with using the function “marginal_effects”, but as a result of my conditions being categorical, marginal_effects displays only the reference category each time, for example

marginal_effects(model7B.model, effects = noise_level:token)

displays the interaction only for the first speaker.

Instead, I’d like to have a figure that displays what the noise_level by token interaction looks like if speakers were lumped/pooled together, rather than just for one category of data (speaker 1 / speaker 2 etc.)

I had two ideas on how I might achieve this:

**Idea 1**

I simply rerun the Bayesian statistics with this model:

model4B.model <- brm(sensitivity_scores ~ noise_level * token + (1|listener), data=mydata, iter = 20000, control = list(adapt_delta = 0.99), init = 0)

And then use

marginal_effects(model4B.model, effects = “noise_level:token”))

To get the interaction plot of noise_level and vowel if speaker was not a factor. I can then do the same for the there interaction.

**Idea 2**

I extracted the posterior draws using:

intercept_draws = model7B.model[[“fit”]]@sim[["samples”]][[1]][[“b_Intercept”]][(4001:20000)]

noise_level_hi_draws = model7B.model[[“fit”]]@sim[["samples”]][[1]][[“b_noise_level”]][(4001:20000)]

token2_draws = model7B.model[[“fit”]]@sim[["samples”]][[1]][[“b_token”]][(4001:20000)]

speaker2_draws = model7B.model[[“fit”]]@sim[["samples”]][[1]][[“b_speaker1”]][(4001:20000)]

speaker3_draws = model7B.model[[“fit”]]@sim[["samples”]][[1]][[“b_speaker2”]][(4001:20000)]

token2_noise_level_hi_draws = c(model7B.model[[“fit”]]@sim[[“samples”]][[1]][[“b_dBlevelhi:voweliy”]][(4001:20000)])

token2_speaker2_draws = c(model7B.model[[“fit”]]@sim[[“samples”]][[1]][[“b_voweliy:speakerm30”]][(4001:20000)])

token2speaker3_draws = c(model7B.model[[“fit”]]@sim[[“samples”]][[1]][[“b_voweliy:speakerm44”]][(4001:20000)])

(In these particular lines I used 4001:20000 to exclude the 20% warmup values. This line only extracts the posterior draw values for the first chain for the purposes of the question.)

I then pooled together the speakers by adding up the posterior draws of each speaker for a particular condition, say low noise level and first token, and dividing by the number of groups I added (so, 3 for the three speakers). See below:

#The next line of code adds together:

#the draws for the the speaker 1, low noise level and token 1 (intercept draws);

#the draws for the the speaker 2, low noise level and token 1 (intercept draws+speaker2_draws);

#the draws for the the speaker 3, low noise level and token 1 (intercept draws+speaker2_draws);

#Then divide by 3, for the three groups added together.

#subsequent 3 lines essentially do the same for different conditionspoolspeakers_lo_noise_t1 = ((intercept_draws) + (intercept_draws + speaker2_draws) + (thedrawsIntercept + speaker3_draws))/3

poolspeak_lo_noise_t2 = ((intercept_draws + thedrawsvoweliy) + (intercept_draws + token2_draws + speaker2_draws + token2_speaker2_draws) + thedrawsIntercept + token2_draws + speaker3_draws + token2speaker3_draws)/3

poolspeak_hi_noise_t1 = ((intercept_draws + noise_level_hi_draws) + (intercept_draws + noise_level_hi_draws + speaker2_draws) + (intercept_draws + noise_level_hi_draws + speaker3_draws))/3

poolspeak_hi_noise_t2 = ((intercept_draws + noise_level_hi_draws + token2_draws + token2_noise_level_hi_draws) + (intercept_draws + noise_level_hi_draws + token2_draws + noise_level_hi_draws + speaker2_draws + token2_speaker2_draws) + (intercept_draws + noise_level_hi_draws + token2_draws + noise_level_hi_draws + speaker3_draws + token2speaker3_draws))/3

I plot the average of the pooled categories and the highest density intervals as error bars.

I’m pretty new to Bayesian statistics, so I’m still trying to wrap my head around everything, but for idea 1, I worry that a different model will experience different “shrinkage” (model4B.model compared to model7B.model).

For idea 2, I’m assuming that the new pooled posterior draws I “create” is a reasonable representation of what the Bayesian model would draw if speakers were not a factor, but then this seems to be essentially what idea 1 will do, and perhaps “shrinkage” will again be an issue?

Essentially I’d like to know if either of these ideas actually yield a valid marginal effects plot for my categorical data that are conditioned on some “average” rather than the reference category.

Any information or suggestions on how to get the correct plot would be greatly appreciated!