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: Figure 8.4 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).


sim_data <- function(n,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)

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)


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)

# -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!

  • 2
    $\begingroup$ I don't think the causal effect is identified with g-computation here. Conditioning on $L$ doesn't remove the selection bias. $\endgroup$
    – Noah
    Oct 9, 2023 at 23:58
  • $\begingroup$ @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. $\endgroup$
    – Lachlan
    Oct 10, 2023 at 0:11
  • 1
    $\begingroup$ 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. $\endgroup$ Oct 10, 2023 at 1:24

1 Answer 1


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.


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
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])

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])

(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.

  • 1
    $\begingroup$ 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! $\endgroup$
    – Lachlan
    Oct 11, 2023 at 4:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.