Modelling and interpreting brms output

Question

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")

Many many thanks in advance!

Answer 1

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.