causens: an R package for causal sensitivity analysis methods
In causens: Perform Causal Sensitivity Analyses Using Various Statistical Methods

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Introduction

Methods to estimate causal quantities often rely no unmeasured confounding. However, in practice, this is often not the case. causens is an R package that provides methods to estimate causal quantities in the presence of unmeasured confounding. In particular, causens implements methods from the three following school of thoughts:

Frequentist methods based on sensitivity functions
Bayesianism via parametric modelling
Monte Carlo sensitivity analysis

All implemented methods are accessible via the causens function with the method parameter. The package also provides a simulate_data function to generate synthetic data on which to test the methods.

Installation

The package can be installed from GitHub via devtools.

library(devtools)
devtools::install_github("Kuan-Liu-Lab/causens")

library(causens)

Methods

Summary of the Unmeasured Confounder Problem

In causal inference, the potential outcome framework is often used to define causal quantities. For instance, if a treatment $Z$ is binary, we respectively define $Y(1)$ and $Y(0)$ as the outcomes under treatment and control as if we can intervene on them a priori. However, in observational settings, this is often not the case, but we can still estimate estimands, e.g. the average treatment effect $\tau = \mathbb{E}[Y(1) - Y(0)]$, using observational data using causal assumptions that relate potential outcome variables to their observational counterpart:

Consistency: $Y = Y(1)Z + Y(0)(1 - Z)$
Conditional exchangeability: $Y(1), Y(0) \perp Z | X$
Positivity: $0 < P(Z = 1 | X) < 1$

where $X$ is a set of observed confounders that can be used to adjust for confounding. In practice, it is often difficult to know whether these assumptions hold, and in particular, whether all confounding variables are contained in $X$.

Hereon, we define $U$, the set of unmeasured confounders, although, in causens, we assume $U$ to be univariate for simplicity.

Simulated Data Mechanism

We posit the following data generating process:

\begin{align} X_1 &\sim \text{Normal}(0, 1) \ X_2 &\sim \text{Normal}(0, 1) \ X_3 &\sim \text{Normal}(0, 1) \ U &\sim \text{Bern}(0.5) \ Z &\sim \text{Bern}(\text{logit}^{-1}(\alpha_{uz} U + \alpha_{xz} X_3)) \ Y1 &\sim \text{Normal}(\beta_{uy} U + 0.5 X_1 + 0.5 X_2 - 0.5 X_3) \ Y0 &\sim \text{Normal}(0.5 X_1 + 0.5 X_2 - 0.5 X_3) \, . \end{align}

Using consistency, we define $Y = Y_1 Z + Y_0 (1 - Z)$. Note that, in simData.R, we provide some options to allow the simulation procedure to be more flexible. Parameters $\alpha_{uz}$ and $\beta_{uy}$ dictate how the unmeasured confounder $U$ affects the treatment $Z$ and the outcome $Y$, respectively; by default, they are set to 0.2 and 0.5.

Frequentist Methods (Brumback et al. 2004, Li et al. 2011)

Building on the potential outcome framework, the primary assumption that requires adjustment is conditional exchangeability. First, the latent conditional exchangeability can be articulated as follows to include $U$:

\begin{align} Y(1), Y(0) \perp Z | X, U \end{align}

Secondly, we can define the sensitivity function $c(z, e)$, with $z$ being a valid treatment indicator and $e$ being a propensity score value.

\begin{align} c(z, e) = \mathbb{E}\Big[Y(z) \mid Z = 1, e\Big] - \mathbb{E}\Big[Y(z) \mid Z = 0, e\Big] \end{align}

When there are no unmeasured confounders, $c(z, e) = 0$ for all $z$ and $e$ since $U = \emptyset$.

Given a valid sensitivity function, we can estimate the average treatment effect via a weighting approach akin to inverse probability weighting.

For parsimony, we provide an API for constant and linear sensitivity functions, but we eventually will allow any valid user-defined sensitivity function.

# Simulate data
data <- simulate_data(
  N = 10000, seed = 123, alpha_uz = 1,
  beta_uy = 1, treatment_effects = 1
)

