---
title: "Causal Inference"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Causal Inference}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r knitr-opts, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library(deli)
```

## Overview

deli provides estimating equations for several causal inference methods.
A key advantage of the M-estimation framework is that all nuisance model
parameters (e.g., propensity scores, outcome models) are estimated
jointly, so the sandwich variance correctly propagates uncertainty from
all stages.

This vignette demonstrates:

- G-formula (g-computation / standardization)
- Inverse probability weighting (IPW)
- Augmented IPW (AIPW / doubly robust)

## Simulated data

We simulate data with a binary treatment `A`, a confounder `W`, and a
continuous outcome `Y`:

```{r simulate-data}
set.seed(42)
n <- 1000
W1 <- rnorm(n)
W2 <- rbinom(n, 1, 0.4)
A <- rbinom(n, 1, plogis(-0.5 + 0.5 * W1 + 0.3 * W2))
# True ATE = 1.5
Y <- 2 + 1.5 * A + W1 - 0.5 * W2 + rnorm(n)
```

## G-formula (g-computation)

The g-formula estimates the average causal effect by fitting an outcome
model and standardizing (averaging predictions over the covariate
distribution):

```{r gformula}
X <- cbind(1, A, W1, W2)      # Observed design matrix
X1 <- cbind(1, 1, W1, W2)     # Counterfactual: all treated
X0 <- cbind(1, 0, W1, W2)     # Counterfactual: all untreated

psi_gformula <- function(theta) {
  ee_gformula(theta, y = Y, X = X, X1 = X1, X0 = X0)
}

# theta: ACE, E[Y(1)], E[Y(0)], beta0, beta1, beta2, beta3
m_gformula <- m_estimate(
  stacked_equations = psi_gformula,
  init = c(0, 0, 0, 0, 0, 0, 0)
)

# ACE estimate (true = 1.5)
m_gformula@theta[1]

# With confidence interval
summary(m_gformula)
```

The first parameter is the average causal effect (ACE).

## Inverse probability weighting (IPW)

IPW reweights the observed data by the inverse of the treatment
probability (propensity score):

```{r ipw}
W_ps <- cbind(1, W1, W2)  # Propensity score design matrix

psi_ipw <- function(theta) {
  ee_ipw(theta, y = Y, A = A, W = W_ps)
}

# theta: ACE, E[Y(1)], E[Y(0)], alpha0, alpha1, alpha2
m_ipw <- m_estimate(stacked_equations = psi_ipw, init = c(0, 0, 0, 0, 0, 0))

# ACE
m_ipw@theta[1]
summary(m_ipw)
```

### Truncating propensity scores

Propensity scores near zero or one produce very large weights, and a
handful of large weights can dominate an IPW estimate. `truncate` clips
the fitted scores to a range before they are inverted. Whether that
changes anything depends on the scores themselves, so look at them
first:

```{r ps-range}
# Propensity scores implied by the alphas in the fitted stack
ps <- plogis(drop(W_ps %*% m_ipw@theta[4:6]))
range(ps)

# How many scores a tight c(0.3, 0.7) range would clip, by tail
c(lower = sum(ps < 0.3), upper = sum(ps > 0.7))
```

Nothing here is close to zero or one, so a conventional range such as
`c(0.1, 0.9)` would leave every score untouched and the fit unchanged. A
deliberately tight range shows what truncation does. `c(0.3, 0.7)` clips
166 of the 1000 scores, all but nine of them in the lower tail:

```{r ipw-truncate}
psi_ipw_trunc <- function(theta) {
  ee_ipw(theta, y = Y, A = A, W = W_ps,
         truncate = c(0.3, 0.7))
}

m_ipw_trunc <- m_estimate(
  stacked_equations = psi_ipw_trunc,
  init = c(0, 0, 0, 0, 0, 0)
)
summary(m_ipw_trunc)
```

Both the ACE and its standard error move relative to the untruncated
fit:

```{r compare-truncation}
data.frame(
  fit = c("Untruncated", "Truncated"),
  estimate = round(c(m_ipw@theta[[1]], m_ipw_trunc@theta[[1]]), 4),
  std_err = round(
    sqrt(c(vcov(m_ipw)[1, 1], vcov(m_ipw_trunc)[1, 1])),
    4
  )
)
```

Truncation is usually motivated as a variance reduction, so the rise in
the standard error is worth dwelling on. Clipping pays for itself only
when the weights it caps are genuinely extreme, which is not the case
here. It also costs precision that joint estimation would otherwise
supply: a clipped score no longer moves with the propensity score
parameters, so the sandwich variance credits less of the gain from
estimating those parameters rather than treating them as known. Compare
a truncated fit against the untruncated one instead of reaching for a
range by default.

## Augmented IPW (doubly robust)

AIPW combines outcome modeling and IPW. It is "doubly robust": the
estimate is consistent if either the outcome model or the propensity
score model is correctly specified. It reuses both sets of design
matrices built above:

```{r aipw}
psi_aipw <- function(theta) {
  ee_aipw(theta, y = Y, A = A,
          W = W_ps,             # Propensity score model
          X = X, X1 = X1, X0 = X0)  # Outcome model
}

# theta: ACE, E[Y(1)], E[Y(0)], alpha (3), beta (4)
m_aipw <- m_estimate(
  stacked_equations = psi_aipw,
  init = c(0, 0, 0, rep(0, 3), rep(0, 4))
)

summary(m_aipw)
```

## Comparing methods

The g-formula, IPW, and AIPW fits above all target the same average
causal effect. Both the outcome model and the propensity score model are
correctly specified in this simulation, so the estimates agree closely
with one another and with the true value:

```{r compare-methods}
data.frame(
  method = c("G-formula", "IPW", "AIPW"),
  estimate = round(
    c(m_gformula@theta[[1]], m_ipw@theta[[1]], m_aipw@theta[[1]]),
    3
  ),
  true_ate = 1.5
)
```

## Further reading

- `?ee_ipw_msm`: Marginal structural models with IPW
- `?ee_gestimation_snmm`: G-estimation for structural nested mean models
- `?ee_iv_causal`: Instrumental variable estimation
- `?ee_2sls`: Two-stage least squares
- `?ee_mean_sensitivity_analysis`: Sensitivity analysis for unmeasured
  confounding