I’ll answer this question in two parts. In the first part, I’ll answer your basic question directly. In the second part, I’ll suggest an alternative workflow.
Part 1: How to do this
I’m not aware of a simple way to convert beta priors to normal priors on the log-odds scale (though see this link for some ggdist wizardry). But a simple iterative way to get there is to:
a) plot the beta priors of interest and compute their basic descriptive statistics,
b) simulate a large number of draws from some normal distribution,
c) convert those normal draws to the log-odds scale with the plogis() function,
d) plot those transformed draws,
e) take the basic descriptive statistics from those transformed draws,
f) compare the results of your simulation with their target beta distribution, and finally
h) iterate by adjusting the mean and/or sd arguments in rnorm() until the results are what you want.
Here’s what this looks like in practice.
a) plot the beta priors and compute sample stats.
Taking your beta(25, 122) prior as an example, here’s how to plot the prior with helper functions from ggdist.
# load packages
library(tidyverse)
library(brms)
library(ggdist)
# plot
prior(beta(25, 122)) %>%
parse_dist() %>%
ggplot(aes(xdist = .dist_obj)) +
stat_halfeye(point_interval = mean_qi, .width = .95) +
scale_y_continuous(NULL, breaks = NULL) +
labs(title = "Beta(25, 122)")
Looks great, doesn’t it? I’m a huge fan of that ggdist and its handy parse_dist() and stat_halfeye() functions. Here’s how to compute the population mean and variance of that beta prior using its two parameters.
# beta parameters
a <- 25
b <- 122
# mean
a / (a + b)
# variance
a * b / ((a + b)^2 * (a + b + 1))
[1] 0.170068
[1] 0.0009536817
That overall shape and those descriptive statistics are our targets. [You could focus on other sample statistics, like quantiles or whatever. Means and variances are just a fine place to start.]
b) simulate a large number of normal draws.
The reason we’re simulating draws from the normal distribution is because it’s a good default distribution for \beta priors in logistic regression models. You could use some other distribution, like a small-\nu Student-t distribution, if you wanted. But anyway, pick some distribution and start simulating. Since we’ve chosen the normal distribution, we simulate with rnorm(). Note the values I’ve set for the mean and sd arguments. Pick different ones and see what happens
# you want some large number of draws
n_draw <- 1e6
# set your seed to make it reproducible
set.seed(1)
# simulate
my_simulation <- tibble(n = rnorm(n = n_draw, mean = -1.5, sd = 0.25))
# what have we done?
glimpse(my_simulation)
Rows: 1,000,000
Columns: 1
$ n <dbl> -1.656613, -1.454089, -1.708907, -1.101180, -1.417623, -1.705117, -1.378143, -1.315419, -1.356055, -1.…
Now we have 1{,}000{,}000 draws to play with.
c) convert those draws with plogis() .
Next we use the plogis() function to put those draws on the log-odds scale. We’ll save the results in a new column called p.
my_simulation <- my_simulation %>%
mutate(p = plogis(n))
# what have we done?
glimpse(my_simulation)
Rows: 1,000,000
Columns: 2
$ n <dbl> -1.656613, -1.454089, -1.708907, -1.101180, -1.417623, -1.705117, -1.378143, -1.315419, -1.356055, -1.…
$ p <dbl> 0.1602171, 0.1893730, 0.1533055, 0.2495189, 0.1950345, 0.1537981, 0.2013074, 0.2115815, 0.2048823, 0.1…
The n column has the values on their original scale, and the p column has the values after transforming them with plogis().
d) plot those transformed draws.
Here’s the plot.
my_simulation %>%
ggplot(aes(x = p)) +
stat_halfeye(point_interval = mean_qi, .width = .95) +
scale_y_continuous(NULL, breaks = NULL) +
coord_cartesian(xlim = c(0, 1)) +
labs(title = expression(logit^{-1}*"[Normal(-1.5, 0.25)]"))
e) more descriptive statistics.
Now compute the sample mean and variance for our transformed simulated draws.
my_simulation %>%
summarise(mean = mean(p), # the target mean is 0.170068
variance = var(p)) # the target variance is 0.0009536817
# A tibble: 1 × 2
mean variance
<dbl> <dbl>
1 0.185 0.00142
f) compare the results.
If you compare the plot and sample statistics of our simulation, you’ll see they’re pretty close to their beta targets, not perfectly there. If this is close enough for you, move forward with the model and your new priors. But…
h) iterate.
If you’d like to get a little closer to the target beta shape and descriptive statistics, try changing the mean and/or sd arguments within rnorm() way back in step b) and just keep iterating until you’re satisfied. Note, there’s no guarantee you’re going to find the normal distribution that perfectly replicates the original beta distribution after that plogis() transformation. So it goes…
Part 2: Should one do this?
Maybe. I do stuff like this sometimes, even as recently as last week. However, I’m not so sure this is actually the best method for your use case. If your goal is just to learn how to work with brm() models, this is a great practice. But if your goal is to analyze real-world data with a scientifically justified model for peer review, I would not try to fit the model two ways, one with the identity link and the other with the conventional logit link. Rather, I would just fit the model with the priors you believe in, which appears to be the one with the identity link. If your colleagues prefer odds ratios, you can just work with the posterior draws to transform the posterior probabilities into an odds ratio metric.
Let’s say your identity-link model is still called m3. Here’s how to pull the posterior draws with as_draws_df(). We’ll save the results as a data frame named draws.
draws <- as_draws_df(m3)
# what?
head(draws)
# A draws_df: 6 iterations, 1 chains, and 4 variables
b_grp20 b_grp21 lprior lp__
1 0.16 0.33 4.8 -2.1
2 0.16 0.22 1.3 -2.9
3 0.16 0.35 4.6 -3.1
4 0.15 0.36 4.3 -3.9
5 0.19 0.31 4.5 -2.2
6 0.16 0.34 4.7 -2.6
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}
The head() function returned the first 6 rows, but the full draws data frame has 4,000 rows by default. The posterior draws for the probabilities of your two groups are in the b_grp20 and b_grp21 columns. In the next step, we’ll give them simpler names, p0 and p1, drop the unneeded columns, and hand compute the odds ratio.
draws <- draws %>%
# rename and drop the unneeded columns
transmute(p0 = b_grp20,
p1 = b_grp21) %>%
# compute the OR
mutate(or = (p1 / (1 - p1)) / (p0 / (1 - p0)))
# what have we done?
head(draws)
# A tibble: 6 × 3
p0 p1 or
<dbl> <dbl> <dbl>
1 0.163 0.326 2.48
2 0.160 0.218 1.46
3 0.164 0.346 2.68
4 0.151 0.356 3.09
5 0.195 0.314 1.89
6 0.161 0.338 2.66
Now you can summarize or display your odds ratio as desired. For example, here it is in a plot.
draws %>%
ggplot(aes(x = or)) +
stat_halfeye(point_interval = mean_qi, .width = .95) +
scale_y_continuous(NULL, breaks = NULL) +
xlab("My odds ratio")
No second model needed. All you need to know is how to convert your posterior probabilities to the metric your colleagues prefer.
