Treatment Effect Estimators as Weighted Outcomes

Application 401(k) - average effects

Replication comments:

• The many required NxN smoother matrices make this notebook relatively memory intensive (32GB RAM should suffice).

• Running the notebook within the replication docker ensures that results perfectly replicate. Otherwise, results might differ depending on which package versions you use.

This notebook runs the application described in Section 5.1 of Knaus (2024). The first part replicates the results presented in the paper, the second part provides supplementary information

1. Replication of paper results

Getting started

First, load packages and set the seed:

if (!require("OutcomeWeights")) install.packages("OutcomeWeights", dependencies = TRUE); library(OutcomeWeights)
if (!require("hdm")) install.packages("hdm", dependencies = TRUE); library(hdm)
if (!require("grf")) install.packages("grf", dependencies = TRUE); library(grf)
if (!require("cobalt")) install.packages("cobalt", dependencies = TRUE); library(cobalt)
if (!require("tidyverse")) install.packages("tidyverse", dependencies = TRUE); library(tidyverse)
if (!require("viridis")) install.packages("viridis", dependencies = TRUE); library(viridis)
if (!require("gridExtra")) install.packages("gridExtra", dependencies = TRUE); library(gridExtra)

set.seed(1234)

Next, load the data. Here we use the 401(k) data of the hdm package. However, you can adapt the following code chunk to load any suitable data of your choice. Just make sure to call the treatment D, covariates X, and instrument Z. The rest of the notebook should run without further modifications.

data(pension) # Find variable description if you type ?pension in console

# Treatment
D = pension$p401
# Instrument
Z = pension$e401
# Outcome
Y = pension$net_tfa
# Controls
X = model.matrix(~ 0 + age + db + educ + fsize + hown + inc + male + marr + pira + twoearn, data = pension)
var_nm = c("Age","Benefit pension","Education","Family size","Home owner","Income","Male","Married","IRA","Two earners")
colnames(X) = var_nm

Run Double ML

In the following we run double ML with default honest random forest (tuning only increases running time without changing the insights in this application). As standard implementations do currently not allow to extract the outcome smoother matrices, the OutcomeWeights package comes with a tailored internal implementation called dml_with_smoother(), which is used in the following.

2-folds

First, we run all estimators with 2-fold cross-fitting:

# 2 folds
dml_2f = dml_with_smoother(Y,D,X,Z,
                           n_cf_folds = 2)
results_dml_2f = summary(dml_2f)

          Estimate      SE      t         p    
PLR        13903.3  1539.3 9.0323 < 2.2e-16 ***
PLR-IV     12931.8  1914.3 6.7554 1.504e-11 ***
AIPW-ATE   11680.5  1185.0 9.8566 < 2.2e-16 ***
Wald-AIPW  11332.9  1661.0 6.8228 9.448e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

plot(dml_2f)

Now, we use the get_outcome_weights() method to extract the outcome weights as described in the paper. To illustrate that the algebraic results provided in the paper are indeed numerical equivalences and no approximations, we check whether the weights multiplied by the outcome vector reproduces the conventionally generated point estimates.

omega_dml_2f = get_outcome_weights(dml_2f)
cat("ω'Y replicates point etimates?", 
    all.equal(as.numeric(omega_dml_2f$omega %*% Y),
      as.numeric(results_dml_2f[,1])
    ))

ω'Y replicates point etimates? TRUE

5-fold

Run double ML also with 5-fold cross-fitting:

# 5 folds
dml_5f = dml_with_smoother(Y,D,X,Z,
                           n_cf_folds = 5)
results_dml_5f = summary(dml_5f)

          Estimate      SE      t         p    
PLR        13898.6  1525.3 9.1123 < 2.2e-16 ***
PLR-IV     13194.8  1902.7 6.9348 4.322e-12 ***
AIPW-ATE   11574.8  1158.6 9.9903 < 2.2e-16 ***
Wald-AIPW  11386.7  1638.4 6.9499 3.887e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

plot(dml_5f)

extract the weights and confirm numerical equivalence:

omega_dml_5f = get_outcome_weights(dml_5f)
cat("ω'Y replicates point etimates?", 
    all.equal(as.numeric(omega_dml_5f$omega %*% Y),
      as.numeric(results_dml_5f[,1])
    ))

ω'Y replicates point etimates? TRUE

Check covariate balancing

We use the infrastructure of the cobalt package to plot Standardized Mean Differences where we need to flip the sign of the untreated outcome weights to make them compatible with the package framework. This is achieved by multiplying the outcome weights by \(2 \times D-1\):

threshold = 0.1

create_love_plot = function(title, index) {
  love.plot(
    D ~ X,
    weights = list(
      "2-fold" = omega_dml_2f$omega[index, ] * (2*D-1),
      "5-fold" = omega_dml_5f$omega[index, ] * (2*D-1)
    ),
    position = "bottom",
    title = title,
    thresholds = c(m = threshold),
    var.order = "unadjusted",
    binary = "std",
    abs = TRUE,
    line = TRUE,
    colors = viridis(3), # color-blind-friendly
    shapes = c("circle", "triangle", "diamond")
  )
}

# Now you can call this function for each plot:
love_plot_plr = create_love_plot("PLR", 1)
love_plot_plriv = create_love_plot("PLR-IV", 2)
love_plot_aipw = create_love_plot("AIPW", 3)
love_plot_waipw = create_love_plot("Wald-AIPW", 4)
love_plot_plr

love_plot_plriv

love_plot_aipw

love_plot_waipw

Create the combined plot that ends up in the paper as Figure 2:

