library(rstanarm)
library(bayestestR)library(dplyr)
library(tidyr)
library(distributional)n=100
simdata = tibble(contextual=numeric(n), arm=numeric(n), outcome=double(n))
simdata <- simdata %>%
mutate(contextual = c(rep(0, n/2), rep(1,n/2)),
arm = rep(c(0,1,0,1), n/4)) %>%
mutate(prob = case_when(
contextual == 1 & arm == 1 ~ 0.6,
contextual == 0 & arm == 1 ~ 0.55,
contextual == 1 & arm == 0 ~ 0.5,
contextual == 0 & arm == 0 ~ 0.45
)) %>%
mutate(dist = dist_bernoulli(prob)) %>%
mutate(outcome = as.numeric(generate(dist,1)))
# Check distribution
simdata %>%
group_by(contextual, arm) %>%
summarise(
n = n(),
mean_outcome = mean(outcome),
.groups = 'drop'
) %>% gt::gt()| contextual | arm | n | mean_outcome |
|---|---|---|---|
| 0 | 0 | 25 | 0.40 |
| 0 | 1 | 25 | 0.64 |
| 1 | 0 | 25 | 0.64 |
| 1 | 1 | 25 | 0.68 |
m1 <- stan_glm(outcome~ factor(arm) + factor(contextual), data=simdata,refresh=0)
describe_posterior(m1,ci=0.9) %>% print_md()| Parameter | Median | 90% CI | pd | ROPE | % in ROPE | Rhat | ESS |
|---|---|---|---|---|---|---|---|
| (Intercept) | 0.45 | [ 0.31, 0.59] | 100% | [-0.05, 0.05] | 0% | 1.000 | 5420.00 |
| factor(arm)1 | 0.14 | [-0.02, 0.30] | 93.38% | [-0.05, 0.05] | 14.68% | 0.999 | 5235.00 |
| factor(contextual)1 | 0.14 | [-0.03, 0.30] | 91.47% | [-0.05, 0.05] | 17.47% | 1.000 | 5475.00 |