Causal ML: Double Selection and Partially Linear Double ML

Application notebook

Goals:

• Illustrate why naive implementations of model selection using Lasso are problematic

Introducing the data

The Application Notebook builds on the dataset that is kindly provided in the hdm package. The data was used in Chernozhukov and Hansen (2004). Their paper investigates the effect of participation in the employer-sponsored 401(k) retirement savings plan (p401) on net assets (net_tfa). Since then the data was used to showcase many new methods. It is not the most comprehensive dataset with basically ten covariates/regressors/predictors:

• age: age

• db: defined benefit pension

• educ: education (in years)

• fsize: family size

• hown: home owner

• inc: income (in US $)

• male: male

• marr: married

• pira: participation in individual retirement account (IRA)

• twoearn: two earners

However, it is publicly available and the relatively few covariates ensure that the programs do not run too long.

# To install the causalDML package uncomment the following two lines
# library(devtools)
# install_github(repo="MCKnaus/causalDML")

# Load the packages required for later
library(hdm)
library(tidyverse)
library(causalDML)
library(grf)
library(estimatr)

set.seed(1234) # for replicability

data(pension)
# Outcome
Y = pension$net_tfa
# Treatment
W = pension$p401
# Create main effects matrix
X = model.matrix(~ 0 + age + db + educ + fsize + hown + inc + male + marr + pira + twoearn, data = pension)

Hand-coded Double Selection

To understand the procedure of Double Selection, we proceed step by step using only the main effects for simplicity:

1. Select variables in the outcome regression without the treatment:

# Select variables in outcome regression
sel_y = rlasso(X,Y)
# Which variables are selected?
which(sel_y$beta != 0)

    age   fsize    hown     inc    pira twoearn 
      1       4       5       6       9      10 

2. Select variables in the treatment regression:

# Select variables in treatment regression
sel_w = rlasso(X,W)
which(sel_w$beta != 0)

     db    hown     inc    pira twoearn 
      2       5       6       9      10 

Note that variable db is now selected which was not selected in step one.

3. Use the union of the in total seven selected variables to run a standard OLS regression with robust standard errors:

# Double selection
X_sel_union = X[,sel_y$beta != 0 | sel_w$beta != 0]
ds_hand = lm_robust(Y ~ W + X_sel_union)
summary(ds_hand)

Call:
lm_robust(formula = Y ~ W + X_sel_union)

Standard error type:  HC2 

Coefficients:
                     Estimate Std. Error  t value  Pr(>|t|)   CI Lower  CI Upper   DF
(Intercept)        -4.268e+04  3298.4608 -12.9386 5.530e-38 -4.914e+04 -36211.86 9906
W                   1.158e+04  1811.3259   6.3953 1.675e-10  8.033e+03  15134.49 9906
X_sel_unionage      6.733e+02    58.9094  11.4298 4.579e-30  5.578e+02    788.80 9906
X_sel_uniondb      -5.465e+03  1359.2053  -4.0208 5.844e-05 -8.129e+03  -2800.73 9906
X_sel_unionfsize   -6.872e+02   301.2848  -2.2810 2.257e-02 -1.278e+03    -96.65 9906
X_sel_unionhown     9.267e+02   950.5966   0.9749 3.296e-01 -9.366e+02   2790.08 9906
X_sel_unioninc      8.904e-01     0.1019   8.7390 2.733e-18  6.907e-01      1.09 9906
X_sel_unionpira     2.849e+04  1819.6634  15.6541 1.400e-54  2.492e+04  32052.14 9906
X_sel_uniontwoearn -1.866e+04  2315.0448  -8.0607 8.465e-16 -2.320e+04 -14122.83 9906

Multiple R-squared:  0.2346 ,   Adjusted R-squared:  0.234 
F-statistic:   135 on 8 and 9906 DF,  p-value: < 2.2e-16

Double Selection with hdm package

In practice we want to have one function that does everything at once. This is the rlassoEffect command of the hdm package.

ds1 = rlassoEffect(X,Y,W)
summary(ds1)

[1] "Estimates and significance testing of the effect of target variables"
   Estimate. Std. Error t value Pr(>|t|)    
d1     11584       1809   6.405 1.51e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

It produces the same point estimate as the hand-coded version:

all.equal(as.numeric(ds1$alpha),as.numeric(ds_hand$coefficients[2]))

[1] TRUE

Only the standard error differs slightly because of different defaults. If you do not like this, applying se_type = "HC1 in the lm_robust() function replicates the standard error of rlassoEffect().

More flexible dictionaries

We can check whether more flexible covariate matrices provide different results:

• X2 with 88 variables: Second order polynomials of the continuous variables age, education and income as well as second order interactions of all variables

• X3 with 567 variables: Third order polynomials of the continuous variables age, education and income as well as third order interactions of all variables

