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.91	[-2.29, -1.53]	100%	[-0.18, 0.18]	0%	1.000	3134.00
f1TRUE	0.73	[ 0.33, 1.14]	99.98%	[-0.18, 0.18]	0%	1.000	3234.00
f2TRUE	0.61	[ 0.23, 1.00]	99.98%	[-0.18, 0.18]	0%	1.000	3654.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.91	[-2.29, -1.53]	100%	[-0.18, 0.18]	0%	1.000	3134.00
f1TRUE	0.73	[ 0.33, 1.14]	99.98%	[-0.18, 0.18]	0%	1.000	3234.00
f2TRUE	0.61	[ 0.23, 1.00]	99.98%	[-0.18, 0.18]	0%	1.000	3654.00

Code

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

Summary of Posterior Distribution
Parameter	Median	95% CI	pd	ROPE	% in ROPE	Rhat	ESS
(Intercept)	-2.09	[-2.43, -1.76]	100%	[-0.18, 0.18]	0%	0.999	3682.00
f1TRUE	0.76	[ 0.45, 1.08]	100%	[-0.18, 0.18]	0%	1.000	3247.00
f2TRUE	-0.16	[-0.47, 0.15]	86.08%	[-0.18, 0.18]	55.68%	1.001	3861.00

Code

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

Summary of Posterior Distribution
Parameter	Median	95% CI	pd	ROPE	% in ROPE	Rhat	ESS
(Intercept)	-1.71	[-2.15, -1.31]	100%	[-0.18, 0.18]	0%	0.999	2912.00
f1TRUE	0.81	[ 0.35, 1.31]	99.95%	[-0.18, 0.18]	0%	1.000	2918.00
f2TRUE	-1.23	[-1.72, -0.76]	100%	[-0.18, 0.18]	0%	1.000	2873.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.91	[-2.29, -1.53]	100%	[-0.18, 0.18]	0%	1.000	3134.00
f1TRUE	0.73	[ 0.33, 1.14]	99.98%	[-0.18, 0.18]	0%	1.000	3234.00
f2TRUE	0.61	[ 0.23, 1.00]	99.98%	[-0.18, 0.18]	0%	1.000	3654.00

Code

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

Summary of Posterior Distribution
Parameter	Median	95% CI	pd	ROPE	% in ROPE	Rhat	ESS
(Intercept)	-2.09	[-2.72, -1.50]	100%	[-0.18, 0.18]	0%	1.000	3707.00
f1TRUE	0.76	[ 0.39, 1.12]	99.98%	[-0.18, 0.18]	0%	1.000	3889.00
f2TRUE	0.34	[-0.02, 0.71]	96.65%	[-0.18, 0.18]	17.47%	1.000	3457.00

Code

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

Summary of Posterior Distribution
Parameter	Median	95% CI	pd	ROPE	% in ROPE	Rhat	ESS
(Intercept)	-1.63	[-2.77, -0.69]	100%	[-0.18, 0.18]	0%	1.000	2544.00
f1TRUE	1.11	[ 0.02, 2.40]	97.78%	[-0.18, 0.18]	3.37%	1.000	2633.00
f2TRUE	-1.55	[-2.89, -0.40]	99.62%	[-0.18, 0.18]	0%	1.000	2614.00