Sensitivity Analysis for the Effect of 401(k) on Net Financial Assets

Victor Chernozhukov, Carlos Cinelli, Whitney Newey, Amit Sharma, and Vasilis Syrgkanis

Introduction

Unmeasured confounding is a pervasive issue in studies that aim to draw causal inferences from observational data. Such studies typically rely on a conditional ignorability (also known as unconfoundedness) assumption, which states that the treatment assignment is independent of potential outcomes given a set of observed covariates. This assumption, however, requires that there are no unobserved confounders influencing both the treatment and the outcome—a requirement that is often difficult to defend in practice.

The R package dml.sensemakr implements the sensitivity analysis tools developed in Chernozhukov, Cinelli, Newey, Sharma, and Syrgkanis (2026), which provides a general theory of omitted variable bias for a wide range of causal parameters estimated with machine learning. The theory applies to nonparametric models, while naturally allowing for semiparametric restrictions such as partial linearity. We show how simple plausibility judgments on the maximum explanatory power of omitted variables are sufficient to bound the magnitude of the bias, thus facilitating sensitivity analysis in otherwise complex, nonlinear models.

Since sensitivity analysis requires contextual knowledge to be properly interpreted, we illustrate the functionality of the package with a real example. Here we reproduce the results found in Section 5 of Chernozhukov et al. (2026), which estimates the average treatment effect of 401(k) eligibility on net financial assets. Further details about this application can be found in Poterba et al. (1994, 1995) and Chernozhukov et al. (2018).

The data and research question

A 401(k) plan is an employer-sponsored tax-deferred savings option that allows individuals to deduct contributions from their taxable income and accrue tax-free interest on investments within the plan. Introduced in the early 1980s as an incentive to increase individual savings for retirement, an important question in the savings literature is precisely to quantify the causal impact of 401(k) eligibility on net financial assets.

The dml.sensemakr package comes with an example dataset, pension, drawn from the 1991 Survey of Income and Program Participation (SIPP). The data can be loaded with the command data("pension").

# loads package
library(dml.sensemakr)

## loads data
data("pension")
y <- pension$net_tfa # net total financial assets
d <- pension$e401 # 401K eligibility
x <- model.matrix(~ -1 + age + inc + educ+ fsize + marr + twoearn + pira + hown, data = pension)

The outcome variable of interest is net_tfa, which measures net total financial assets. The treatment variable is e401, an indicator for whether the individual is eligible to enroll in a 401(k) plan. The vector of observed covariates consists of: age (in years), inc (income), educ (years of education), fsize (family size), marr (a binary indicator of marital status), twoearn (a two-earner household indicator), pira (an IRA participation indicator), and hown (a homeownership indicator). For further details, see ?pension.

A naive comparison of net financial assets between those with and without 401(k) eligibility suggests a positive and large impact:

mean(y[d==1]) - mean(y[d==0])
#> [1] 19559.34

This naive difference amounts to approximately $19,559. The problem, however, is that 401(k) plans can be obtained only by those individuals who work for a firm that offers such a savings option—and employment decisions are far from randomized. As an attempt to overcome this lack of random assignment, Poterba et al. (1994, 1995), and more recently Chernozhukov et al. (2018), leveraged the 1991 SIPP data to adjust for potential confounding factors between 401(k) eligibility and net financial assets. As explained in Poterba et al. (1994), at least around the time 401(k) plans initially became available, people were unlikely to make employment decisions based on whether an employer offered a 401(k) plan; instead, their main focus was on salary and other aspects of the job. Thus, as a first approximation, whether one is eligible for a 401(k) plan could be taken as ignorable once we condition on income and other covariates related to job choice.

More precisely, we consider that the decision to work for a firm that offers a 401(k) plan depends both on the observed covariates $X$, but also on unobserved firm characteristics, denoted by $U$; moreover, $X$, $U$, and $D$ are jointly affected by a set of latent factors $L$. Most importantly, note the assumption of absence of direct arrows, both from $U$ and $L$, to $Y$. Under such assumption, conditional ignorability holds adjusting for $X$ only. The story represented by this causal diagram is one way of rationalizing the identification strategy used in earlier papers.

Estimates under conditional ignorability

We now estimate the average treatment effect (ATE) of 401(k) eligibility on net financial assets under this conditional ignorability assumption. For this, we use the dml() function, which estimates causal parameters using debiased machine learning (DML) with 5-fold cross-fitting. The procedure is repeated 5 times to reduce the variance that stems from sample splitting, and estimates are combined using the median as the final estimate. The argument model = "plm" specifies a partially linear model, in which the treatment enters linearly but the nuisance functions (the conditional expectations of $Y$ and $D$ given $X$) are estimated nonparametrically using Random Forests.

