Bayesian Aggregation Simulatons

Author

Ilya Musabirov

Code

library(rstanarm)
library(bayestestR)
library(dplyr)
library(insight)

N1 = 600
baseline_probability1 = 0.1
main_effect1 = 0.3
f11 <- sample(c(T,F), N1, replace=TRUE)
f12 <- sample(c(T,F), N1, replace=TRUE)
y1 <- rbinom(N1, 1, baseline_probability1)
y1[f11==T] <- rbinom(length(y1[f11==T]), 
                   1, main_effect1)
y1[f12==T] <- rbinom(length(y1[f12==T]), 
                   1, main_effect1)

sim1 = data.frame(y = y1, f1 = f11, f2 = f12)

m1 <- stan_glm(y ~ f1+f2, data=sim1, family="binomial", refresh=-1)

Posterior for study 1 (weakly informative priors)

Code

describe_posterior(m1)  %>% print_md()

Summary of Posterior Distribution
Parameter	Median	95% CI	pd	ROPE	% in ROPE	Rhat	ESS
(Intercept)	-1.90	[-2.29, -1.54]	100%	[-0.18, 0.18]	0%	1.000	3375.00
f1TRUE	0.77	[ 0.40, 1.17]	100%	[-0.18, 0.18]	0%	1.000	3762.00
f2TRUE	0.74	[ 0.35, 1.13]	99.92%	[-0.18, 0.18]	0%	1.000	3586.00

Code

posteriors1 <- get_parameters(m1)

posteriors1 %>% 
  summarise(across(.fns=median)) %>% 
  as.numeric() -> posterior1_location

posteriors1 %>% 
  summarise(across(.fns=mad)) %>% 
  as.numeric() -> posterior1_scale

Code

N2 = 600
baseline_probability2 = 0.2
main_effect21 = 0.4
main_effect22 = 0.1

f21 <- sample(c(T,F), N2, replace=TRUE)
f22 <- sample(c(T,F), N2, replace=TRUE)
y2 <- rbinom(N2, 1, baseline_probability2)
y2[f21==T] <- rbinom(length(y2[f21==T]), 
                   1, main_effect21)
y2[f22==T] <- rbinom(length(y2[f22==T]), 
                   1, main_effect22)

sim2 = data.frame(y = y2, f1 = f21, f2 = f22)

sim3 = sim2 %>% sample_n(100)

m2.1 <- stan_glm(y ~ f1+f2, data = sim2,
                 family="binomial", 
               prior = normal(
                 location = posterior1_location[2:3], 
                 scale = posterior1_scale[2:3]))

m2.2 <- stan_glm(y ~ f1+f2, 
                 data = sim2, 
                 family="binomial")

m3.1 <- stan_glm(y ~ f1+f2, data = sim3,
                 family="binomial", 
               prior = normal(
                 location = posterior1_location[2:3], 
                 scale = posterior1_scale[2:3]))
m3.2 <- stan_glm(y ~ f1+f2, data = sim3,
                 family="binomial")

Study 1 (weakly informative priors) vs Study 2 (N = 600) (priors from study 1), Study 2 (WIP)

Code

describe_posterior(m1) %>% print_md()

Summary of Posterior Distribution
Parameter	Median	95% CI	pd	ROPE	% in ROPE	Rhat	ESS
(Intercept)	-1.90	[-2.29, -1.54]	100%	[-0.18, 0.18]	0%	1.000	3375.00
f1TRUE	0.77	[ 0.40, 1.17]	100%	[-0.18, 0.18]	0%	1.000	3762.00
f2TRUE	0.74	[ 0.35, 1.13]	99.92%	[-0.18, 0.18]	0%	1.000	3586.00

Code

describe_posterior(m2.1) %>% print_md()

Summary of Posterior Distribution
Parameter	Median	95% CI	pd	ROPE	% in ROPE	Rhat	ESS
(Intercept)	-1.32	[-1.59, -1.06]	100%	[-0.18, 0.18]	0%	0.999	3628.00
f1TRUE	0.71	[ 0.46, 0.97]	100%	[-0.18, 0.18]	0%	1.001	3947.00
f2TRUE	-0.24	[-0.50, 0.03]	96.00%	[-0.18, 0.18]	33.53%	1.000	3812.00

Code

describe_posterior(m2.2) %>% print_md()

Summary of Posterior Distribution
Parameter	Median	95% CI	pd	ROPE	% in ROPE	Rhat	ESS
(Intercept)	-0.91	[-1.23, -0.59]	100%	[-0.18, 0.18]	0%	1.000	3683.00
f1TRUE	0.67	[ 0.28, 1.07]	99.83%	[-0.18, 0.18]	0%	1.000	3037.00
f2TRUE	-1.21	[-1.62, -0.80]	100%	[-0.18, 0.18]	0%	1.001	3279.00

Study 1 (weakly informative priors) vs Study 2 (N = 100) (priors from study 1), Study 2 (WIP)

Code

describe_posterior(m1) %>% print_md()

Summary of Posterior Distribution
Parameter	Median	95% CI	pd	ROPE	% in ROPE	Rhat	ESS
(Intercept)	-1.90	[-2.29, -1.54]	100%	[-0.18, 0.18]	0%	1.000	3375.00
f1TRUE	0.77	[ 0.40, 1.17]	100%	[-0.18, 0.18]	0%	1.000	3762.00
f2TRUE	0.74	[ 0.35, 1.13]	99.92%	[-0.18, 0.18]	0%	1.000	3586.00

Code

describe_posterior(m3.1) %>% print_md()

Summary of Posterior Distribution
Parameter	Median	95% CI	pd	ROPE	% in ROPE	Rhat	ESS
(Intercept)	-1.71	[-2.27, -1.19]	100%	[-0.18, 0.18]	0%	1.000	3709.00
f1TRUE	0.78	[ 0.45, 1.14]	100%	[-0.18, 0.18]	0%	0.999	4141.00
f2TRUE	0.39	[ 0.06, 0.73]	98.70%	[-0.18, 0.18]	9.11%	0.999	3854.00

Code

describe_posterior(m3.2) %>% print_md()

Summary of Posterior Distribution
Parameter	Median	95% CI	pd	ROPE	% in ROPE	Rhat	ESS
(Intercept)	-0.92	[-1.86, -0.08]	98.50%	[-0.18, 0.18]	2.11%	1.001	3412.00
f1TRUE	0.86	[-0.18, 1.94]	94.83%	[-0.18, 0.18]	7.61%	1.001	3078.00
f2TRUE	-1.79	[-3.00, -0.77]	100%	[-0.18, 0.18]	0%	0.999	2762.00