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?