figure2 = grid.arrange(
  love_plot_plr, love_plot_aipw,
  love_plot_plriv,love_plot_waipw,
  nrow = 2
)

2. Supplementary results

10-folds

Alternatively, we can also apply 10 instead of 5-fold cross-fitting:

dml_10f = dml_with_smoother(Y,D,X,Z,n_cf_folds = 10)
results_dml_10f = summary(dml_10f)

          Estimate      SE      t         p    
PLR        13530.4  1516.6 8.9216 < 2.2e-16 ***
PLR-IV     12609.8  1885.9 6.6864 2.411e-11 ***
AIPW-ATE   11298.4  1158.4 9.7537 < 2.2e-16 ***
Wald-AIPW  11080.6  1647.1 6.7271 1.826e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

omega_dml_10f = get_outcome_weights(dml_10f)
cat("ω'Y replicates point etimates?",
    all.equal(as.numeric(omega_dml_10f$omega %*% Y),
      as.numeric(results_dml_10f[,1])
    ))

ω'Y replicates point etimates? TRUE

rm(dml_10f)

We observe that it does sometimes improve over 5-fold and sometimes is worse. 5-folds seem to provide a good compromise between balancing quality and computational speed.

create_love_plot = function(title, index) {
  love.plot(
    D ~ X,
    weights = list(
      "2-fold" = omega_dml_2f$omega[index, ] * (2*D-1),
      "5-fold" = omega_dml_5f$omega[index, ] * (2*D-1),
      "10-fold" = omega_dml_10f$omega[index, ] * (2*D-1)
    ),
    position = "bottom",
    title = title,
    thresholds = c(m = threshold),
    var.order = "unadjusted",
    binary = "std",
    abs = TRUE,
    line = TRUE,
    colors = viridis(4), # color-blind-friendly
  )
}

# Now you can call this function for each plot:
love_plot_plr = create_love_plot("PLR", 1)
love_plot_plriv = create_love_plot("PLR-IV", 2)
love_plot_aipw = create_love_plot("AIPW", 3)
love_plot_waipw = create_love_plot("Wald-AIPW", 4)
love_plot_plr

love_plot_plriv

love_plot_aipw

love_plot_waipw

Weights descriptives

Finally, we investigate descriptives of the weights. Given the results of the paper it is most interesting to observe that the “Sum of weights” of PLR(-IV) and Wald-AIPW indeed are only scale-normalized. However, they all sum to values very close to one.

cat("2-fold: \n")

2-fold: 

summary(omega_dml_2f)

Weight summary for estimator PLR 
                          Control    Treated    
Minimum weight            -epsilon*  0.0001     
Maximum weight            0.0005     0.0006     
% Negative                0.0194     0.0000     
Sum largest 10%           0.2322     0.1394     
Sum of weights            0.9963     0.9963     
Sum of absolute weights   0.9985     0.9963     
* epsilon <  1e-04 

Weight summary for estimator PLR-IV 
                          Control    Treated    
Minimum weight            -0.0007    epsilon*   
Maximum weight            0.0007     0.0007     
% Negative                0.1527     0.0000     
Sum largest 10%           0.3573     0.1649     
Sum of weights            0.9953     0.9953     
Sum of absolute weights   1.8597     0.9953     
* epsilon <  1e-04 

Weight summary for estimator AIPW-ATE 
                          Control    Treated    
Minimum weight            epsilon*   epsilon*   
Maximum weight            0.0003     0.0030     
% Negative                0.0000     0.0000     
Sum largest 10%           0.1562     0.2878     
Sum of weights            1.0000     1.0000     
Sum of absolute weights   1.0000     1.0000     
* epsilon <  1e-04 

Weight summary for estimator Wald-AIPW 
                          Control    Treated    
Minimum weight            -0.0028    -epsilon*  
Maximum weight            0.0009     0.0030     
% Negative                0.1485     0.0004     
Sum largest 10%           0.3057     0.2920     
Sum of weights            0.9993     0.9993     
Sum of absolute weights   1.8643     0.9994     
* epsilon <  1e-04 

cat("\n\n5-fold: \n")

5-fold: 

summary(omega_dml_5f)

Weight summary for estimator PLR 
                          Control    Treated    
Minimum weight            -epsilon*  0.0001     
Maximum weight            0.0005     0.0006     
% Negative                0.0022     0.0000     
Sum largest 10%           0.2287     0.1388     
Sum of weights            0.9978     0.9978     
Sum of absolute weights   0.9978     0.9978     
* epsilon <  1e-04 

Weight summary for estimator PLR-IV 
                          Control    Treated    
Minimum weight            -0.0007    epsilon*   
Maximum weight            0.0007     0.0007     
% Negative                0.1524     0.0000     
Sum largest 10%           0.3569     0.1640     
Sum of weights            0.9969     0.9969     
Sum of absolute weights   1.8769     0.9969     
* epsilon <  1e-04 

Weight summary for estimator AIPW-ATE 
                          Control    Treated    
Minimum weight            epsilon*   0.0001     
Maximum weight            0.0004     0.0033     
% Negative                0.0000     0.0000     
Sum largest 10%           0.1551     0.2881     
Sum of weights            1.0000     1.0000     
Sum of absolute weights   1.0000     1.0000     
* epsilon <  1e-04 

Weight summary for estimator Wald-AIPW 
                          Control    Treated    
Minimum weight            -0.0040    0.0001     
Maximum weight            0.0011     0.0042     
% Negative                0.1486     0.0000     
Sum largest 10%           0.3058     0.2893     
Sum of weights            0.9999     0.9999     
Sum of absolute weights   1.8772     0.9999     

# cat("\n\n10-fold: \n")
# summary(omega_dml_10f)