# Correcting for selection bias with standardisation/g-computation

Two sets of methods for correcting for selection bias are g-computation (standardisation) and inverse probability of censoring weighting (IPCW). I'm having a difficult time understanding how to apply the former.

As a demonstration, I generate some data from the following DAG - Figure 8.4 from Hernan and Robins Causal Inference:

$$A$$ is a binary treatment, $$Y$$ is a continuous outcome, $$C$$ indicates censoring (e.g. loss to follow-up), and $$U$$ is an unmeasured variable. The true causal effect is zero (there are no paths from A into Y).

library(tidyverse)

sim_data <- function(n,seed){
set.seed(seed)
A <- rbinom(n, 1, 0.5)
U <- rnorm(n)
L <- rbinom(n, 1, plogis(-0.5 + A + U))
C <- rbinom(n, 1, plogis(-1 + 3*L))
Y <- U + rnorm(n)
tibble(A,L,C,Y)
}

df <- sim_data(n=100000,seed=123)


Recovering the true treatment effect requires no backdoor paths between $$Y$$ and $$A$$ (which holds unconditional on $$L$$, a collider) and no backdoor paths between $$Y$$ and $$C$$ (which requires adjustment for $$L$$).

Recovering the true treatment effect seems pretty straightforward by IPCW:

# IPC weighting
C_model <-
glm(C==0 ~ L, data = df, family = binomial)

df$ipc_weights <- 1/predict(C_model, newdata = df, type = "response") Y_model <- lm(Y ~ A, data = filter(df, C == 0), weights = filter(df, C == 0)$ipc_weights)

coef(Y_model)["A"]
#0.0001701217


Another option is g-computation. The formula given in this paper, as I've (mis)interpreted it and adapted to this example, is: $$\sum_lE[Y|A=1,L=l,C=0]-E[Y|A=0,L=l,C=0]Pr[L=l]$$

When I apply this approach, however, I'm consistently getting the wrong answer. E.g.

# G-computation
df_c <- filter(df, C == 0)

treated_L0 <-
mean(df_c[df_c$$A == 1 & df_c$$L==0,]$$Y) treated_L1 <- mean(df_c[df_c$$A == 1 & df_c$$L==1,]$$Y)
untreated_L0 <-
mean(df_c[df_c$$A == 0 & df_c$$L==0,]$$Y) untreated_L1 <- mean(df_c[df_c$$A == 0 & df_c$$L==1,]$$Y)

Pr_L0 <- sum(df$L==0)/nrow(df) Pr_L1 <- 1 - Pr_L0 out <- ((treated_L0 - untreated_L0)*Pr_L0) + ((treated_L1 - untreated_L1)*Pr_L1) out # -0.1629754  Where am I going wrong here? My instinct is that I should not be summing over $$Pr[L=l]$$ but instead over $$Pr[L=l|A=a]$$ (which does seem to give the right answer), but I don't really understand why, if that is the case. Any help much appreciated! • I don't think the causal effect is identified with g-computation here. Conditioning on$L\$ doesn't remove the selection bias.
– Noah
Oct 9, 2023 at 23:58
• @Noah You may well be right. I had the intuitive impression that because IP weighting and standardisation are mathematically equivalent, there couldn't be a situation in one worked but not the other. But there seems to be no way around conditioning on L with g-computation. Oct 10, 2023 at 0:11
• Yep, Hernan and Robins actually address the exact case you're mentioning later on (see Ch. 8.5) -- in short, for IPCW, you're not really "directly" conditioning on L when you come up with your final estimate of the ATE. Oct 10, 2023 at 1:24

I appreciate both comments to the question, but standardization/g-computation is identified here. The issue is that the formula in this causal structure differs from the usual formula.

In the paper by Breskin et al. 2018 they consider a DAG similar to the one provided above (except they have an arrow from $$A \rightarrow Y$$).

