## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  cache = TRUE,
  comment = "#>",
  fig.width=7, 
  fig.height=4
)

## ----message = FALSE, warning = FALSE-----------------------------------------
library(dplyr)
library(ggplot2)
library(lmtest)
library(sandwich)
library(prepost)

set.seed(02134)
theme_set(theme_classic())


data(delponte)

delponte |> select(t_commonality, itaid_bin, angry_bin) |> sample_n(5)

## -----------------------------------------------------------------------------
lm(angry_bin ~ t_commonality * itaid_bin, data = delponte) |> 
  coeftest(vcov = vcovHC)

## -----------------------------------------------------------------------------
lm(itaid_bin ~ t_commonality, data = delponte) |> 
  coeftest(vcov = vcovHC)

## -----------------------------------------------------------------------------
delponte_subset <- delponte |>
  filter(sopscale >= quantile(sopscale, 0.25, na.rm = TRUE),
         Corriere == 1)

lm(itaid_bin ~ t_commonality, data = delponte_subset) |>
  coeftest(vcov = vcovHC)

## -----------------------------------------------------------------------------
bounds_1 <- post_bounds(formula = angry_bin ~ t_commonality,
                        data = delponte_subset,
                        moderator = ~ itaid_bin)

print(bounds_1)

## -----------------------------------------------------------------------------
bounds_2 <- post_bounds(
  formula = angry_bin ~ t_commonality,
  data = delponte_subset,
  moderator = ~ itaid_bin,
  moderator_mono = 1,
  stable_mod = TRUE
)

print(bounds_2)

## -----------------------------------------------------------------------------
gibbs_1 <- prepost_gibbs_nocovar(
  formula = angry_bin ~ t_commonality,
  data = delponte_subset,
  prepost = ~1,
  moderator = ~ itaid_bin,
  iter = 250
)

## -----------------------------------------------------------------------------
gibbs_2 <- prepost_gibbs(
  formula = angry_bin ~ t_commonality,
  data = delponte_subset,
  prepost = ~1,
  moderator = ~ itaid_bin,
  covariates = ~ north + satisf,
  iter = 250
)

## -----------------------------------------------------------------------------
delta_list <- list(gibbs_1$delta.1,
                   gibbs_1$delta.2,
                   gibbs_2$delta.1,
                   gibbs_2$delta.2)

get_ci <- function(x){
  tribble(~mean, ~lowerCI, ~upperCI,
          mean(x), quantile(x, 0.025), quantile(x, 0.975))
}

lapply(delta_list, get_ci) |>
  bind_rows() |>
  mutate(model = c(
    "In-sample delta, w/o covars",
    "Population delta, w/o covars",
    "In-sample delta, with covars",
    "Population delta, with covars"
  )) |>
  as_tibble() |>
  ggplot(aes(
    x = mean,
    y = model,
    xmin = lowerCI,
    xmax = upperCI
  )) +
  geom_point() +
  geom_errorbar(width = 0.2) +
  geom_vline(xintercept = 0, lty = "dotted") +
  labs(x = "", 
       y = "", 
       title = "Bayesian credible intervals")