X2 = model.matrix(~ 0 + (fsize + marr + twoearn + db + pira + hown + male +
                           poly(age,2) + poly(educ,2) + poly(inc,2))^2, data = pension)
X3 = model.matrix(~ 0 + (fsize + marr + twoearn + db + pira + hown + male +
                           poly(age,3) + poly(educ,3) + poly(inc,3))^3, data = pension)

Indeed, the effects are by more than $2000 Dollars higher, but going from two to three order terms makes basically no difference:

ds2 = rlassoEffect(X2,Y,W)
summary(ds2)

[1] "Estimates and significance testing of the effect of target variables"
   Estimate. Std. Error t value Pr(>|t|)    
d1     13893       1484    9.36   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

ds3 = rlassoEffect(X3,Y,W)
summary(ds3)

[1] "Estimates and significance testing of the effect of target variables"
   Estimate. Std. Error t value Pr(>|t|)    
d1     13839       1491   9.281   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Hand-coded Double ML for partially linear model

If we are not willing to assume a linear model and use for example random forest to estimate the nuisance parameters of a partially linear model, we need to predict the nuisance parameters out-of-sample. The easiest way to do this is via two-fold cross-fitting:

• Split the sample in two random subsamples, S1 and S2

• Form prediction models in S1, use it to predict in S2

• Form prediction models in S2, use it to predict in S1

• Run residual-on-residual regression with the combined predictions

# Initialize nuisance vectors
n = length(Y)
mhat = ehat = rep(NA,n)
# Draw random indices for sample 1
index_s1 = sample(1:n,n/2)
# Create S1
x1 = X[index_s1,]
w1 = W[index_s1]
y1 = Y[index_s1]
# Create sample 2 with those not in S1
x2 = X[-index_s1,]
w2 = W[-index_s1]
y2 = Y[-index_s1]
# Model in S1, predict in S2
rf = regression_forest(x1,w1)
ehat[-index_s1] = predict(rf,newdata=x2)$predictions
rf = regression_forest(x1,y1)
mhat[-index_s1] = predict(rf,newdata=x2)$predictions
# Model in S2, predict in S1
rf = regression_forest(x2,w2)
ehat[index_s1] = predict(rf,newdata=x1)$predictions
rf = regression_forest(x2,y2)
mhat[index_s1] = predict(rf,newdata=x1)$predictions
# RORR
res_y = Y-mhat
res_w = W-ehat
pl_2f = lm_robust(res_y ~ 0+res_w)
summary(pl_2f)

Call:
lm_robust(formula = res_y ~ 0 + res_w)

Standard error type:  HC2 

Coefficients:
      Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper   DF
res_w    13944       1517   9.191 4.68e-20    10970    16918 9914

Multiple R-squared:  0.01129 ,  Adjusted R-squared:  0.01119 
F-statistic: 84.47 on 1 and 9914 DF,  p-value: < 2.2e-16

Double ML for partially linear model with causalDML package

2-fold cross-fitting is easy to implement by hand but especially in small sample sizes, using only 50% of the data to estimate the nuisance parameters might lead to unstable predictions.

Thus, we use the DML_partial_linear function of the causalDML package to run 5-fold cross-fitting. This package requires to create the methods that we use because it allows for ensemble methods (for a more detailed intro see the GitHub page). For now, we focus again on the random forest.

With 5-fold cross-fitting, the program splits the sample in 5 folds and uses 4 folds (80% of the data) to predict the left out fold (20% of the data). It iterates such that every fold is left out once.

# 5-fold cross-fitting with causalDML package
# Create learner
forest = create_method("forest_grf",args=list(tune.parameters = "all"))
# Run partially linear model
pl_5f = DML_partial_linear(Y,W,X,ml_w=list(forest),ml_y=list(forest),cf=5)
summary(pl_5f)

     Coefficient      SE      t         p    
[1,]     13756.4  1514.7 9.0818 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Comparison of results

We can now compare all the different methods. Besides Double Selection with only the main effects all methods agree on an effect of 401(k) participation of wealth of about $14k:

# Collect the results
Coefficient = c(ds1$alpha,ds2$alpha,ds3$alpha,pl_2f$coefficients,pl_5f$result[1])
se = c(ds1$se,ds2$se,ds3$se,pl_2f$std.error,pl_5f$result[2])
data.frame(Coefficient,se,
                Method = c("DS1","DS2","DS3","PL 2-fold","PL 5-fold"),
                cil = Coefficient - 1.96*se,
                ciu = Coefficient + 1.96*se)  %>% 
  ggplot(aes(x=Method,y=Coefficient,ymin=cil,ymax=ciu)) + geom_point(size=2.5) + geom_errorbar(width=0.15)  +
  geom_hline(yintercept=0)