1
$\begingroup$

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
$\endgroup$

1 Answer 1

1
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – RobertF
    Commented May 18, 2021 at 1:04

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.