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

## ----library------------------------------------------------------------------
library(causalsim)

## ----dgp----------------------------------------------------------------------
dgp <- causalsim_dgp(
  n = 500,
  n_confounders = 1,
  effect = 2,
  propensity = "moderate",
  baseline = "moderate"
)
dgp

## ----draw---------------------------------------------------------------------
dat <- causalsim_draw(dgp, seed = 1L)
head(dat)

## ----confounding-check--------------------------------------------------------
aggregate(W ~ A, data = dat, FUN = mean)

## ----estimators---------------------------------------------------------------
# Naive: regresses Y on A only, omits the confounder
naive_est <- function(data) {
  fit <- lm(Y ~ A, data = data)
  est <- coef(fit)[["A"]]
  se <- sqrt(vcov(fit)["A", "A"])
  c(estimate = est, ci_lower = est - 1.96 * se, ci_upper = est + 1.96 * se)
}

# OLS: adjusts for the observed confounder W
ols_est <- function(data) {
  fit <- lm(Y ~ A + W, data = data)
  est <- coef(fit)[["A"]]
  se <- sqrt(vcov(fit)["A", "A"])
  c(estimate = est, ci_lower = est - 1.96 * se, ci_upper = est + 1.96 * se)
}

## ----estimator-list, eval = FALSE---------------------------------------------
# ols_est <- function(data) {
#   fit <- lm(Y ~ A + W, data = data)
#   ci <- confint(fit)["A", ]
#   list(
#     estimate = coef(fit)["A"],
#     ci_lower = ci[1],
#     ci_upper = ci[2]
#   )
# }

## ----eval-naive---------------------------------------------------------------
eval_naive <- causalsim_eval(dgp, naive_est, reps = 300L, seed = 1L)
eval_naive

## ----eval-ols-----------------------------------------------------------------
eval_ols <- causalsim_eval(dgp, ols_est, reps = 300L, seed = 1L)
eval_ols

## ----summary-ols--------------------------------------------------------------
summary(eval_ols)

## ----plot-ols, fig.cap = "Distribution of OLS estimates over 300 replications. Solid line: true ATE. Dashed line: mean estimate."----
plot(eval_ols)

## ----grid-n-------------------------------------------------------------------
grid_n <- causalsim_eval_grid(
  dgp = dgp,
  estimator = ols_est,
  vary = list(n = c(100L, 250L, 500L, 1000L)),
  reps = 300L,
  metrics = c("bias", "rmse"),
  seed = 1L
)
grid_n

## ----grid-confounding---------------------------------------------------------
conf_levels <- list(propensity = c("low", "moderate", "high"))

grid_naive <- causalsim_eval_grid(dgp, naive_est,
                             vary = conf_levels,
                             reps = 300L,
                             metrics = "bias",
                             seed = 1L)

grid_ols <- causalsim_eval_grid(dgp, ols_est,
                           vary = conf_levels,
                           reps = 300L,
                           metrics = "bias",
                           seed = 1L)

comparison <- rbind(
  cbind(estimator = "naive", grid_naive$results),
  cbind(estimator = "ols",   grid_ols$results)
)
comparison <- comparison[order(comparison$propensity, comparison$estimator), ]
rownames(comparison) <- NULL
comparison

## ----het-dgp------------------------------------------------------------------
het_dgp <- causalsim_dgp(
  n = 4000,
  covariates = list(
    W = causalsim_covar("normal", role = "confounder"),
    V = causalsim_covar("binary", role = "effect_modifier", prob = 0.5)
  ),
  effect = function(V) 2 + 3 * V,   # effect is 2 when V = 0, 5 when V = 1
  propensity = function(W) plogis(0.5 * W),
  baseline = function(W) W
)
het_dgp

## ----het-truth----------------------------------------------------------------
d <- causalsim_draw(het_dgp, seed = 1L)
tapply(d$.tau, d$V, mean)   # 2 for V = 0, 5 for V = 1

## ----het-estimators-----------------------------------------------------------
overall <- lm(Y ~ A + W, data = d)         # assumes a constant effect
interact <- lm(Y ~ A * V + W, data = d)     # allows the effect to vary with V

c(
  average = coef(overall)[["A"]],
  subgroup_v0 = coef(interact)[["A"]],
  subgroup_v1 = coef(interact)[["A"]] + coef(interact)[["A:V"]]
)

## ----session-info-------------------------------------------------------------
sessionInfo()