Table 2 Fallacy - Example

table2_fallacy_ex.R

# Title: Table 2 Fallacy - Example ----

# Description: When interpreting results, if the goal is causal inference, 
# it's important to consider pathways between the exposure and outcome. 
# If the exposure of interest changes, so do the pathways. This is known as 
# the Table 2 Fallacy. 

# This code demonstrates the Table 2 Fallacy. Highly recommend the paper 
# The table 2 fallacy: presenting and interpreting confounder and modifier coefficients
# by Westreich & Greenland 

# Setup ----

#... Libraries ----

library(tidyverse) # ol' faithful
library(ggdag) # for drawing DAG 
library(WeightIt) # for IP weighting
library(broom) # tidying output from model

# Simulating Data ----

# Two confounders: z1 & z2
# Additional variable that is a confounder for the effect of sleep on happiness (z3)
# Binary Exposure 
# Continuous outcome 

# Note: variable name to DAG
# z1 - sleep
# z2 - sunshine
# z3 - lifting
# x - coffee
# y - happiness

set.seed(456) # for reproducibility

n = 250 # arbitrary sample size 

df <- data.frame(
  z3 = rnorm(n = n, mean = 5, sd = 2), 
  z2 = rbinom(n = n, size = 1, prob = 0.45)
) %>% 
  dplyr::mutate(
    z1 = rnorm(n = n, mean = 0.5*z3 + 1, sd = 2),
    prob = plogis(0.05*z1 + 0.2*z2), 
    x = rbinom(n = n, size = 1, prob = prob), 
    y = 1.5*x + 2*z1 + 0.5*z2 + 1.3*z3 + rnorm(n = n, mean = 0, sd = 1),
  )

# Adjusting for Confounders ----

# Using a doubly robust estimator, where the covariates in the propensity score 
# model are included in the outcome model. The target estimand will be the 
# average treatment effect. 

#... IP Weighting ----

ps.mod <- WeightIt::weightit(x ~ z1 + z2, 
                             data = df, 
                             method = "glm", 
                             estimand = "ATE",
                             stabilize = TRUE)

#... Outcome Model ----

glm_weightit(y ~ x + z1 + z2,
             data = df, 
             family = gaussian(), 
             weights = ps.mod$weights) %>% 
  broom::tidy(conf.int = TRUE, conf.level = 0.95) %>% 
  select(term, estimate, conf.low, conf.high)

# Effect of Z1 on Happiness ----- 

# If we were to estimate the effect of z1 on happiness, 
# we would need to adjust for z3. Demonstrated below 

# For simplicity, we will use only an outcome model. 

glm(y ~ z1 + z3, 
    data = df, 
    family = gaussian()) %>% 
  broom::tidy(conf.int = TRUE)
  
# From the above, you can see this estimate is closer to the "true" effect for Z1