# Modelling and interpreting brms output

I do apologize in advance for this might be very basic questions. I am not really familiar with Bayesian statistics and too, unfortunately, this is the very first time I am analysing data in general.

My participants had to rate three players (A, B & C) on a 0 - 100% rating scale in steps of 10, so I got proportional data with several 0s and 1s. We had two treatmend conditions (Experimental & Control) and put participants into a group (low, medium, high). Participants did the experiment 2 times, at one session being in the control and at the other session being in the experimental condition but the group stayed the same. Each player was rated 2 times per session, so I got repeated measures.

I found out about the zero and one inflated beta regression in the brms package from Paul Bürkner and tried to fit a model, as the package is really nice and straight forward. I do not have any prior knowledge from other studies and anyway, even after reading several articles, vignettes and some sections of books I have no clue what I'd have to do to get the prior distribution for my data. Therefore I decided to stick with the default setting in the brms package, which uses weakly informative priors. My code looks like this:

rating = brm(Rating ~ when*Player*Condition +(1|ID), conf.df, family = zero_one_inflated_beta(), cores = 4,  save_ranef = F)


Now, I am not sure if I fit the models correctly or made some serious mistakes. I guess it might be the latter, as I get error messages when using coef(), plot(), icc(), pp_check, loo, WAIC and so on. Most of them say "Subscript out of bounds" or "attempt to set an attribute on NULL". plot() says "object 'Value' not found".

And even if it was correct I am stuck with interpreting the output. My estimates are given on the logit scale as I can see at the top of the summary() output. Regarding to this exp(estimate) should give me the OR relative to my reference category which here is [First rating - Player B - Control].

Family: zero_one_inflated_beta
Links: mu = logit; phi = identity; zoi = identity; coi = identity
Formula: Rating ~ when * Player * Condition + (1 | ID)
Data: conf.df (Number of observations: 1080)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000

Group-Level Effects:
~ID (Number of levels: 45)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)     0.30      0.04     0.23     0.39        798 1.00

Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept                                 b0 =  0.23      0.08     0.07     0.38       1203 1.00
whenlast                                  b1 = -0.02      0.09    -0.21     0.17       1276 1.00
PlayerA                                   b2 = -0.14      0.09    -0.32     0.05       1564 1.00
PlayerC                                   b3 = -0.04      0.09    -0.22     0.14       1558 1.00
ConditionExperimental                     b4 = -0.13      0.09    -0.31     0.06       1242 1.00
whenlast:PlayerA                          b5 =  0.59      0.13     0.32     0.86       1348 1.00
whenlast:PlayerC                          b6 = -0.23      0.13    -0.49     0.03       1518 1.00
whenlast:ConditionExperimental            b7 =  0.20      0.13    -0.08     0.46       1138 1.00
PlayerA:ConditionExperimental             b8 =  0.19      0.13    -0.07     0.44       1433 1.00
PlayerC:ConditionExperimental             b9 =  0.06      0.13    -0.21     0.31       1512 1.00
whenlast:PlayerA:ConditionExperimental   b10 = -0.29      0.19    -0.66     0.07       1275 1.00
whenlast:PlayerC:ConditionExperimental   b11 = -0.13      0.19    -0.49     0.25       1402 1.00

Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
phi    10.67      0.48     9.76    11.66       4714 1.00
zoi     0.10      0.01     0.08     0.12       5922 1.00
coi     0.80      0.04     0.72     0.87       5219 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).


For example for "whenlast:PlayerA" the effect is different from Null with an estimate of 0.59. But how do I mathematically get the percentage of change of this Player over time in regard to the reference? And too, if you have a look at "whenlast:PlayerC" the result tells me that the effect is not different from Null. In the graph however you can clearly see that, in the Control condition, the CIs of Player B at the first rating and Player C at the last rating do not overlap too much, indicating that there actually is a difference in rating - or am I totally wrong?

So, in short:

• Did I fit my model right? And if not, what is wrong?
• If it's correct, how do I interpret the results in terms of proportional change?

I really hope anyone can help me and explain all this in a very basic way.

Here is a reproducible example of my data:

structure(list(ID = structure(c(14L, 15L, 17L, 21L, 9L, 31L,
4L, 29L, 14L, 42L), .Label = c("1", "2", "3", "4", "5", "6",
"7", "8", "9", "10", "11", "12", "14", "15", "16", "17", "18",
"19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29",
"30", "31", "33", "34", "35", "36", "37", "38", "39", "40", "41",
"42", "43", "44", "45", "46", "47"), class = "factor"), Player = c("A",
"A", "B", "A", "A", "C", "B", "B", "C", "C"), Rating = c(100,
60, 60, 50, 80, 30, 60, 100, 50, 60), Condition = c("Experimental",
"Control", "Experimental", "Control", "Experimental", "Control",
"Experimental", "Control", "Control", "Experimental"), Group =
structure(c(2L,
3L, 3L, 2L, 1L, 3L, 1L, 1L, 2L, 3L), .Label = c("Low", "Medium",
"High"), class = "factor"), when = structure(c(2L, 1L, 1L, 2L,
2L, 1L, 2L, 2L, 1L, 1L), .Label = c("first", "last"), class = "factor")),
.Names = c("ID",
"Player", "Rating", "Condition", "Group", "when"), row.names = c(382L,
49L, 411L, 96L, 674L, 183L, 232L, 508L, 41L, 273L), class = "data.frame")