# run DML
dml.401k.plm <- dml(y, d, x, model = "plm", cf.folds = 5, cf.reps = 5)
#> Debiased Machine Learning
#> 
#>  Model: Partially Linear 
#>  Target: ate (unconditional) 
#>  Cross-Fitting: 5 folds, 5 reps 
#>  ML Method: outcome (yreg0:ranger, yreg1:ranger), treatment (ranger)
#>  Tuning: dirty 
#> 
#> 
#> ====================================
#> Tuning parameters using all the data
#> ====================================
#> 
#> - Tuning Model for D.
#> -- Best Tune:
#>   mtry min.node.size splitrule
#> 1    2             5  variance
#> 
#> - Tuning Model for Y (partially linear).
#> -- Best Tune:
#>   mtry min.node.size splitrule
#> 1    2             5  variance
#> 
#> 
#> ======================================
#> Repeating 5-fold cross-fitting 5 times
#> ======================================
#> 
#> -- Rep 1 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 2 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 3 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 4 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 5 -- Folds: 1  2  3  4  5

The summary() method reports the estimated ATE, along with the $R^2$ of the outcome and treatment models:

# results under conditional ignorability
summary(dml.401k.plm)
#> 
#> Debiased Machine Learning
#> 
#>  Model: Partially Linear 
#>  Cross-Fitting: 5 folds, 5 reps 
#>  ML Method: outcome (yreg0:ranger, yreg1:ranger, R2 = 27.269%), treatment (ranger, R2 = 11.551%)
#>  Tuning: dirty 
#> 
#> Average Treatment Effect : 
#> 
#>     Estimate Std. Error t value  P(>|t|)    
#> ate     8966       1314   6.825 8.79e-12 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Note: DML estimates combined using the median method.
#> 
#> Verbal interpretation of DML procedure:
#> 
#> -- Average treatment effects were estimated using DML with 5-fold cross-fitting. In order to reduce the variance that stems from sample splitting, we repeated the procedure 5 times. Estimates are combined using the median as the final estimate, incorporating variation across experiments into the standard error as described in Chernozhukov et al. (2018). The outcome regression uses Random Forest from the R package ranger; the treatment regression uses Random Forest from the R package ranger.

As we can see, after flexibly taking into account observed confounding factors, the estimate of the effect of 401(k) eligibility on net financial assets is substantially attenuated relative to the naive difference, but it is still large, positive, and statistically significant (approximately $9,000).

Sensitivity analysis

The threat of unobserved confounding

The previous estimate requires the assumption of no unobserved confounders for unbiasedness. This assumption, however, may be violated. Employers often offer a benefit in which they “match” a proportion of an employee’s contribution to their 401(k) up to 5% of the employee’s salary. This “matched amount,” denoted by $M$, is determined by unobserved firm characteristics $U$, observed worker characteristics $X$, and by 401(k) eligibility $D$. In this model, adjustment for $X$ alone is not sufficient for control of confounding. Instead, we now need to condition both on observed covariates $X$and unobserved confounders $U$ for ignorability to hold. How strong would the omitted firm characteristics $U$ have to be in order to overturn our previous conclusions? And how plausible are the strengths revealed to be problematic? In what follows, we use our sensitivity analysis results to address these questions.

Robustness values and the main confounding scenario

A useful summary of how robust the estimate is to unmeasured confounding is the robustness value (RV). Whereas standard errors, t-values or p-values communicate how robust the estimate is to sampling errors, the idea of robustness values is to quickly communicate how robust the estimate is to systematic errors due to residual confounding. More precisely, the $\text{RV}_{\theta = 0, \alpha = 0.05}$ measures the minimal strength on both confounding factors such that the estimated confidence bound for the ATE would include zero, at the 5% significance level. Any postulated confounding scenario that does not meet this minimal criterion of strength cannot overturn the results of the original study.

