Posterior distribution for difference of two propotions with brms in R

Modeling
,
1

I would like to produce a posterior distribution for the difference of two groups. I can do this in JAGS in R, but am hoping to move into this century and replicate this with Stan. I am using brms.

My two questions are:

  1. Have I correctly produced the posterior distribution for the effect size?
  2. How can I account for a prior effect size that could be negative? I imagine this requires using a non-beta prior, but am not sure what would be preferable.

The below code specifies a control and intervention group, with a control and intervention group event rate. It generates beta distributions for each group (step 1), prior distributions (step 2), builds a binomial model (step 3), and posteriors (step 4).

I am reasonably confident what I’ve done in the above steps is correct. I am a bit uncertain as to what the values generated by as_draws_df are - namely the b_Intercept_p vs. b_groupint_p. Further reading on these forums suggested subtracting b_groupint_p from b_Intercept_p; which does produce a plot that could plausibly be the posterior interval (labelled new_p when plotting the alt.graph object), but I am pretty uncertain on this point and would appreciate any clarification.

# Setup
## Packages
library(tidyverse)
library(ProbBayes)
library(brms)
library(tidybayes)

## Options
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())



# ***
# Setup simulation
# ***

# Trial data
con.n = 150
con.event = 50
int.n = 150
int.event = 30


# Distributions
## x1 is the expected control event rate at probability p.x1
## x2 is the control event rate at probability p.x2
x1 = 0.4
p.x1 = 0.5
x2 = 0.6
p.x2 = 0.9

## As above, but for the difference between intervention and control group
y1 = 0.05
p.y1 = 0.5
y2 = 0.1
p.y2 = 0.9


# Model setup
## For brm()
n.chains = 3
n.iter = 4000
n.warmup = 500

## For beta.select(); based on model so the number of observations are equal
n.draws = n.iter - n.warmup



# ***
# Simulation starts here
# ***



# Define Functions
## Log transformation
fun_log.trans = function(x) {
  log.trans = log(x / (1 - x))
  log.trans
}


## Inverse log transformation
fun_invlog.trans = function(x) {
  invlog.trans = exp(x) / (1 + exp(x))
  invlog.trans
}



# Run things
# 0. Put data into dataframe
data = data.frame(group = c("int"    , "con"),
                  n =     c(int.n    , con.n),
                  event = c(int.event, con.event))


# 1. Generate prior distributions
## Beta prior for the control group event rate
beta0.val = beta.select(list(x = x1, p = p.x1),
                        list(x = x2, p = p.x2))
p0.sim = rbeta(n.draws, beta0.val[1], beta0.val[2])
### Log transform it
theta0.sim = fun_log.trans(p0.sim)


## Prior distribution for the difference in logit-p for each group
beta1.val = beta.select(list(x = y1, p = p.y1),
                        list(x = y2, p = p.y2))
p1.sim = rbeta(n.draws, beta1.val[1], beta1.val[2])
### Log transform
theta1.sim = fun_log.trans(p1.sim)


## 2. Create a matrix to store priors
priors = get_prior(family = binomial,
                   event | trials(n) ~ group,
                   data = data)

### Get characteristics of the normal distribution for the priors generated in step 1
theta0.val = c(mean(theta0.sim), sd(theta0.sim))
theta1.val = c(mean(theta1.sim), sd(theta1.sim))

### Input the these characteristics into the prior matrix generated at the start of step 2
priors$prior[3] = paste("normal(", theta0.val[1], ", ", theta0.val[2], ")", sep = "")
priors$prior[2] = paste("normal(", theta1.val[1], ", ", theta1.val[2], ")", sep = "")



# 3. Generate model with stan
model = brm(family = binomial,
            event | trials(n) ~ group,
            data = data,
            prior = priors,
            iter = n.iter,
            warmup = n.warmup,
            chains = n.chains,
            refresh = 0)



## Plot model
plot(model)
print(model)

# 4. Generate posteriors
posteriors = as_draws_df(model)

## Inverse log transform function on theta to get p again
posteriors = posteriors %>%
  mutate(across(starts_with("b_"), .f = fun_invlog.trans, .names = "{.col}_p")) %>%
  rename_with(~paste0(., "_theta"), .cols = starts_with("b_") & !ends_with("_p"))

## 95% posterior interval estimate
quantile(posteriors$b_Intercept_p, c(0.025, 0.975))


# 5. Plot posterior densities
## Take the posterior interval data and bind it with the priors
## Ideally, n.draws = number of iterations - warmup
graph = posteriors %>%
  select(ends_with("_p")) %>%
  cbind(p0.sim, p1.sim) %>%
  pivot_longer(cols = everything(),
               names_to = "distribution",
               values_to = "value")

alt.graph = posteriors %>%
  mutate(new_p = b_Intercept_p - b_groupint_p) %>%
  select(new_p) %>%
  cbind(p0.sim, p1.sim)

quantile(alt.graph$new_p, c(0.025, 0.975))

alt.graph = alt.graph %>%
  pivot_longer(cols = everything(),
               names_to = "distribution",
               values_to = "value")

## Plot it
alt.graph %>%
  ggplot(aes(x = value)) +
  geom_density(aes(fill = distribution),
               alpha = 0.5) +
  theme_light()