How to get the output from brms binomial regression model in risk ratio rather than log odds ratio scale directly?

Interfaces brms
1

I would like to directly estimate the adjusted risk ratio from my brms binomial regression as recommended here.

Alternatively, I would appreciate if anyone could help implement this formula by McNutt in brms in order to get the estimates in risk ratio scale rather than log odds ratio:

RRi=P(Y|Ea, x2, …, xki)/P(Y|Eb, x2i, …, xki)

where:

Y is the outcome factor of interest (dependent variable), E is the exposure of interest, and x 2 , …, x k are confounders.

The reason for that is that as pointed out here, for common outcomes as it is in my case, “the odds ratio always overstates the relative risk, sometimes dramatically” and the formula proposed by Zhang and McNutt to convert Odds ratio into risk ratio as implemented here yield biased results.

Below are sample data and a binomial regression model for which I would like to get the

Estimate Est.Error l-95%CI u-95%CI 95%HDI

in RR scale when I do

summary(model)

2

Try out

pl <- posterior_linpred(
  mod, newdata = data.frame(n = 1, treat = c(1, 0), ...),
  transform = TRUE
)
posterior_summary(pl[, 1] / pl[, 2])

I don’t have time to test this myself right now, but I hope the doc of posterior_linpred can help you further.

1 Like
3

Thanks @paul.buerkner. Just tried out.

But it gives the error:

Error in get(g, data) : object ‘group’ not found

Thanks in advance for this.

4

Please read the documentation of posterior_linpred carefully. In your case, you may have to set re_formula = NA if you want to exclude varying effects from the prediction.

5

Thanks, @paul.buerkner. The documentation was helpful.

I ran the code

lp <- posterior_linpred(mod, newdata = data.frame(n = 1, treat = c(1, 0), group = c(1, 0), c2=c(0, 1), c1=c(0, 1)), transform = TRUE, re_formula = NA)
posterior_summary(lp[, 1] / lp[, 2])

and got the results:

     Estimate Est.Error         Q2.5    Q97.5
[1,] 0.1530741 0.2798745 0.0002440068 0.987347

Please kindly help me get the results for the other predictors. I mean, how should I proceed to get results that are in a format and order like those in log odds ratio scale from summary(mod), that is, one row for each variable?

Thanks in advance for your help.

6

I think there is a typo. It should be posterior_summary(lp[, 1] / lp[, 2]).

Since you have interactions, you have to think carefully to which single value (per predictor) you fix each predictor as this affects the efficiency of the treatment as predicted by your model.

The same problem would occur when just looking at odds-ratios, but you perhaps didn’t realize this complexity because you were just looking at the main effect of treat. If you want to get the predictions for the main effect, set other predictor variables to 0.

In any case, you won’t just get the same output structure and format as in the summary output. Instead, you have to think carefully what exactly you want to compute and then use posterior_linpred for it.