We now proceed to construct a particular confounding scenario, based on the contextual details of the problem. We start with the assumption that $U$ explains as much variation in net financial assets as the total variation of the maximal matched amount of income (5%) over the period of three years (roughly the period over which the effect is measured). In the worst case scenario, this would lead to an additional 3% of total variation explained, resulting in a partial $R^2$ of outcome with omitted firm characteristics $U$ of $C_Y \approx 4\%$. Following similar reasoning, and more conservatively, we posit that omitted firm characteristics can explain an additional 2.5% of the variation in 401(k) eligibility, resulting in $C_D \approx 3\%$ for the partially linear model. Since both $C_Y \approx 4\%$ and $C_D \approx 3\%$ are below the robustness value, we immediately conclude that such a confounding scenario is not capable of bringing the lower limit of the confidence bound of the ATE to zero.

We apply the sensemakr() method to the dml object, passing the confounding scenario through the arguments cf.y (the postulated partial $R^2$ of omitted variables with the outcome) and cf.d (the postulated partial $R^2$ of omitted variables with the treatment):

# sensitivity analysis
sens.401k <- sensemakr(dml.401k.plm, cf.y = 0.04, cf.d = 0.03,)

# summary
summary(sens.401k)
#> ==== Original Analysis ====
#> 
#> Debiased Machine Learning
#> 
#>  Model: Partially Linear 
#>  Cross-Fitting: 5 folds, 5 reps 
#>  ML Method: outcome (yreg0:ranger, yreg1:ranger, R2 = 27.269%), treatment (ranger, R2 = 11.551%)
#>  Tuning: dirty 
#> 
#> Average Treatment Effect : 
#> 
#>     Estimate Std. Error t value  P(>|t|)    
#> ate     8966       1314   6.825 8.79e-12 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Note: DML estimates combined using the median method.
#> 
#> Verbal interpretation of DML procedure:
#> 
#> -- Average treatment effects were estimated using DML with 5-fold cross-fitting. In order to reduce the variance that stems from sample splitting, we repeated the procedure 5 times. Estimates are combined using the median as the final estimate, incorporating variation across experiments into the standard error as described in Chernozhukov et al. (2018). The outcome regression uses Random Forest from the R package ranger; the treatment regression uses Random Forest from the R package ranger.
#> 
#> ==== Sensitivity Analysis ====
#> 
#> Null hypothesis: theta = 0 
#> Signif. level: alpha = 0.05 
#> 
#> Robustness Values:
#>     RV (%) RVa (%)
#> ate 7.2641  5.4829
#> 
#> Verbal interpretation of robustness values:
#> 
#> -- Robustness Value for the Bound (RV): omitted variables that explain more than RV% of the residual variation both of the treatment (cf.d) and of the outcome (cf.y) are sufficiently strong to make the estimated bounds include 0. Conversely, omitted variables that do not explain more than RV% of the residual variation of both the treatment and the outcome are not sufficiently strong to do so.
#> 
#> -- Robustness Value for the Confidence Bound (RVa): omitted variables that explain more than RV% of the residual variation both of the treatment (cf.d) and of the outcome (cf.y) are sufficiently strong to make the confidence bounds include 0, at the  significance level of alpha = 0.05. Conversely, omitted variables that do not explain more than RV% of the residual variation of both the treatment and the outcome are not sufficiently strong to do so.
#> 
#> Bounds on omitted variable bias (postulated scenario):
#> 
#>  target          bound.label cf.y cf.d rho2 theta.minus theta.plus      lwr
#>     ate Confounding Scenario 0.04 0.03    1    4785.561    13147.2 2593.483
#>       upr
#>  15396.16
#> 
#> -- Bounds on each target parameter under latent variables with the postulated sensitivity parameters cf.y, cf.d and rho2. theta.minus and theta.plus are the point bounds; lwr and upr are one-sided 95% confidence bounds treating those parameters as fixed.

# contour plots
oldpar <- par(mfrow = c(1,2))
plot(sens.401k, which.bound = "lwr", bound.label = "Max match (3 years)", col.contour = "blue")
plot(sens.401k, which.bound = "upr", bound.label = "Max match (3 years)", col.contour = "blue")

par(oldpar)

The summary() method reports the robustness value, the bias, the bounds, and the confidence bounds for the ATE under the postulated scenario. The plot() method produces contour plots for the lower and upper limits of the confidence bounds.

Sensitivity contour plots

