Simulating potential outcomes with a binary outcome

Question

I want to create some simple simulations of potential outcomes to explore issues of confounding. I start with a binary confounder X and a binary treatment A. When my outcome is continuous, I can create the data generating model for the outcome Yc, and then simulate the counterfactuals by replacing A with 0 or 1. In this case, Y=A*Y1.c + (1-A)*Y0.c is always true (the potential outcome under a given treatment A=a is always equal to the observed outcome when A=a).

library(tidyverse)

set.seed(0324)

x<- rbinom(100, 1, 0.3)
A<- rbinom(100, 1, (0.1 + 0.2*x))
Yc<- 4 + 3*x + 4*A
Y1.c<- 4 + 3*x + 4*1
Y0.c <-  4 + 3*x + 4*0

However, when I have a binary outcome, I'm not sure how to execute the same process. I could start by simulating Y as a bernoulli random variable with probability equal to some some linear combination of my treatment and confounder, where when x=0 and A=0 the probability is 0.1, when x=1 and A=0 the probability is 0.15, etc.

Yb <- rbinom(100, 1, (0.1 +0.05*x + 0.2*A))

However, this model is not deterministic like the linear equation was. As a result, Y=A*Y1.c + (1-A)*Y0.c won't always be true (and in the below example, it's not true 27% of the time).

Y1.b <- rbinom(100, 1, (0.1 + 0.05*x + 0.2*1))
Y0.b <-rbinom(100, 1, (0.1 + 0.05*x + 0.2*0))

df <- data.frame(x,A,Yc, Y1.c, Y0.c, Yb, Y1.b, Y0.b) %>%
  mutate(discordant.b=ifelse(A==1 & Y1.b != Yb, 1,
                         ifelse(A==0 & Y0.b != Yb, 1, 0))) %>%
  mutate(discordant.c=ifelse(A==1 & Y1.c != Yc, 1,
                         ifelse(A==0 & Y0.c != Yc, 1, 0)))

summarize(df, mean(discordant.b), mean(discordant.c))

#  mean(discordant.b) mean(discordant.c)
1               0.27                  0

I can simulate it the other way, first by creating the potential outcomes and then by simply defining Y based the potential outcome

df <- df %>% mutate (Ynew= ifelse(A=1, Yb.1, Yb.0))

But this feels "backwards" to me. Am I missing something? Is there a similar deterministic approach to a binary outcome? Or is the continuous outcome a special case where it's possible?

Answer 1

You are indeed simulating the data backward. You should simulate the potential outcomes first and then the outcomes from them. The "consistency equation" $Y = AY_1 + (1-A)Y_0$ describes the data-generating process for $Y$, the observed outcomes. The potential outcomes "occur" before the treatment is assigned in the formulation of causal inference that involves them. Potential outcomes exist for each individual prior to treatment assignment, and then treatment "reveals" one of the potential outcomes via the consistency equation. See my answer here about this very topic. This is also how you should generate your data.

The reason it appears to work in the reverse-ordered way you generated the potential outcomes for the continuous case is that you omitted a stochastic component to the potential outcomes, which is unrealistic. That is, if you generated your observed and potential outcomes more realistically as

Yc   <- 4 + 3*x + 4*A + rnorm(100)
Y1.c <- 4 + 3*x + 4*1 + rnorm(100)
Y0.c <- 4 + 3*x + 4*0 + rnorm(100)

you would also find disagreement between the generated observed outcome and the outcome implied by the consistency equation applied to the generated potential outcomes. Instead, you should generate your outcomes as follows:

Y1.c <- 4 + 3*x + 4*1 + rnorm(100)
Y0.c <- 4 + 3*x + 4*0 + rnorm(100)
Yc   <- A*Y1.c + (1-A)*Y0.c

which is more in line with the data-generating process implied by the consistency equation. For binary outcomes, you should do the same. Note that if you considered $p_Y = P(Y_i = 1)$ as the outcome in the binary case (i.e., omitting the step where you turn the probabilities into zeroes and ones using rbinom()), the probabilities would satisfy the consistency equation even in your reversed order because there is no stochastic component in generating them.