# How to Compute Mean Ratios and Their 95% Confidence Intervals in a Bayesian Model

I am working on a Bayesian model using the brm function from the brms package in R, and I am interested in comparing mean responses of different groups. Specifically, I would like to calculate the mean ratios of responses, along with their 95% confidence intervals. I know how to calculate the odds ratios and mean differences, but I'm not sure how to go about calculating the mean ratios.

The Model and Data
Here's the model I used for demonstration purposes:

library(brms)

# Create a variable as a proportion of the max Sepal.Length with interval [0, 1[
iris$$Sepal.Length.Percent <- iris$$Sepal.Length / max(iris$Sepal.Length + 0.0000001) # Fit a Bayesian model with a Beta distribution mdl <- brm(bf(Sepal.Length.Percent ~ Species), family = Beta(), data = iris)  What I've Tried I used the emmeans package to get the expected marginal means and odds ratios, as well as 95% confidence intervals for those estimates. library(emmeans) > emmeans(mdl, pairwise ~ Species, type = 'response')$emmeans
Species    response lower.HPD upper.HPD
setosa        0.628     0.598     0.658
versicolor    0.744     0.718     0.773
virginica     0.867     0.848     0.888

$contrasts contrast odds.ratio lower.HPD upper.HPD setosa / versicolor 0.580 0.473 0.697 setosa / virginica 0.257 0.207 0.319 versicolor / virginica 0.444 0.346 0.547 > emmeans(mdl, pairwise ~ Species)$emmeans
Species    emmean lower.HPD upper.HPD
setosa      0.522     0.399     0.655
versicolor  1.066     0.917     1.206
virginica   1.878     1.716     2.066

$contrasts contrast estimate lower.HPD upper.HPD setosa - versicolor -0.545 -0.735 -0.349 setosa - virginica -1.360 -1.573 -1.143 versicolor - virginica -0.813 -1.036 -0.586  What I'm Looking For I want to find a way to calculate (or extract) the mean ratios, not odds ratios or mean differences, along with their 95% confidence intervals. Any guidance on how to achieve this would be much appreciated! Edit: Ideally, I would like to find the least cumbersome solution, for instance with a function from emmeans. • I am not sure what mean ratio is in this setting. Do you perhaps mean relative risk? (For groups A and B, this is the ratio of the probability {Y=1} in group A to the probability {Y=1} in group B which is defined for a binary Y.) In this case the outcome is a fraction, so E{Y|A} / E{Y|B}? Sep 23 at 10:29 • The idea is similar, but note in this question we are not estimating binomial probabilities, we are estimating mean via beta regression. So we want ratios of means, not ratios of probabilities. Sep 25 at 18:46 • HPD stands for highest posterior density interval, it is not the confidence interval. Sep 25 at 19:03 • Possibly a change of the link function from 'logit' to 'log', might establish what you want. The model will compute the contrasts on the log scale and convert this back to the exponential scale where the contrasts now represent a ratio. Much like Russ Lenth's answer, but the conversion back is done before computing the results (the mean and interval boundaries). Sep 25 at 21:24 • @SextusEmpiricus could you please submit a response with the R code that would achieve what you suggest? – mat Sep 26 at 18:43 ## 1 Answer A short answer is that you can add regrid = "log" to the emmeans call: > emmeans(mdl, pairwise ~ Species, regrid = "log", type = 'response')$emmeans
Species    response lower.HPD upper.HPD
setosa        0.627     0.596     0.657
versicolor    0.744     0.717     0.773
virginica     0.868     0.845     0.887

Point estimate displayed: median
Results are back-transformed from the log scale
HPD interval probability: 0.95

$contrasts contrast ratio lower.HPD upper.HPD setosa / versicolor 0.843 0.792 0.895 setosa / virginica 0.723 0.686 0.763 versicolor / virginica 0.857 0.817 0.892 Point estimate displayed: median Results are back-transformed from the log scale HPD interval probability: 0.95  This converts the means to the $$\log(\mu)$$ scale, so that when you summarize with type = "response", we get back-transformed estimates of $$\mu$$ and $$\exp(\log(\mu_i)-\log(\mu_j)) = \mu_i/\mu_j$$. • Is the HPD interval of the contrasts$log(\mu_i)-\log(\mu_j)$the same as the HPD interval of$\mu_i/\mu_j\$ by transforming the boundaries? Sep 25 at 19:05
• @SextusEmpiricus I compared the emmeans intervals results to the one based on posterior draws and they do match
– mat
Sep 29 at 8:56
• @SextusEmpiricus The intervals we obtain this way are HPD intervals for the difference of logs, back-transformed, as noted in the annotations. However, those intervals are not the HPD intervals for the ratios, because HPD is characterized by finding the shortest interval with that content. To get the actual HPD intervals for the ratios, obtain post<-exp(as.mcmc(EMM[[2]])) where EMM is the rersult of the above emmeans() call. Then call coda::HPDinterval(post) Oct 16 at 17:50