The key is to think about relative risks (or whatever is of interest) as an implication of the generative process you are trying to model instead of necessarily being a parameter in that model. In the simple case, you could estimate a Bernoulli or a binomial model with a good link function (not log
), obtain the posterior distribution of the probabilities under different scenarios, and then divide to obtain a relative risk ratio. In rstanarm (or brms) it would look like
library(rstanarm)
post <- stan_glm(y ~ treatment + x, data = dataset, family = binomial(link = "logit"), QR = TRUE)
nd <- dataset
nd$treatment <- 0
Pr_0 <- posterior_linpred(post, newdata = nd, transform = TRUE)
nd$treatment <- 1
Pr_1 <- posterior_linpred(post, newdata = nd, transform = TRUE)
RR <- Pr_1 / Pr_0
dim(RR)
where dim(RR)
will be number of simulations by number of observations and you can summarize that however you want.