Researchers may wish to refine their analysis by visually exploring the whole range of possible estimates that confounders with different strengths could cause. The contour plots above serve this purpose. The horizontal axis describes the fraction of residual variation of the treatment explained by unobserved confounders, whereas the vertical axis describes the share of residual variation of the outcome explained by unobserved confounders. The contour lines show the lower (left panel) and upper (right panel) limits of the confidence bounds for the ATE, given a pair of hypothesized values of partial $R^2$. The red dashed line marks the zero contour—the critical threshold beyond which the estimate changes sign. Note $\text{RV}_{\theta = 0, \alpha = 0.05}$ is simply the point where the zero contour crosses the 45-degree line. The labeled point Max match (3 years) shows the bounds on the partial $R^2$ as previously discussed, confirming that this confounding scenario is not sufficient to overturn the original conclusions. Any other scenario can also be contemplated by simply picking different pairs of $R^2$ values on the plot.

Benchmarking with observed covariates

Another approach to construct confounding scenarios is to use observed covariates to bound the plausible strength of unobserved covariates. For instance, in our empirical example, we know that employment decisions are largely driven by salary considerations. Similarly, salary is clearly an important determinant of net financial assets. One could therefore argue that it is implausible to imagine other latent firm characteristics that would be even a fraction as strong as the observed income of individuals, in terms of explanatory power in predicting 401(k) eligibility and net financial assets. Whenever such claims of relative importance can be made, they can be used to set plausible bounds on the strength of unmeasured confounding. The function dml_benchmark() automates this comparison:

# benchmarks
bench.plm <- dml_benchmark(dml.401k.plm, benchmark_covariates = c("inc", "pira", "twoearn"))
#> 
#> === Computing benchmarks using covariate: inc  ===
#> 
#> Debiased Machine Learning
#> 
#>  Model: Partially Linear 
#>  Target: ate (unconditional) 
#>  Cross-Fitting: 5 folds, 5 reps 
#>  ML Method: outcome (yreg0:ranger, yreg1:ranger), treatment (ranger)
#>  Tuning: dirty 
#> 
#> 
#> ====================================
#> Tuning parameters using all the data
#> ====================================
#> 
#> - Tuning Model for D.
#> -- Best Tune:
#>   mtry min.node.size splitrule
#> 1    2             5  variance
#> 
#> - Tuning Model for Y (partially linear).
#> -- Best Tune:
#>   mtry min.node.size splitrule
#> 1    2             5  variance
#> 
#> 
#> ======================================
#> Repeating 5-fold cross-fitting 5 times
#> ======================================
#> 
#> -- Rep 1 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 2 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 3 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 4 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 5 -- Folds: 1  2  3  4  5  
#> 
#> 
#> === Computing benchmarks using covariate: pira  ===
#> 
#> Debiased Machine Learning
#> 
#>  Model: Partially Linear 
#>  Target: ate (unconditional) 
#>  Cross-Fitting: 5 folds, 5 reps 
#>  ML Method: outcome (yreg0:ranger, yreg1:ranger), treatment (ranger)
#>  Tuning: dirty 
#> 
#> 
#> ====================================
#> Tuning parameters using all the data
#> ====================================
#> 
#> - Tuning Model for D.
#> -- Best Tune:
#>   mtry min.node.size splitrule
#> 1    2             5  variance
#> 
#> - Tuning Model for Y (partially linear).
#> -- Best Tune:
#>   mtry min.node.size splitrule
#> 1    2             5  variance
#> 
#> 
#> ======================================
#> Repeating 5-fold cross-fitting 5 times
#> ======================================
#> 
#> -- Rep 1 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 2 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 3 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 4 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 5 -- Folds: 1  2  3  4  5  
#> 
#> 
#> === Computing benchmarks using covariate: twoearn  ===
#> 
#> Debiased Machine Learning
#> 
#>  Model: Partially Linear 
#>  Target: ate (unconditional) 
#>  Cross-Fitting: 5 folds, 5 reps 
#>  ML Method: outcome (yreg0:ranger, yreg1:ranger), treatment (ranger)
#>  Tuning: dirty 
#> 
#> 
#> ====================================
#> Tuning parameters using all the data
#> ====================================
#> 
#> - Tuning Model for D.
#> -- Best Tune:
#>   mtry min.node.size splitrule
#> 1    2             5  variance
#> 
#> - Tuning Model for Y (partially linear).
#> -- Best Tune:
#>   mtry min.node.size splitrule
#> 1    2             5  variance
#> 
#> 
#> ======================================
#> Repeating 5-fold cross-fitting 5 times
#> ======================================
#> 
#> -- Rep 1 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 2 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 3 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 4 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 5 -- Folds: 1  2  3  4  5
bench.plm
#> Covariate benchmarks:
#> 
#>         gain.Y se.gain.Y gain.D se.gain.D     rho se.rho      delta se.delta
#> inc     0.1512    0.0485 0.0481    0.0051 -0.3365 0.0754 -3370.8812 490.2191
#> pira    0.0515    0.0227 0.0018    0.0026 -0.1327 0.4074  -265.6566 412.6432
#> twoearn 0.0376    0.0164 0.0072    0.0034  0.4850 0.2200   816.1287 312.6554
#> 
#> Combined across cross-fitting repetitions using the median method.

