Distribution to match an example with collider bias?

Question

I am doing exercises from "Causal Inference in Statistics: A Primer", by Pearl et al (2016). In chapter 1.2 there is a training challenge that goes like:

In an attempt to estimate the effectiveness of a new drug, a randomized experiment is conducted. In all, 50% of the patients are assigned to receive the new drug and 50% to receive a placebo. A day before the actual experiment, a nurse hands out lollipops to some patients who show signs of depression, mostly among those who have been assigned to treatment the next day (i.e., the nurse’s round happened to take her through the treatment-bound ward). [..] Although the drug proved beneficial to the population as a whole, drug takers were less likely to recover than nontakers, among both lollipop receivers and lollipop nonreceivers. Assuming that lollipop sucking in itself has no effect whatsoever on recovery [..] Is the drug beneficial to the population as a whole or harmful?

Trying to understand the situation, I've drawn a DAG on www.dagitty.net like below. daggity identifies that lollies is a confounder. I've drawn the conclusion that I should not control for lollies. The drug is good on the whole.

But I want to have deeper understanding of this Simpsons Paradox! With numeric experimentation I can generate a distribution that gives a positive effect on the population as a whole, and a negative effect on the lollie-recievers. But I cannot come up with any distribution that matches the situation described in the question (positive effect of drug on outcome as a whole, but negative for both groups with lollie recievers and nonrecievers)

My question: Is there any such distribution, or is there any provable reason why I cannot come up ith it?

Answer 1

You probably just need to make the collider bias stronger in your simulations. To help you get started, here's a simple example that obtains the desired reversal:

set.seed(10)
n     <- 1e6
x     <- rbinom(n, 1, 0.5)
dep   <- rbinom(n, 1, 0.5)
lolli <- rbinom(n, 1, 0.5*x + 0.5*dep)
y     <- rbinom(n, 1, 0.1*x + 0.8*dep)

mean(y[x==1]) - mean(y[x==0]) # 0.1
mean(y[x==1 & lolli==1]) - mean(y[x==0 & lolli==1]) # -0.17
mean(y[x==1 & lolli==0]) - mean(y[x==0 & lolli==0]) # -0.17

As you can see, the correct causal effect of the drug (0.1) is obtained when looking at the aggregated data. If you, however, look at each subpopulation, you would wrongly conclude the drug is harmful (-0.17).

You can find more about Simpson's paradox in this answer too.