1
$\begingroup$

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.

Improve this question
$\endgroup$
5
  • $\begingroup$ 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}? $\endgroup$
    – dipetkov
    Sep 23 at 10:29
  • $\begingroup$ 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. $\endgroup$
    – Russ Lenth
    Sep 25 at 18:46
  • $\begingroup$ HPD stands for highest posterior density interval, it is not the confidence interval. $\endgroup$
    – Sextus Empiricus
    Sep 25 at 19:03
  • $\begingroup$ 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). $\endgroup$
    – Sextus Empiricus
    Sep 25 at 21:24
  • $\begingroup$ @SextusEmpiricus could you please submit a response with the R code that would achieve what you suggest? $\endgroup$
    – mat
    Sep 26 at 18:43

1 Answer 1

Reset to default
3
+50
$\begingroup$

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$.

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ 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? $\endgroup$
    – Sextus Empiricus
    Sep 25 at 19:05
  • $\begingroup$ @SextusEmpiricus I compared the emmeans intervals results to the one based on posterior draws and they do match $\endgroup$
    – mat
    Sep 29 at 8:56
  • $\begingroup$ @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) $\endgroup$
    – Russ Lenth
    Oct 16 at 17:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.