The IVDML package implements an instrumental variable (IV) estimator for potentially heterogeneous treatment effects in the presence of endogeneity as presented in Scheidegger, Guo and Bühlmann (2025). The estimator is based on double/debiased machine learning (DML) (Chernozhukov et al., 2018) and uses efficient machine learning instruments (MLIV) and kernel smoothing.
We assume that we observe $N$ i.i.d. copies of the model $$Y_i = \beta(A_i)D_i + g(X_i) + \epsilon_i.$$ $Y_i$ is the response variable and $D_i$ is the treatment variable. $X_i$ are the observed covariates and $A_i$ is a univariate continuous covariate with respect to which we want to consider heterogeneity (usually $A_i$ is a component of $X_i$). $\epsilon_i$ is an error term that satisfies $\mathbb E[\epsilon_i|X_i] = 0$, but we have endogeneity, meaning that $\mathbb E[\epsilon_i|D_i, X_i]\neq 0$. To deal with the endogeneity, we assume that we have access to an instrumental variable $Z_i$ that satisfies $\mathbb E[\epsilon_i|Z_i, X_i] = 0$. To goal is to conduct inference on the heterogeneous treatment effect $\beta(a)$ for some specific value $a$ of $A_i$.
We quickly describe the main idea of our estimator. Assume for a moment that the treatment effect $\beta(\cdot) = \beta$ is constant (i.e. does not depend on $A_i$). Then, $\beta$ is identified via $$\beta = \frac{\mathbb E[(Y_i - \mathbb E[Y_i|X_i])(\mathbb E[D_i|Z_i, X_i] - \mathbb E[D_i|X_i])]}{\mathbb E[(D_i - \mathbb E[D_i|X_i])(\mathbb E[D_i|Z_i, X_i] - \mathbb E[D_i|X_i])]}.$$ Hence, an estimate $\hat\beta$ of $\beta$ can be obtained by estimating the so-called nuisance functions $f(Z_i, X_i) = \mathbb E[D_i|Z_i, X_i]$, $\phi(X_i) = \mathbb E[D_i|X_i]$ and $l(X_i) = \mathbb E[Y_i|X_i]$ using arbitrary user-chosen machine-learning and calculating a sample version of the identification given above. In practice, this is done using a cross-fitting scheme (Chernozhukov et al. 2018).
If we allow for the treatment effect $\beta(\cdot)$ to depend on the univariate continuous quantity $A_i$, and $A_i$ is a component of $X_i$, we can make the identification $$\beta(a) = \frac{\mathbb E[(Y_i - \mathbb E[Y_i|X_i])(\mathbb E[D_i|Z_i, X_i] - \mathbb E[D_i|X_i])|A_i = a]}{\mathbb E[(D_i - \mathbb E[D_i|X_i])(\mathbb E[D_i|Z_i, X_i] - \mathbb E[D_i|X_i])|A_i = a]}.$$ conditional on $A_i = a$. In the sample version, we then use a kernel function $K(\cdot)$ (for example the Gaussian kernel $K(t) = \exp(-t^2/2)/\sqrt{2 \pi}$ ) to estimate the conditional expectations given $A_i = a$.
For a detailed discussion of the method, we refer to Scheidegger, Guo and Bühlmann (2025). We now demonstrate, how the IVDML package is used in practice.
You can install the released CRAN version of IVDML with
install.packages("IVDML")
You can install the development version of IVDML from GitHub with
devtools::install_github("cyrillsch/IVDML")
This is a basic example presenting the functionality of the IVDML package for the estimation and inference for homogeneous treatment effects. We first simulate a dataset with $N = 200$ observations.
set.seed(1)
N <- 200
Z <- rnorm(N)
X <- Z + rnorm(N)
H <- rnorm(N) # hidden confounding
D <- sin(2 * Z) - 2 * exp(-X^2) + H + 0.5 * rnorm(N)
Y <- -2 * D + tanh(X) - H + 0.5 * rnorm(N)
We see that the treatment effect $\beta = -2$ is homogeneous. Moreover, there is a hidden confounding variable $H$, which is why an instrument $Z$ is needed. An IVDML model can be fit using the following code.
library(IVDML)
fitted_ivdml <- fit_IVDML(Y = Y, D = D, Z = Z, X = X, ml_method = "gam", S_split = 10)
Here, ml_method = "gam"
indicates that the nuisance functions should
be estimated uing a generalized additive model. S_split = 10
indicates
that the cross-fitting is repeated 10 times with 10 random sample
splits. We can obtain a coefficient estimate and its associated standard
error using the coef()
and se()
methods.
print(coef(fitted_ivdml, iv_method = "mlIV"))
#> [1] -2.062829
print(se(fitted_ivdml, iv_method = "mlIV"))
#> [1] 0.1256957
Indicating iv_method = "mlIV"
gives the coefficient and standard error
estimate for the estimator based on the identification presented above.
Alternatively, we can also use an estimator that does not estimate the
nuisance function $f(Z_i, X_i) = \mathbb E[D_i|Z_i, X_i]$ and uses the
instruments linearly (i.e. using $Z_i - \mathbb E[Z_i|X_i]$ instead of
$\mathbb E[D_i|Z_i, X_i] - \mathbb E[D_i|X_i])]$). This estimator was
considered in Chernozhukov et al. (2018) and in Emmenegger and Bühlmann
(2021).
print(coef(fitted_ivdml, iv_method = "linearIV"))
#> [1] -2.192405
print(se(fitted_ivdml, iv_method = "linearIV"))
#> [1] 0.1991072
Observe that using iv_method = "mlIV"
leads to smaller estimated
standard error than iv_method = "linearIV"
. There are also methods for
confidence intervals for the estimated treatment effect, both standard
confidence intervals based on the point estimate and its standard error
and a robust confidence interval that is more robust with respect to
weak instrumental variables.
# mlIV
print(standard_confint(fitted_ivdml, iv_method = "mlIV", level = 0.95))
#> $CI
#> lower upper
#> -2.309188 -1.816470
#>
#> $level
#> [1] 0.95
#>
#> $heterogeneous_parameters
#> NULL
print(robust_confint(fitted_ivdml, iv_method = "mlIV", level = 0.95, CI_range = c(-10, 10)))
#> $CI
#> lower upper
#> -2.258353 -1.761513
#>
#> $level
#> [1] 0.95
#>
#> $message
#> [1] "The interval is contained in CI_range."
#>
#> $heterogeneous_parameters
#> NULL
# linearIV
print(standard_confint(fitted_ivdml, iv_method = "linearIV", level = 0.95))
#> $CI
#> lower upper
#> -2.582648 -1.802162
#>
#> $level
#> [1] 0.95
#>
#> $heterogeneous_parameters
#> NULL
print(robust_confint(fitted_ivdml, iv_method = "linearIV", level = 0.95, CI_range = c(-10, 10)))
#> $CI
#> lower upper
#> -2.533527 -1.596524
#>
#> $level
#> [1] 0.95
#>
#> $message
#> [1] "The interval is contained in CI_range."
#>
#> $heterogeneous_parameters
#> NULL
Here, CI_range = c(-10, 10)
specifies an a priori range for the robust
confidence interval.
We now present a basic example on how to estimate a heterogeneous treatment effect using IVDML. We first simulate a dataset with $N = 200$ observations.
set.seed(1)
N <- 200
Z <- rnorm(N)
X <- Z + rnorm(N)
A <- X
H <- rnorm(N) # hidden confounding
D <- sin(2 * Z) - 2 * exp(-X^2) + H + 0.5 * rnorm(N)
Y <- -2 * cos(A) * D + tanh(X) - H + 0.5 * rnorm(N)
We see that the treatment effect $\beta(A_i) = -2 \cos(A_i)$ is heterogeneous. We fit an IVDML model using the following code.
library(IVDML)
fitted_ivdml <- fit_IVDML(Y = Y, D = D, Z = Z, X = X, A = A, ml_method = "gam", S_split = 10)
The arguments of the fit_IVDML()
functions are the same as before with
the addition that we now also specify A = A
. This is however not
strictly necessary, as A
can also be provided when the coefficient,
standard error and confidence intervals are estimated, provided that A
is a deterministic function of X
(usually a component). The coef()
,
se()
, standard_confint()
and robust_confint()
methods now need
additional arguments a
, kernel_name
and bandwidth
. The normal
reference rule bandwidth according to Silverman (1986) can be obtained
using
h_normal <- bandwidth_normal(A = A)
print(h_normal)
#> [1] 0.4319743
For asymptotically valid inference, one needs to choose a slightly smaller bandwidth (undersmoothing). We use the following code for point estimate, standard errors and confidence intervals for the heterogeneous treatment effect $\beta(0)$ evaluated at $a = 0$.
print(coef(fitted_ivdml, iv_method = "mlIV", a = 0, A = A, kernel_name = "gaussian", bandwidth = h_normal/N^0.1))
#> [1] -2.003535
print(se(fitted_ivdml, iv_method = "mlIV", a = 0, A = A, kernel_name = "gaussian", bandwidth = h_normal/N^0.1))
#> [1] 0.218914
print(standard_confint(fitted_ivdml, iv_method = "mlIV", a = 0, A = A, kernel_name = "gaussian", bandwidth = h_normal/N^0.1), level = 0.95)
#> $CI
#> lower upper
#> -2.432598 -1.574471
#>
#> $level
#> [1] 0.95
#>
#> $heterogeneous_parameters
#> $heterogeneous_parameters$a
#> [1] 0
#>
#> $heterogeneous_parameters$kernel_name
#> [1] "gaussian"
#>
#> $heterogeneous_parameters$bandwidth
#> [1] 0.254305
print(robust_confint(fitted_ivdml, iv_method = "mlIV", a = 0, A = A, kernel_name = "gaussian", bandwidth = h_normal/N^0.1, level = 0.95, CI_range = c(-10, 10)))
#> $CI
#> lower upper
#> -2.362179 -1.318388
#>
#> $level
#> [1] 0.95
#>
#> $message
#> [1] "The interval is contained in CI_range."
#>
#> $heterogeneous_parameters
#> $heterogeneous_parameters$a
#> [1] 0
#>
#> $heterogeneous_parameters$kernel_name
#> [1] "gaussian"
#>
#> $heterogeneous_parameters$bandwidth
#> [1] 0.254305
As for the homogeneous treatment effect estimator, we can also calculate
theses quantities for iv_method = "linearIV"
instead of
iv_method = "mlIV"
.
print(coef(fitted_ivdml, iv_method = "linearIV", a = 0, A = A, kernel_name = "gaussian", bandwidth = h_normal/N^0.1))
#> [1] -1.765916
print(se(fitted_ivdml, iv_method = "linearIV", a = 0, A = A, kernel_name = "gaussian", bandwidth = h_normal/N^0.1))
#> [1] 0.4572597
print(standard_confint(fitted_ivdml, iv_method = "linearIV", a = 0, A = A, kernel_name = "gaussian", bandwidth = h_normal/N^0.1), level = 0.95)
#> $CI
#> lower upper
#> -2.6621283 -0.8697031
#>
#> $level
#> [1] 0.95
#>
#> $heterogeneous_parameters
#> $heterogeneous_parameters$a
#> [1] 0
#>
#> $heterogeneous_parameters$kernel_name
#> [1] "gaussian"
#>
#> $heterogeneous_parameters$bandwidth
#> [1] 0.254305
print(robust_confint(fitted_ivdml, iv_method = "linearIV", a = 0, A = A, kernel_name = "gaussian", bandwidth = h_normal/N^0.1, level = 0.95, CI_range = c(-10, 10)))
#> $CI
#> lower upper
#> -2.4715430 0.5675936
#>
#> $level
#> [1] 0.95
#>
#> $message
#> [1] "The interval is contained in CI_range."
#>
#> $heterogeneous_parameters
#> $heterogeneous_parameters$a
#> [1] 0
#>
#> $heterogeneous_parameters$kernel_name
#> [1] "gaussian"
#>
#> $heterogeneous_parameters$bandwidth
#> [1] 0.254305
More examples can be found in Scheiegger, Guo and Bühlmann (2025) and the associated GithHb repository IVDML_Application.
Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney Newey, and James Robins. Double/debiased machine learning for treatment and structural parameters. The Econometrics Journal, 21(1): C1–C68, 2018.
Corinne Emmenegger and Peter Bühlmann. Regularizing double machine learning in partially linear endogenous models. Electronic Journal of Statistics, 15(2):6461–6543, 2021.
Cyrill Scheidegger, Zijian Guo and Peter Bühlmann. Inference for heterogeneous treatment effects with efficient instruments and machine learning. Preprint, arXiv:2503.03530, 2025.
Bernard W. Silverman. Density estimation for statistics and data analysis. Chapman & Hall/CRC monographs on statistics and applied probability. Chapman and Hall, London, 1986.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.