# Correct methodology using $g$-computation to estimate Average Treatment Effect on the Treated ($ATT$)?

I have a question about the $$g$$-computation methodology for estimating the Average Treatment Effect on the Treated ($$ATT$$) in the following article.

The authors recommend estimating the $$ATT$$ by first performing a regression model for $$Y$$ on the binary treatment $$A$$ and covariates $$X$$, then finding $$E(Y^1|A=1)$$ from the observed values of $$Y|A=1$$ and estimating the counterfactual $$E(Y^0|A=1)$$ by imputing the predicted value of $$Y$$ from the regression formula by substituting $$A=0$$ into the regression formula, but only for the observed $$A=1$$ records. So far so good.

For the next step, Wang et al. stack the records with the $$E(Y^1|A=1)$$ and $$E(Y^0|A=1)$$ estimates into a single dataset coupled with their $$A$$ values of 1 and 0, respectively. They then regress $$Y$$ on $$A$$, where $$ATT$$ is estimated from the regression coefficient on $$A$$.

Rather than regress $$Y^*$$ on $$A$$ why not simply compute the average of the difference $$E(Y^1|A=1) - E(Y^0|A=1)$$ for all records where $$A=1$$?

When I run a simulation the two methods return different results for the ATT point estimates - see R code below on a simulated dataset.

# Simulation data from "Targeted Maximum Likelihood Estimation for Causal Inference in    #
# Observational Studies", Schuler & Rose, 2016.                                           #
x1 <- rbinom(n=10000, size=1, prob=0.55)
x2 <- rbinom(n=10000, size=1, prob=0.3)

# Binary treatment variable
A <- rbinom(n=10000, size=1, prob=exp(-.5 + .75*x1 + x2)/(1 + exp(-.5 + .75*x1 + x2)))

# Outcome variable
Y <- rnorm(n=10000, mean=24 - 3*A + 3*x1 - 4*x2 + 7*x1*x2 + 5*A*x1 - 10*A*x2 + 15*A*x1*x2, sd=4.5)

# Create data frame
df <- data.frame(cbind(Y, A, x1, x2))

### Run linear regression
lm_test <- lm(Y ~ . + x1*x2 + A*x1 + A*x2 + A*x1*x2, data=df)

### Calculate ATT using g-computation. ###
# ATT = Expected difference between observed and counterfactual Y values at record level. #
df_trmt <- df[df$$A==1,] df_a0 <- df_trmt; df_a0$$A = 0;
y_a0 = predict(lm_test, newdata=df_a0)
y_a1 = df_trmt$Y df_preds <- data.frame(cbind(df_trmt, y_a0, y_a1)) df_preds$$att <- y_a1 - y_a0 mean(df_preds$$att) # [1] -0.284822 <==== ATT estimate sqrt(var(df_preds$$att)/length(df_preds$$att)) # [1] 0.08707805 # Alternate g-computation calculation of ATT (based on Wang et al. (2017) article). # df_a0_2 = data.frame(cbind(Y=y_a0, A=df_a0$$A)) df_a1_2 = data.frame(cbind(Y=y_a1, A=df_trmt$$A)) df_preds_2 <- data.frame(rbind(df_a0_2, df_preds_2)) lm_test_2 <- lm(Y ~ A, data=df2) summary(lm_test_2) # Call: # lm(formula = Y ~ A, data = df2) # # Residuals: # Min 1Q Median 3Q Max # -32.536 -4.337 0.592 5.399 24.712 # # Coefficients: # Estimate Std. Error t value Pr(>|t|) # (Intercept) 25.1946 0.1265 199.207 < 2e-16 *** # A 0.5096 0.1723 2.957 0.00311 ** <==== ATT estimate # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # # Residual standard error: 8.591 on 9998 degrees of freedom # Multiple R-squared: 0.0008738, Adjusted R-squared: 0.0007739 # F-statistic: 8.744 on 1 and 9998 DF, p-value: 0.003113  ## 1 Answer I think you have a bug in your code. I don't remember what I changed, but your code doesn't even run in a fresh R session. Logic demands that the results should be equal, and that is what I get: set.seed(666) # Simulation data from "Targeted Maximum Likelihood Estimation for Causal Inference in # # Observational Studies", Schuler & Rose, 2016. # x1 <- rbinom(n=10000, size=1, prob=0.55) x2 <- rbinom(n=10000, size=1, prob=0.3) # Binary treatment variable A <- rbinom(n=10000, size=1, prob=exp(-.5 + .75*x1 + x2)/(1 + exp(-.5 + .75*x1 + x2))) # Outcome variable Y <- rnorm(n=10000, mean=24 - 3*A + 3*x1 - 4*x2 + 7*x1*x2 + 5*A*x1 - 10*A*x2 + 15*A*x1*x2, sd=4.5) # Create data frame df <- data.frame(cbind(Y, A, x1, x2)) ### Run linear regression lm_test <- lm(Y ~ . + x1*x2 + A*x1 + A*x2 + A*x1*x2, data=df) ### Calculate ATT using g-computation. ### # ATT = Expected difference between observed and counterfactual Y values at record level. # df_trmt <- df[df$$A==1,] df_a0 <- df_trmt; df_a0$$A = 0; y_a0 = predict(lm_test, newdata=df_a0) y_a1 = df_trmt$Y

df_preds <- data.frame(cbind(df_trmt, y_a0, y_a1))

df_preds\$att <- y_a1 - y_a0

# Alternate g-computation calculation of ATT (based on Wang et al. (2017) article). #
df_a0_2 = data.frame(cbind(Y=y_a0, A=0))
df_a1_2 = data.frame(cbind(Y=y_a1, A=1))
df_preds_2 <- data.frame(rbind(df_a0_2, df_a1_2))
lm_test_2 <- lm(Y ~ A, data=df_preds_2)

summary(lm_test_2)
# Call:
#   lm(formula = Y ~ A, data = df_preds_2)
#
# Residuals:
#   Min      1Q  Median      3Q     Max
# -33.951  -2.485   1.001   4.143  24.921
#
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)
# (Intercept)  26.0230     0.1060 245.597   <2e-16 ***
#   A            -0.3081     0.1498  -2.056   0.0398 *
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 7.911 on 11146 degrees of freedom
# Multiple R-squared:  0.0003792,   Adjusted R-squared:  0.0002896
# F-statistic: 4.229 on 1 and 11146 DF,  p-value: 0.03977
mean(y_a1 - y_a0)
# [1] -0.3081427


The standard errors don't agree, but of course neither standard error is correct, so this isn't much of a concern to me.

• That worked, thank you! Not sure where I messed up in my code. Yea the standard errors aren't correct - in the paper the authors used bootstrapping to estimate the SE. Commented May 18, 2021 at 1:04