They give the following equality (translated to the variable names given here): $$E[Y^a] = \sum_l E[Y | A=a,L=l,C=0] \Pr(L=l | A=a)$$ This is distinct from the usual g-computation formula we usually see. See the proof in the paper for the math and ideas that lead to this form. They also demonstrate this with simulations.

A similar proof is given in Example 1 of my paper here. It considers a modification of the DAG shown above, but shares some similarities.

Update:

I had some extra time so I am returning to this question to show how the corresponding g-computation algorithm can be implemented. I am going to do this with Python and my package delicatessen for fitting the models. The data generating mechanism comes from the Breskin et al. paper referenced above.

import numpy as np
import pandas as pd
from delicatessen import MEstimator
from delicatessen.estimating_equations import ee_regression, ee_glm
from delicatessen.utilities import inverse_logit

# Generating the data
np.random.seed(20231107)
n = 1000000
d = pd.DataFrame()
d['A'] = np.random.binomial(n=1, p=0.5, size=n)
d['U'] = np.random.binomial(n=1, p=0.3, size=n)
d['W'] = np.random.binomial(n=1, p=0.2 + 0.4*d['A'] + 0.3*d['U'], size=n)
d['S'] = np.random.binomial(n=1, p=0.1 + 0.8*d['W'], size=n)
d['Y'] = np.random.binomial(n=1, p=0.3 + 0.5*d['U'] - 0.2*d['A'], size=n)
d['I'] = 1
ds = d.loc[d['S'] == 0].copy()


As a comparison, we can use the full data on Y to assess the causal effect

# Truth
X = np.asarray(d[['I', 'A']])
y = np.asarray(d['Y'])

def psi(theta):
return ee_glm(theta, X=X, y=y, distribution='binomial', link='identity')

estr = MEstimator(psi, init=[0.5, -0.2])
estr.estimate()


which gives the estimated causal effect as -0.20066583. When ignoring the selection bias, we get an answer of -0.23594802 and when adjusting for W via a model we get -0.25898981. The following code provides an implementation of g-computation for this type of selection bias problem as an estimating equation

# G-computation
d['Y'] = np.where(d['S'] == 0, d['Y'], -9999)
d['AW'] = d['A'] * d['W']
X = np.asarray(d[['I', 'A', 'W', 'AW']])
y = np.asarray(d['Y'])
a = np.asarray(d['A'])
s = np.asarray(d['S'])
da = d.copy()
da['A'] = 1
da['AW'] = d['W']
X1 = np.asarray(da[['I', 'A', 'W', 'AW']])
da['A'] = 0
da['AW'] = 0
X0 = np.asarray(da[['I', 'A', 'W', 'AW']])

def psi(theta):
# Subsetting parameters out
rd, r1, r0 = theta[0:3]
beta = theta[3:]

# Outcome nuisance model, Pr(Y | A,W,S=0)
ee_out = ee_regression(beta, X=X, y=y, model='logistic')
ee_out = ee_out * (1-s)  # subsetting to S=0
y1hat = inverse_logit(np.dot(X1, beta))  # Predicted Y^{1}
y0hat = inverse_logit(np.dot(X0, beta))  # Predicted Y^{0}

# Risk functions
ee_r1 = a*(y1hat - r1)         # Mean Y^{1} among A=1
ee_r0 = (1-a)*(y0hat - r0)     # Mean Y^{0} among A=0
ee_rd = np.ones(y.shape) * (r1 - r0 - rd)  # ACE

return np.vstack([ee_rd, ee_r1, ee_r0, ee_out])

estr = MEstimator(psi, init=[0., 0.5, 0.5, -0.3, -1., 0.6, 0.3])
estr.estimate()


(starting values are mildly informative to speed up the fitting process). Here, the procedure estimates the average causal effect as -0.19898929 which is close to the truth.

• My intuition was that we should be averaging over Pr[L=l|A=a] and for once it's turned out to be right... I'm reading the cited papers now for the details now. Thanks! Oct 11, 2023 at 4:43