---
title: "A Complete Simulation Study with causalsim"
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{A Complete Simulation Study with causalsim}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

```{css, echo = FALSE}
.callout {
  background-color: #f0f4f8;
  border-left: 4px solid #3a7ebf;
  padding: 0.75em 1em;
  margin: 1.25em 0;
  border-radius: 0 4px 4px 0;
}

.callout-warn {
  background-color: #fff8e1;
  border-left: 4px solid #f0a500;
  padding: 0.75em 1em;
  margin: 1.25em 0;
  border-radius: 0 4px 4px 0;
}

table {
  width: 100%;
  border-collapse: collapse;
  font-size: 0.9em;
  margin: 1em 0;
}

th {
  background-color: #3a7ebf;
  color: white;
  padding: 0.5em 0.75em;
  text-align: left;
}

td {
  padding: 0.4em 0.75em;
  border-bottom: 1px solid #dde;
}

tr:nth-child(even) td {
  background-color: #f5f7fa;
}
```

## Overview

Evaluating causal estimators requires knowing the answer. `causalsim` gives you
that by making the data-generating process explicit: you specify the structural
model, the package simulates data from it, and you measure how well any estimator
recovers the truth you specified.

If you only need a single dataset with known ground truth, `causalsim()` does it
in one call:

```r
data <- causalsim(n = 500, n_confounders = 1, effect = 2, seed = 1L)
head(data)
```

The columns `.tau` and `.p` carry the ground truth for whatever analysis you run
on the data. This vignette goes further: defining a DGP explicitly, evaluating
estimator performance over many replications, and sweeping across a parameter
grid.

### What this vignette covers

| Step | Function | What it does |
|---|---|---|
| 0 | `causalsim()` | Simulate one dataset in a single call |
| 1 | `causalsim_dgp()` | Define the structural model and true ATE |
| 2 | `causalsim_draw()` | Simulate one dataset and inspect it |
| 3 | `causalsim_eval()` | Measure estimator performance over many replications |
| 4 | `causalsim_eval_grid()` | Sweep over sample sizes and confounding levels |

```{r library}
library(causalsim)
```

---

## Step 1: Define the DGP

The structural model for this study is:

$$
W \sim N(0,1), \quad
A \mid W \sim \text{Bernoulli}\!\left(\text{logistic}(0.5\,W)\right), \quad
Y = 2A + 0.5W + \varepsilon, \quad
\varepsilon \sim N(0,1)
$$

In `causalsim_dgp()` terms: one standard-normal confounder, a constant effect of
2, moderate propensity confounding (logistic coefficient 0.5), and a moderate
baseline shift.

```{r dgp}
dgp <- causalsim_dgp(
  n = 500,
  n_confounders = 1,
  effect = 2,
  propensity = "moderate",
  baseline = "moderate"
)
dgp
```

<div class="callout">
The true ATE is exact when `effect` is a scalar; no Monte Carlo approximation
is needed. For function-valued effects, `causalsim_dgp()` approximates the ATE
via 10,000 Monte Carlo draws at construction time.
</div>

---

## Step 2: Inspect a Draw

`causalsim_draw()` simulates one dataset from the DGP. The columns `.tau` and
`.p` are the individual causal effect and propensity score, diagnostic metadata
that is not available in real observational data.

```{r draw}
dat <- causalsim_draw(dgp, seed = 1L)
head(dat)
```

With moderate confounding, the treated and control groups differ on the
pre-treatment covariate:

```{r confounding-check}
aggregate(W ~ A, data = dat, FUN = mean)
```

<div class="callout-warn">
That difference is exactly what makes naive regression biased. Any estimator
that omits `W` will absorb part of its association with `Y` into the treatment
coefficient.
</div>

---

## Step 3: Define Estimators

An estimator is any function that accepts the data frame returned by
`causalsim_draw()` and returns a named numeric vector. The `ci_lower` and
`ci_upper` fields are optional but enable coverage and power metrics.

```{r 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)
}
```

Named lists are also accepted, so the following is equivalent:

```{r 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]
  )
}
```

---

## Step 4: Evaluate Each Estimator

`causalsim_eval()` runs `reps` independent replications and returns a tidy
summary of bias, RMSE, coverage, and power with Monte Carlo standard errors.

```{r eval-naive}
eval_naive <- causalsim_eval(dgp, naive_est, reps = 300L, seed = 1L)
eval_naive
```

```{r eval-ols}
eval_ols <- causalsim_eval(dgp, ols_est, reps = 300L, seed = 1L)
eval_ols
```

The naive estimator's bias is substantial: treatment is positively correlated
with $W$, which also raises $Y$ through the baseline, so the unadjusted
coefficient absorbs part of that association. OLS eliminates the bias by
conditioning on $W$. Coverage for the naive estimator falls well below the
nominal 95% because the confidence intervals are centered on the wrong value.

`summary()` adds the full distribution of per-replication estimates, and
`plot()` shows it as a histogram:

```{r summary-ols}
summary(eval_ols)
```

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

---

## Step 5: Vary Sample Size with `causalsim_eval_grid()`

`causalsim_eval_grid()` evaluates an estimator over the Cartesian product of the
supplied parameter values, returning a tidy data frame of metrics for each cell.
Here we vary `n` across four levels to track how the OLS estimator's precision
improves with more data.

```{r 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
```

RMSE roughly halves as $n$ quadruples, consistent with $\sqrt{n}$-rate
convergence for OLS in a correctly specified model. Bias stays near zero at
every sample size.

---

## Step 6: Vary Confounding Strength

Varying the `propensity` preset shows how bias scales with confounding. Because
`causalsim_eval_grid()` accepts one estimator at a time, we run it separately and
combine the results.

```{r 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
```

Naive bias grows proportionally with confounding strength. OLS remains near zero
across all three levels because $W$ is observed and included in the model.

---

## Step 7: Heterogeneous Effects (Effect Modifiers)

So far the treatment effect has been constant. Real effects often vary across
subgroups. Declare a covariate with `role = "effect_modifier"` and reference it
in a function passed to `effect`:

```{r 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
```

The function passed to `effect` is what activates the modifier. A covariate
labelled `effect_modifier` but never referenced by `effect` is inert, and
`causalsim_dgp()` warns when that happens — so the role can never silently do
nothing.

Ground truth is carried per unit in `.tau`, so the subgroup effects are known
exactly:

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

An estimator that ignores the modifier recovers only the overall average
effect, while one that interacts treatment with `V` recovers the subgroup
effects:

```{r 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"]]
)
```

The constant-effect model lands near the true ATE (`r round(het_dgp$true_ate, 2)`),
but hides the heterogeneity; the interaction model recovers both subgroup
effects. This is the setup for benchmarking CATE / heterogeneous-effect
estimators (causal forests, meta-learners): because the true subgroup effects
are known, any estimator's recovery of them can be scored.

---

## Where to go next

This workflow (define, evaluate, grid) scales to more complex settings. A few
directions:

| Goal | How |
|---|---|
| Just generate data | Use `causalsim()` for a single dataset in one call |
| Heterogeneous effects | See Step 7 — declare an `effect_modifier` and pass a function to `effect` |
| Non-normal covariates | Use `causalsim_covar("binary")` or `causalsim_covar("uniform")` in `covariates` |
| Multiple confounders | Set `n_confounders = 3` or pass named `covariates` |
| Custom covariate structure | Mix `n_confounders` with explicit `covariates = list(...)` |

See `?causalsim`, `?causalsim_dgp`, and `?causalsim_covar` for the full API.

---

## Session info

```{r session-info}
sessionInfo()
```