ploting_conjoint_interactions.R

rm(list = ls())

library(tidyverse)
library(estimatr)

N <- 1000

dat <-
  data.frame(
    Y = rnorm(N),
    X1 = rbinom(N, 1, .5),
    X2 = rbinom(N, 1, .5),
    X3 = rbinom(N, 1, .5),
    X4 = rbinom(N, 1, .5),
    X5 = rbinom(N, 1, .5),
    Z = rbinom(N, 1, .5) # thing we want to split by
  )


# fit three regressions
fit_Z0 <- lm_robust(Y ~ X1 + X2 + X3 + X4 + X5, data = filter(dat, Z == 0))
fit_Z1 <- lm_robust(Y ~ X1 + X2 + X3 + X4 + X5, data = filter(dat, Z == 1))
fit_int <- lm_robust(Y ~ (X1 + X2 + X3 + X4 + X5)*Z, data = dat)

# prepare output for ggplot
fit_Z0_df <- fit_Z0 %>% tidy() %>% filter(term != "(Intercept)") %>% mutate(model = "Z0")
fit_Z1_df <- fit_Z1 %>% tidy() %>% filter(term != "(Intercept)") %>% mutate(model = "Z1")
fit_int_df <- fit_int %>% tidy() %>% filter(grepl(pattern = ":", x = term)) %>% mutate(model = "Difference")
fit_int_df$term <- fit_Z0_df$term # fix up the label

gg_df <- bind_rows(fit_Z0_df, fit_Z1_df, fit_int_df) %>%
  mutate(model = factor(model, levels = c("Z0", "Z1", "Difference")))

# plot
ggplot(data = gg_df, aes(estimate, term)) +
  geom_point() +
  geom_errorbarh(aes(xmin = ci.lower, xmax = ci.upper), height = 0) +
  facet_wrap( ~ model) +
  geom_vline(xintercept = 0) +
  theme_bw()