The output shows the bounds on the strength of the latent variable $U$ if it were as strong as (i) income (inc), (ii) whether a worker has an individual retirement account (pira), and (iii) whether the worker’s family has a two-earner status (twoearn). Note that, apart from income, latent variables as strong as these covariates would result in a weaker confounding scenario than the one we have previously considered (Max match). As for income, the worst-case bound indicates that omitted firm characteristics as important as income could be sufficiently strong to overturn the original results. However, one could argue such a scenario to be implausible, as it is hard to imagine latent firm characteristics that would explain more variation in job choice than income itself. A more realistic, but still conservative, scenario is thus provided by the benchmark point 1/4 x Income, which shows the bound on the strength of $U$ if it were 25% as strong as income in predicting treatment and outcome variation.

Nonparametric model

The partially linear model imposes a constant treatment effect. With a binary treatment, however, a fully nonparametric model allows for more flexible estimation and, importantly, for exploring heterogeneous treatment effects. We now stratify by income quartiles and estimate the group average treatment effects (GATEs), in addition to the overall ATE. This is done by passing a groups argument to dml(), along with model = "npm" for the nonparametric model:

## compute income quartiles
g1 <- cut(x[,"inc"], quantile(x[,"inc"], c(0, 0.25,.5,.75,1), na.rm = TRUE),
          labels = c("q1", "q2", "q3", "q4"), include.lowest = TRUE)

## Nonparametric model
dml.401k.npm <- dml(y, d, x, groups = g1, model = "npm", cf.folds = 5, cf.reps = 5)
#> Debiased Machine Learning
#> 
#>  Model: Nonparametric 
#>  Target: ate (unconditional) 
#>  Cross-Fitting: 5 folds, 5 reps 
#>  ML Method: outcome (yreg0:ranger, yreg1:ranger), treatment (ranger)
#>  Tuning: dirty 
#> 
#> 
#> ====================================
#> Tuning parameters using all the data
#> ====================================
#> 
#> - Tuning Model for D.
#> -- Best Tune:
#>   mtry min.node.size splitrule
#> 1    2             5  variance
#> 
#> - Tuning Model for Y (non-parametric).
#> -- Best Tune:
#>   mtry min.node.size splitrule
#> 1    2             5  variance
#>   mtry min.node.size splitrule
#> 1    2             5  variance
#> 
#> 
#> ======================================
#> Repeating 5-fold cross-fitting 5 times
#> ======================================
#> 
#> -- Rep 1 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 2 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 3 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 4 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 5 -- Folds: 1  2  3  4  5

The summary() method now reports both the overall ATE and the group-specific treatment effects:

# results under conditional ignorability
summary(dml.401k.npm)
#> 
#> Debiased Machine Learning
#> 
#>  Model: Nonparametric 
#>  Cross-Fitting: 5 folds, 5 reps 
#>  ML Method: outcome (yreg0:ranger, yreg1:ranger, R2 = 26.641%), treatment (ranger, R2 = 11.432%)
#>  Tuning: dirty 
#> 
#> Average Treatment Effect : 
#> 
#>     Estimate Std. Error t value  P(>|t|)    
#> ate     7904       1174   6.734 1.65e-11 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Group Average Treatment Effect: 
#> 
#>          Estimate Std. Error t value  P(>|t|)    
#> g.ate.q1   4515.8      799.2   5.651 1.60e-08 ***
#> g.ate.q2   2546.1     1292.5   1.970 0.048851 *  
#> g.ate.q3   6802.5     1824.1   3.729 0.000192 ***
#> g.ate.q4  17681.5     4015.9   4.403 1.07e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Note: DML estimates combined using the median method.
#> 
#> Verbal interpretation of DML procedure:
#> 
#> -- Average treatment effects were estimated using DML with 5-fold cross-fitting. In order to reduce the variance that stems from sample splitting, we repeated the procedure 5 times. Estimates are combined using the median as the final estimate, incorporating variation across experiments into the standard error as described in Chernozhukov et al. (2018). The outcome regression uses Random Forest from the R package ranger; the treatment regression uses Random Forest from the R package ranger.

