set.seed(9)
#number of sims
N = 1000
# Generate data
U <- rnorm(N) # Unobserved confounder
X <- rnorm(N, mean = 0.5 * U) # Treatment affected by U
Y <- rnorm(N, mean = 0.8 * U) # Outcome affected by U
Z <- rnorm(N, mean = 0.6 * U) # Observed variable that captures U
d <- data.frame(X, Y, Z)
# Fit models
flist1 <- alist(
Y ~ dnorm(mu, sigma),
mu <- a + bX*X,
a ~ dnorm(0, .5),
bX ~ dnorm(0, .25),
sigma ~ dexp(1)
)
m1 <- brm(
data=d,
family=gaussian,
Y ~ 1 + X,
prior = c( prior(normal(0, .50), class=Intercept),
prior(normal(0, .25), class=b),
prior(exponential(1), class=sigma)),
iter=2000, warmup=1000, seed=3, chains=1
)
posterior_summary(m1)
m2 <- brm(
data=d,
family=gaussian,
Y ~ 1 + X + Z,
prior = c( prior(normal(0, .50), class=Intercept),
prior(normal(0, .25), class=b),
prior(exponential(1), class=sigma)),
iter=2000, warmup=1000, seed=3, chains=1
)
posterior_summary(m2)
post.1 <- as_draws_df(m1)
post.2 <- as_draws_df(m2)
results_df = data.frame(naive = post.1$b_X,
adjusted = post.2$b_X)
results_df %>%
pivot_longer(everything()) %>%
ggplot(aes(x = value, fill = name)) +
geom_density(alpha = .5) +
geom_vline(aes(xintercept = 0), linetype = "dashed")