# Treatment model is incorrect since U is "missing"
causens_sf(Z ~ X.1 + X.2 + X.3, "Y", data = data,
           c1 = 0.25, c0 = 0.25)$estimated_ate

plot_causens(Z ~ X.1 + X.2 + X.3, data, "Y",
             c1_upper = 0.5, c1_lower = 0, r = 1, by = 0.01)

Bayesian Methods

In Bayesianism, we can adjust for unmeasured confounding by explicitly modelling $U$ and its relationship with $Z$ and $Y$. Using a JAGS backend to carry out the MCMC procedure, we can estimate the average treatment effect by then marginalizing over $U$. We assume the following data-generating mechanism and Bayesian model in {causens}:

\begin{align} X_1, X_2, X_3 &\stackrel{iid}{\sim} \text{Normal}(0, 1) \ U \mid X &\sim \operatorname{Bern}\left(\text{logit}^{-1}\gamma_{ux} X\right) \ Z \mid U, X &\sim \text{Bern}\left(\text{logit}^{-1}(\alpha_{uz} U + \alpha_{xz} X)\right) \ Y \mid U, X, Z &\sim \text{Normal}(\beta_{uy} U + \beta_{xy} X + \tau Z) \ \end{align}

``` {r, eval = FALSE} data <- simulate_data( N = 1000, alpha_uz = 0.5, beta_uy = 0.2, seed = 123, treatment_effects = 1, y_type = "continuous" )

bayesian_causens( Z ~ X.1 + X.2 + X.3, Y ~ X.1 + X.2 + X.3, U ~ X.1 + X.2 + X.3, data )

## Monte Carlo Approach to Causal Sensitivity Analysis

**Warning** Only works for binary outcomes, for now.

The Monte Carlo approach proceeds through the following steps:

- Produce a "naive" model where $U$ is ignored: $$\text{logit}\Big(P(Y = 1 \mid Z, X)\Big) = \beta_0 + \beta_Z Z + \beta_X X$$.
- With estimates $\widehat{\beta}_Z$ and $\widehat{\text{SD}}(\widehat{\beta}_Z)$ of $\beta_Z$ and $\text{SD}(\widehat{\beta}_Z)$ obtained from fitting the naive model,
we can sample $\beta_U, \gamma_0, \gamma_Z$ from an assumed prior distribution:
\begin{align*}
  \beta_U &\sim \text{Uniform}(-2, 2) \\
  \gamma_0 &\sim \text{Uniform}(-5, 5) \\
  \gamma_Z &\sim \text{Uniform}(-2, 2) \, .
\end{align*}
- For each sample $\left(\beta_U^{(r)}, \gamma_0^{(r)}, \gamma_Z^{(r)}\right)$ indexed by $r = 1, \dots, R$, we can compute bias adjusted $\widehat{\beta}_Z^{(r)}$ estimates:
$$\widehat{\beta}_Z^{(r)} = \widehat{\beta}_Z - \log \left(\frac{(1 + \exp(\beta_U^{(r)} + \gamma_0^{(r)} + \gamma_X^{(r)}))(1 + \exp(\gamma_0^{(r)}))}{(1 + \exp(\beta_U^{(r)} + \gamma_0^{(r)}))(1 + \exp(\gamma_0^{(r)} + \gamma_X^{(r)}))}\right)$$
- Finally, we can sample the bias adjusted average treatment effect estimate
$$\beta^{(r)} \sim \text{Normal}\left(\widehat{\beta}_Z^{(r)}, \widehat{SD}\left(\widehat{\beta}_Z\right)^2\right)$$ to account for errors incurred in previous steps. With $\beta^{(1)}_Z, \dots, \beta^{(R)}_Z$, we can compute the Monte Carlo estimate of the average treatment effect by taking their average.

```r
data <- simulate_data(
  N = 1000, alpha_uz = 0.2, beta_uy = 0.5,
  seed = 123, treatment_effects = 1, y_type = "binary",
  informative_u = FALSE
)

causens_monte_carlo("Y", "Z", c("X.1", "X.2", "X.3"), data)$estimated_ate