The ATE varies substantially across groups, with effects ranging from approximately $4,000 (first quartile) to almost $18,000 (last quartile). We can visualize the heterogeneous treatment effects by plotting the point estimates along with their confidence intervals for each income quartile:

# load ggplot2
library(ggplot2)

# coefficient plot under conditional ignorability
group.names <- paste0("g.ate.q",1:4)
df   <- data.frame(groups = 1:4, estimate = coef(dml.401k.npm)[group.names])
cis  <- confint(dml.401k.npm)[group.names, ]
cis  <- setNames(as.data.frame(cis), c("lwr.ci", "upr.ci"))
df   <- cbind(df, cis)
ggplot(df, aes(x = groups, y = estimate)) + geom_line() +
  geom_ribbon(aes(ymin = lwr.ci, ymax = upr.ci), alpha = 0.1, col = "blue", fill = "blue") +
  theme_bw() +
  xlab("Income Groups by Quartiles") +
  ylab("ATE")

Sensitivity analysis for the nonparametric model

We now repeat the sensitivity analysis using the same confounding scenario ($C_Y = 4\%$, $C_D = 3\%$). The interpretation of the contours is the same as before, with the main difference being that the horizontal axis now describes gains in precision rather than gains in variance explained (see Remark 3 of Chernozhukov et al., 2026). We adopt the same scenario for the NPM, with the understanding that $C_D \approx 3\%$ now corresponds to gains in precision.

# sensitivity analysis
sens.401k.npm <- sensemakr(dml.401k.npm, cf.y = 0.04, cf.d = 0.03)

# summary
summary(sens.401k.npm)
#> ==== Original Analysis ====
#> 
#> Debiased Machine Learning
#> 
#>  Model: Nonparametric 
#>  Cross-Fitting: 5 folds, 5 reps 
#>  ML Method: outcome (yreg0:ranger, yreg1:ranger, R2 = 26.641%), treatment (ranger, R2 = 11.432%)
#>  Tuning: dirty 
#> 
#> Average Treatment Effect : 
#> 
#>     Estimate Std. Error t value  P(>|t|)    
#> ate     7904       1174   6.734 1.65e-11 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Group Average Treatment Effect: 
#> 
#>          Estimate Std. Error t value  P(>|t|)    
#> g.ate.q1   4515.8      799.2   5.651 1.60e-08 ***
#> g.ate.q2   2546.1     1292.5   1.970 0.048851 *  
#> g.ate.q3   6802.5     1824.1   3.729 0.000192 ***
#> g.ate.q4  17681.5     4015.9   4.403 1.07e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Note: DML estimates combined using the median method.
#> 
#> Verbal interpretation of DML procedure:
#> 
#> -- Average treatment effects were estimated using DML with 5-fold cross-fitting. In order to reduce the variance that stems from sample splitting, we repeated the procedure 5 times. Estimates are combined using the median as the final estimate, incorporating variation across experiments into the standard error as described in Chernozhukov et al. (2018). The outcome regression uses Random Forest from the R package ranger; the treatment regression uses Random Forest from the R package ranger.
#> 
#> ==== Sensitivity Analysis ====
#> 
#> Null hypothesis: theta = 0 
#> Signif. level: alpha = 0.05 
#> 
#> Robustness Values:
#>           RV (%) RVa (%)
#> ate       6.0267  4.5051
#> g.ate.q1 10.6673  7.6517
#> g.ate.q2  4.1205  0.5864
#> g.ate.q3  6.8353  3.6177
#> g.ate.q4  8.8498  5.5311
#> 
#> Verbal interpretation of robustness values:
#> 
#> -- Robustness Value for the Bound (RV): omitted variables that explain more than RV% of the residual variation of the outcome (cf.y) and generate an additional RV% of variation on the Riesz Representer (cf.d) are sufficiently strong to make the estimated bounds include 0. Conversely, omitted variables that do not explain more than RV% of the residual variation of the outcome nor generate an additional RV% of variation on the Riesz Representer are not sufficiently strong to do so.
#> 
#> -- Robustness Value for the Confidence Bound (RVa): omitted variables that explain more than RV% of the residual variation of the outcome (cf.y) and generate an additional RV% of variation on the Riesz Representer (cf.d) are sufficiently strong to make the confidence bounds include 0, at the  significance level of alpha = 0.05. Conversely, omitted variables that do not explain more than RV% of the residual variation of the outcome nor generate an additional RV% of variation on the Riesz Representer are not sufficiently strong to do so. 
#> 
#>  The interpretation of sensitivity parameters can be further refined for each target quantity. See more below.
#> 
#> Bounds on omitted variable bias (postulated scenario):
#> 
#>    target          bound.label cf.y cf.d rho2 theta.minus theta.plus        lwr
#>       ate Confounding Scenario 0.04 0.03    1   3432.4239  12376.117  1442.4523
#>  g.ate.q1 Confounding Scenario 0.04 0.03    1   3108.4768   5949.525  1797.6232
#>  g.ate.q2 Confounding Scenario 0.04 0.03    1    418.0099   4674.275 -2151.5788
#>  g.ate.q3 Confounding Scenario 0.04 0.03    1   3423.8196  10181.098   180.6957
#>  g.ate.q4 Confounding Scenario 0.04 0.03    1  10888.0291  24475.065  4227.4794
#>        upr
#>  14392.168
#>   7286.242
#>   6552.709
#>  13190.199
#>  31192.240
#> 
#> -- Bounds on each target parameter under latent variables with the postulated sensitivity parameters cf.y, cf.d and rho2. theta.minus and theta.plus are the point bounds; lwr and upr are one-sided 95% confidence bounds treating those parameters as fixed.
#> 
#> 
#> Interpretation of sensitivity parameters:
#> 
#> -- cf.y: percentage of the residual variation of the outcome explained by latent variables.
#> -- cf.d: percentage gains in the variation of the Riesz Representer generated by latent variables:
#>    ATE: cf.d measures the percentage gains in the average precision on the treatment regression.
#> -- For Group Average Treatment Effects (GATE), parameters are conditional on the relevant group.

