Debiased Machine Learning — dml • dml.sensemakr

Estimates a target parameter of interest, such as an average treatment effect (ATE), using Debiased Machine #Learning (DML).

The function dml_gate is a convenience function that adds groups to a dml object after the model is fit.

dml(
  y,
  d,
  x,
  model = c("plm", "npm"),
  target = "ate",
  groups = NULL,
  cf.folds = 5,
  cf.reps = 1,
  cf.seed = NULL,
  ps.trim = 0.01,
  reg = "ranger",
  yreg = reg,
  dreg = reg,
  dirty.tuning = TRUE,
  save.models = FALSE,
  y.class = FALSE,
  d.class = FALSE,
  conditional = NULL,
  centered = FALSE,
  verbose = TRUE,
  warnings = FALSE
)

dml_gate(dml.fit, groups, ...)

Arguments

y: numeric vector with the outcome.
d: numeric vector with the treatment. If the treatment is binary, it needs to be encoded as: zero = absence of treatment, one = presence of treatment.
x: numeric vector or matrix with covariates. We suggest constructing x using model.matrix.
model: specifies the model. Current available options are plm for a partially linear model, and npm for a fully non-parametric model.
target: specifies the target causal quantity of interest. Available options are ate (ATE - average treatment effect), att (ATT - average treatment effect on the treated), and atu (ATU - average treatment effect on the untreated). Note that for the partially linear model with a continuous treatment the ATE also equals the average causal derivative (ACD). For the nonparametric model, these are only available for binary treatments.
groups: a factor or numeric vector indicating group membership. Groups must be a deterministic function of x.
cf.folds: number of cross-fitting folds. Default is 5.
cf.reps: number of cross-fitting repetitions. Default is 1.
cf.seed: optional integer. A random seed for reproducibility of fold assignments.
ps.trim: trims propensity scores lower than ps.trim and greater than 1-ps.trim, in order to obtain more stable estimates. Alternatively, a named list with elements lower and upper specifying the lower and upper bounds for trimming. The same bounds are used for main and group-specific estimates, including groups added with dml_gate(). This is only relevant for the case of a binary treatment and when model = "npm".
reg: details of the machine learning method to be used for estimating the nuisance parameters (e.g, regression functions of the treatment and the outcome). Currently, this should be specified using the same arguments as caret's train function. The default is random forest using ranger. The default method is fast and usually works well for many applications.
yreg: same as reg, but specifies arguments for the outcome regression alone. Default is the same value of reg. Alternatively, a named list with elements yreg0 and yreg1 specifying separate methods for each.
dreg: same as reg, but specifies arguments for the treatment regression alone. Default is the same value of reg.
dirty.tuning: should the tuning of the machine learning method happen within each cross-fit fold ("clean"), or using all the data ("dirty")? Default is dirty tuning (dirty.tuning = TRUE). As long as the number of choices for the tuning parameters is not too big, dirty tuning is faster and should not affect the asymptotic guarantees of DML.
save.models: should the fitted models of each iteration be saved? Default is FALSE. Note that setting this to true could end up using a lot of memory.
y.class: when y is binary, should the outcome regression be treated as a classification problem? Default is FALSE. Note that for DML we need the class probabilities, and regression gives us that. If you change to classification, you need to make sure the method outputs class probabilities.
d.class: when d is binary, should the treatment regression be treated as a classification problem? Default is FALSE. Note that for DML we need the class probabilities, and regression gives us that. If you change to classification, you need to make sure the method outputs class probabilities.
conditional: logical or NULL. If TRUE, estimates the conditional version of the target, in which the outcome regression is fit on a single treatment arm and used to impute the counterfactual mean for the other arm. Only supported for model = "npm" with a single target of "att" or "atu". When NULL (default), it is set to TRUE automatically for those single-target npm cases and FALSE otherwise. The conditional parameterization is the recommended one (Wang et al., 2026); conditional = FALSE gives the unconditional parameterization of Chernozhukov et al. (2026), retained so those results can be reproduced.
centered: logical. Controls the parameterization of the Riesz-representer imbalance nu2.s in the conditional model. When FALSE (default), nu2.s is the full second moment of the Riesz representer, \(\chi^2 + 1\) (uncentered); when TRUE, it is the centered version, the \(\chi^2\) divergence (that moment minus 1). Only relevant when conditional = TRUE. The uncentered parameterization is the recommended one (Wang et al., 2026); the centered version corresponds to the parameterization of Huang and Pimentel (2025), retained so those results can be reproduced. Estimates of nu2.s are bounded below by 1 in the uncentered version and by 0 in the centered version.
verbose: if TRUE (default) prints steps of the fitting procedure.
warnings: should caret's warnings be printed? Default is FALSE. Note caret has many inconsistent and unnecessary warnings.
dml.fit: an object of class dml.
...: arguments passed to other methods.

Value

An object of class dml with the results of the DML procedure. The object is a list containing:

data: A list with the data used.
call: The original call used to fit the model.
info: A list with general information and arguments of the DML fitting procedure.
fits: A list with the predictions of each repetition.
results: A list with the results (influence functions and estimates) for each repetition.
coefs: A list with the estimates and standard errors for each repetition.

Details

Estimates of nu2.s are constrained in each cross-fitting repetition. For unconditional binary-treatment models the lower bound is \(1/[\hat p(1-\hat p)]\), where \(\hat p\) is the sample treatment share. For partially linear models it is the reciprocal of the empirical treatment variance, computed with divisor \(n\). Group-specific estimates use the corresponding within-group share or variance. The original influence functions are retained; sensitivity scales use the constrained estimates. In the centered conditional model, a zero sensitivity scale gives zero point-estimate bias. When confounding is positive, bound standard errors requiring the square-root delta method are unavailable at that scale.

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. doi:10.1162/REST.a.1705

Wang, J., Sant'Anna, P. H. C., Chernozhukov, V., and Cinelli, C. (2026). "Omitted Variable Bias in Difference-in-Differences Designs." Working Paper. (Basis of the conditional ATT/ATU parameterization.)

Examples

# loads package
library(dml.sensemakr)

## loads data
data("pension")

# set the outcome
y <- pension$net_tfa  # net total financial assets

# set the treatment
d <- pension$e401    # 401K eligibility

# set the covariates (a matrix)
x <- model.matrix(~ -1 + age + inc  + educ+ fsize + marr + twoearn + pira + hown, data = pension)

## compute income quartiles for group ATE.
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)
# run DML (nonparametric model)
## 2 folds (change as needed)
## 1 repetition (change as needed)
dml.401k <- dml(y, d, x, model = "npm", groups = g1, cf.folds = 2, cf.reps = 1)
#> Debiased Machine Learning
#> 
#>  Model: Nonparametric 
#>  Target: ate (unconditional) 
#>  Cross-Fitting: 2 folds, 1 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 2-fold cross-fitting 1 times
#> ======================================
#> 
#> -- Rep 1 -- Folds: 1  2  
#> 

# summary of results with median method (default)
summary(dml.401k, combine.method = "median")
#> 
#> Debiased Machine Learning
#> 
#>  Model: Nonparametric 
#>  Cross-Fitting: 2 folds, 1 reps 
#>  ML Method: outcome (yreg0:ranger, yreg1:ranger, R2 = 26.33%), treatment (ranger, R2 = 10.805%)
#>  Tuning: dirty 
#> 
#> Average Treatment Effect : 
#> 
#>     Estimate Std. Error t value  P(>|t|)    
#> ate     8354       1171   7.133 9.82e-13 ***
#> ---
#> 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     4700        800   5.875 4.24e-09 ***
#> g.ate.q2     2208       1371   1.610 0.107351    
#> g.ate.q3     6753       1865   3.621 0.000293 ***
#> g.ate.q4    19753       3984   4.958 7.13e-07 ***
#> ---
#> 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 2-fold cross-fitting. In order to reduce the variance that stems from sample splitting, we repeated the procedure 1 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.

# coef median method (default)
coef(dml.401k, combine.method = "median")
#>       ate  g.ate.q1  g.ate.q2  g.ate.q3  g.ate.q4 
#>  8353.976  4699.765  2208.199  6752.608 19752.677 

# se median method (default)
se(dml.401k, combine.method = "median")
#>       ate  g.ate.q1  g.ate.q2  g.ate.q3  g.ate.q4 
#> 1171.1649  800.0017 1371.3711 1864.8223 3984.2351 

# confint median method
confint(dml.401k, combine.method = "median")
#>               2.5 %    97.5 %
#> ate       6058.5352 10649.417
#> g.ate.q1  3131.7905  6267.739
#> g.ate.q2  -479.6385  4896.038
#> g.ate.q3  3097.6232 10407.592
#> g.ate.q4 11943.7193 27561.634

## Sensitivity Analysis

### Robustness Values
robustness_value(dml.401k, alpha = 0.05)
#>        ate   g.ate.q1   g.ate.q2   g.ate.q3   g.ate.q4 
#> 0.04824061 0.07683029 0.00000000 0.03536087 0.06785112 

### Confidence Bounds
confidence_bounds(dml.401k, cf.y = 0.03, cf.d = 0.04, level = 0.95)
#>                  lwr         upr
#> ate       1862.87448 14944.23750
#> g.ate.q1  1905.72981  7548.05070
#> g.ate.q2 -2629.90349  6343.04541
#> g.ate.q3    68.59443 13233.29542
#> g.ate.q4  6632.14260 33210.57632
#> 
#> Confidence level: point = 95%; region = 90%.
#> Sensitivity parameters: cf.y = 0.03; cf.d = 0.04; rho2 = 1.

### Contour Plots
ovb_contour_plot(dml.401k, cf.y = 0.03, cf.d = 0.04,
                bound.label = "Max Match (3x years)")