The model summary results you shared here via the summary() output refer to the logit-transfomed (estimated value of the) expected rating. In contrast, the plot shows the (estimated value of the) expected rating directly (not its logit transformation) as a function of the predictor variables included in the model. So it might be a bit difficult to map directly the summary() output and the plot.

You can use the summary() output to write down what the logit-transformed (estimated value of the) expected rating looks like:

log(mu_hat/(1-mu_hat)) = b0 + b1*whenlast +
b2*PlayerA + b3*PlayerC +
b4*ConditionExperimental +
b5*whenlast*PlayerA +
b6*whenlast*PlayerC +
b7*whenlast*ConditionExperimental +
b8*PlayerA*ConditionExperimental +
b9*PlayerC*ConditionExperimental +
b10*whenlast*PlayerA*ConditionExperimental +
b11*whenlast*PlayerC*ConditionExperimental


Since all the predictor variables in your model are dummy variables, you can use that to express the value of mu_hat for various combinations of values of the original variables when, Player and Condition. Here, the b coefficients represent the estimated model coefficients listed in the Estimate column of the summary() output. The dummy variables are defined as follows:

whenlast is equal to 1 if when = last and equal to 0 if when = first
PlayerA is equal to 1 for player A and equal to 0 for players B and C
PlayerC is equal to 1 for player C and equal to 0 for players A and B
ConditionExperimental is equal to 1 if Condition = Experimental and 0 if Condition = Control


The above equation shows that, if when = first, Player = B and Condition = Control, then whenlast = 0, PlayerA = 0, PlayerC = 0 and ConditionExperimental = 0, hence:

log(mu_hat/(1-mu_hat)) = b0


On the other hand, if when = first, Player = B and Condition = Experimental, then whenlast = 0, PlayerA = 0, PlayerC = 0 and ConditionExperimental = 1, hence:

log(mu_hat/(1-mu_hat)) = b0 + b4


So b4 represents the difference in the logit-transformed (estimated) expected values of rating between subjects in the the conditions who rated player A on the first occasion. Computing [1 - exp(b4)]*100% will give you the percent change in the value of mu_hat/(1-mu_hat) associated with changing Condition from Control to Experimental, but keeping when = first and Player = B. Note that mu_hat/(1-mu_hat) is just the ratio of the estimated expected rating to (1 - estimated expected rating), so not exactly odds in the traditional sense.

If you wanted to compute the (estimated) expected value of the rating for the combination when = first, Player = B and Condition = Control, you can just use the inverse logit transform:

mu_hat = exp(b0)/[1 + exp(b0)]


Similarly, the (estimated) expected value of the rating for the combination when = first, Player = B and Condition = Experimental is:

mu_hat = exp(b0 + b4)/[1 + exp(b0 + b4)]


You can proceed in a similar fashion for other combination of values of your three original predictor variables.

Not sure why your model does not include the group variable at all?

This link may come in handy: https://bayesat.github.io/lund2018/slides/andrey_anikin_slides.pdf.

• Thank you for this very clear and detailed answer! Just to be sure that I understood this right; so if I wanted to report the percentage of change in the rating for when = first, Player = B, condition = experimental (compared to placebo) than I’d have to calculate [1-exp(-0.13)]*100, correct? I could then report “Compared to the placebo condition, the rating of Player B at the first rating timepoint is most credibly 12% smaller in the experimental group, (mu = -0.13, 95% CI [-0.30 0.05])”? – joy2709 Mar 18 at 12:59
• I did not include the group variable because it decreased the n_eff of all other variables considerably, resulting in much broader CIs. When running more iterations the CIs get smaller, but the n_effs are still very low. Do you happen to know what I could do to increase the effective sample size? – joy2709 Mar 18 at 13:00
• For your first comment, you are correct since you are computing [1 - exp(b4)]*100% and b4 = -0.13; just make sure to multiply by 100% rather than just 100. I annotated your question to show the values of b0, b1, ..., b11 for ease of reference, which also makes it easy to see the value of b4. You would report the percentage change in the ratio of the ('expected' rating) to (1 - 'expected' rating) between the experimental group and the control group. That is, the percentage change in the value of mu_hat/(1 - mu_hat). – Isabella Ghement Mar 19 at 3:00
• Beta regression models don't really have an easy interpretation. See here: stats.stackexchange.com/questions/297659/…. In your case, it's actually easy to compute a point estimate of the difference in the expected rating between the experimental and placebo conditions for when = first and Player = B. Just take the difference between exp(b0 + b4)/[1 + exp(b0 + b4)] and exp(b0)/[1 + exp(b0)] and that will give you directly the most credible difference in expected ratings between the two conditions assuming when = first and Player = B. – Isabella Ghement Mar 19 at 3:09
• Thank you so much! This was so much help and I even though I actually did exactly what you explained above, I now finally understand the principles behind it and am sure that I did it right, so I marked your answer as my solution! For my thesis I think this analysis is sufficient, anyway, if anybody knows why adding the group variable makes my model fit worse any explanation is appreciated! – joy2709 Mar 19 at 10:02