# sensitivity contour plots
oldpar <- par(mfrow = c(1,2))
plot(sens.401k.npm, which.bound = "lwr", bound.label = "Max match (3 years)", col.contour = "blue")
plot(sens.401k.npm, which.bound = "upr", bound.label = "Max match (3 years)", col.contour = "blue")

par(oldpar)

Confidence bounds for group-wise ATEs can also be computed. The confidence_bounds() function returns the lower and upper bounds on the treatment effect under the confounding scenario; setting level = 0 gives the point-wise bounds (no statistical uncertainty, only the bias from confounding), whereas level = 0.95 gives the confidence bounds (accounting for both bias and sampling uncertainty). The results are shown in the figure below, where the red band shows the bounds and the blue band shows the confidence bounds.

# confidence bounds plot
bds   <- confidence_bounds(dml.401k.npm, cf.y = 0.04, cf.d = 0.03, level = 0)
bds   <- setNames(as.data.frame(bds), c("lwr.bound", "upr.bound"))
cbds  <- confidence_bounds(dml.401k.npm, cf.y = 0.04, cf.d = 0.03, level = .95)
cbds  <- setNames(as.data.frame(cbds), c("lwr.cbound", "upr.cbound"))
df2   <- cbind(df, bds[-1,], cbds[-1, ])
ggplot(df2, aes(x = groups, y = estimate)) + geom_line() +
  geom_ribbon(aes(ymin = lwr.bound, ymax = upr.bound),   alpha = 0.1, col = "red", fill = "red") +
  geom_ribbon(aes(ymin = lwr.cbound, ymax = upr.cbound), alpha = 0.1, col = "blue", fill = "blue") +
  theme_bw() +
  xlab("Income Groups by Quartiles") +
  ylab("ATE")

Note how the bounds are still largely positive, with only a small excursion into the negative in the case of the second quartile group. These results suggest that the main qualitative findings reported in earlier studies are relatively robust to plausible violations of unconfoundedness, such as the one specified by our confounding scenario.

Finally, we provide benchmarks for the nonparametric model:

# benchmarks
bench.npm <- dml_benchmark(dml.401k.npm, benchmark_covariates = c("inc", "pira", "twoearn"))
#> 
#> === Computing benchmarks using covariate: inc  ===
#> 
#> Debiased Machine Learning
#> 
#>  Model: Nonparametric 
#>  Target: ate (unconditional) 
#>  Cross-Fitting: 5 folds, 5 reps 
#>  ML Method: outcome (yreg0:ranger, yreg1:ranger), treatment (ranger)
#>  Tuning: dirty 
#> 
#> 
#> ====================================
#> Tuning parameters using all the data
#> ====================================
#> 
#> - Tuning Model for D.
#> -- Best Tune:
#>   mtry min.node.size splitrule
#> 1    2             5  variance
#> 
#> - Tuning Model for Y (non-parametric).
#> -- Best Tune:
#>   mtry min.node.size splitrule
#> 1    2             5  variance
#>   mtry min.node.size splitrule
#> 1    2             5  variance
#> 
#> 
#> ======================================
#> Repeating 5-fold cross-fitting 5 times
#> ======================================
#> 
#> -- Rep 1 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 2 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 3 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 4 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 5 -- Folds: 1  2  3  4  5  
#> 
#> 
#> === Computing benchmarks using covariate: pira  ===
#> 
#> Debiased Machine Learning
#> 
#>  Model: Nonparametric 
#>  Target: ate (unconditional) 
#>  Cross-Fitting: 5 folds, 5 reps 
#>  ML Method: outcome (yreg0:ranger, yreg1:ranger), treatment (ranger)
#>  Tuning: dirty 
#> 
#> 
#> ====================================
#> Tuning parameters using all the data
#> ====================================
#> 
#> - Tuning Model for D.
#> -- Best Tune:
#>   mtry min.node.size splitrule
#> 1    2             5  variance
#> 
#> - Tuning Model for Y (non-parametric).
#> -- Best Tune:
#>   mtry min.node.size splitrule
#> 1    2             5  variance
#>   mtry min.node.size splitrule
#> 1    2             5  variance
#> 
#> 
#> ======================================
#> Repeating 5-fold cross-fitting 5 times
#> ======================================
#> 
#> -- Rep 1 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 2 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 3 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 4 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 5 -- Folds: 1  2  3  4  5  
#> 
#> 
#> === Computing benchmarks using covariate: twoearn  ===
#> 
#> Debiased Machine Learning
#> 
#>  Model: Nonparametric 
#>  Target: ate (unconditional) 
#>  Cross-Fitting: 5 folds, 5 reps 
#>  ML Method: outcome (yreg0:ranger, yreg1:ranger), treatment (ranger)
#>  Tuning: dirty 
#> 
#> 
#> ====================================
#> Tuning parameters using all the data
#> ====================================
#> 
#> - Tuning Model for D.
#> -- Best Tune:
#>   mtry min.node.size splitrule
#> 1    2             5  variance
#> 
#> - Tuning Model for Y (non-parametric).
#> -- Best Tune:
#>   mtry min.node.size splitrule
#> 1    2             5  variance
#>   mtry min.node.size splitrule
#> 1    2             5  variance
#> 
#> 
#> ======================================
#> Repeating 5-fold cross-fitting 5 times
#> ======================================
#> 
#> -- Rep 1 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 2 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 3 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 4 -- Folds: 1  2  3  4  5  
#> 
#> -- Rep 5 -- Folds: 1  2  3  4  5
bench.npm
#> Covariate benchmarks:
#> 
#>         gain.Y se.gain.Y gain.D se.gain.D     rho se.rho      delta se.delta
#> inc     0.1445    0.0435 0.1338    0.0186 -0.2267 0.0487 -3731.2221 468.2531
#> pira    0.0363    0.0165 0.0012    0.0167 -0.0486 0.1682  -640.6138 369.1035
#> twoearn 0.0067    0.0120 0.0000    0.0140  0.0000 0.0000   581.0005 278.4045
#> 
#> Combined across cross-fitting repetitions using the median method.

References

Chernozhukov, V., Cinelli, C., Newey, W., Sharma, A. and Syrgkanis, V. (2026), “Long Story Short: Omitted Variable Bias in Causal Machine Learning.” Review of Economics and Statistics. ( link )

Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W. and Robins, J. (2018), “Double/debiased machine learning for treatment and structural parameters.” The Econometrics Journal, 21(1), C1–C68.

Cinelli, C. and Hazlett, C. (2020), “Making Sense of Sensitivity: Extending Omitted Variable Bias.” Journal of the Royal Statistical Society, Series B (Statistical Methodology), 82(1), 39–67. ( link )

Imbens, G. (2003), “Sensitivity to Exogeneity Assumptions in Program Evaluation.” American Economic Review, 93(2), 126–132.

Poterba, J., Venti, S. and Wise, D. (1994), “401(k) Plans and Tax-Deferred Saving.” In D. Wise (Ed.), Studies in the Economics of Aging, pp. 105–142, University of Chicago Press.

Poterba, J., Venti, S. and Wise, D. (1995), “Do 401(k) Contributions Crowd Out Other Personal Saving?” Journal of Public Economics, 58(1), 1–32.