[question] multiple treatment variables and final model

[ieriii]

Let's say I have the following set up:

Y = outcome variable
T = list of treatment variables [T1, T2, T3]
X = confounders

The functional form is:

Y = theta_1 * T1 + theta_2 * T2 + theta_3 * T3 + g(X); and
T1 = m(x) + e
T2 = m(x) + e
T3 = m(x) + e

I would like to use double ML to do the following:

use ML to estimate y_hat = g_hat(X) + u
use ML to estimate T1_hat = m_hat(X) + v1
use ML to estimate T2_hat = m_hat(X) + v2
use ML to estimate T3_hat = m_hat(X) + v3
regress the residuals u against v1, v2 and v3 using linear regression.

The intuition of this approach is to (i) net off the confounders X from Y, T1, T2 and T3, (ii) take the residuals and (iii) estimate the impact of the treatment variables on Y.

Is this something that can be implemented with doubleML? If so, how?
I have looked at the docs and source code and I think I can do point 1-5 above (correct me if you think I'm wrong) but cannot identify whether point 6 could be possible.
Thank you.

[MalteKurz]

Some details on the multiple treatment case can be found in the user guide: https://docs.doubleml.org/stable/guide/sim_inf.html. This should then be applicable to your functional form. Note however, that it is not explicitly imposing the form Y = theta_1 * T1 + theta_2 * T2 + theta_3 * T3 + g(X) as you intended. The code to estimate joint confidence intervals for the effect of (T1, T2, T3) on Y as described in the user guide is given below:

import doubleml as dml
import numpy as np
from sklearn.base import clone
from sklearn.linear_model import LassoCV, LinearRegression

np.random.seed(1234)
n_obs = 1000
dim_x = 100
X = np.random.normal(size=(n_obs, dim_x))

theta = np.array([1., 1.5, 2.25])

beta = [1 / (k**2) for k in range(1, dim_x + 1)]
gamma = [1 / (k**2) for k in range(1, dim_x + 1)]

T1 = np.dot(X, gamma) + np.random.normal(size=(n_obs,))
T2 = np.dot(X, gamma) + np.random.normal(size=(n_obs,))
T3 = np.dot(X, gamma) + np.random.normal(size=(n_obs,))

T = np.vstack([T1, T2, T3]).T

y = np.dot(T, theta) + np.dot(X, beta) + np.random.standard_normal(size=(n_obs,))

dml_data = dml.DoubleMLData.from_arrays(X, y, T)
learner = LassoCV()
ml_g = clone(learner)
ml_m = clone(learner)

dml_plr = dml.DoubleMLPLR(dml_data, ml_g, ml_m)

print(dml_plr.fit().bootstrap().confint(joint=True))
print(dml_plr.p_adjust())
print(dml_plr.p_adjust(method='bonferroni'))

       2.5 %    97.5 %
d1  0.969056  1.110728
d2  1.383672  1.533448
d3  2.244911  2.381489
        coef  pval
d1  1.039892   0.0
d2  1.458560   0.0
d3  2.313200   0.0
        coef           pval
d1  1.039892  9.642752e-261
d2  1.458560   0.000000e+00
d3  2.313200   0.000000e+00

The approach you suggested, is not directly implemented in DoubleML. I haven't checked the theoretical details. However, there is the option to export predictions and with that I assume that you can do the final stage (in the style suggested by you) by hand with the following code:

dml_data = dml.DoubleMLData.from_arrays(X, y, T, use_other_treat_as_covariate=False)
learner = LassoCV()
ml_g = clone(learner)
ml_m = clone(learner)

dml_plr = dml.DoubleMLPLR(dml_data, ml_g, ml_m)
dml_plr.fit(store_predictions=True)

g_hat = dml_plr.predictions['ml_g'][:,0,0]
m_hat = dml_plr.predictions['ml_m'][:,0,:]

u_hat = y - g_hat
v_hat = T - m_hat

reg = LinearRegression(fit_intercept=False).fit(v_hat, u_hat)
reg.coef_

array([1.01816006, 1.43210268, 2.29182062])

[ieriii]

Thank you, MalteKurz. This is very helpful.

Note however, that it is not explicitly imposing the form Y = theta_1 * T1 + theta_2 * T2 + theta_3 * T3 + g(X) as you intended.

I understand this means that your first code snippet estimates each theta individually, i.e. using 3 separate OLS regressions (u_hat against v1_hat; u_hat against v2_hat; u_hat against v3_hat).

In your second code snippet, we can impose the linear functional. However, I'm not sure about the validity of each coefficient standard errors and confidence interval in the regression summary - any suggestions welcomed.

In summary, the theta estimated using the first code snippet might be different compared to the second snipped since in the first approach we are running a regression without 'controlling' for other residuals.

Let me know if any of the point above is not correct. Appreciate your help.

[MalteKurz]

Note however, that it is not explicitly imposing the form Y = theta_1 * T1 + theta_2 * T2 + theta_3 * T3 + g(X) as you intended.

I understand this means that your first code snippet estimates each theta individually, i.e. using 3 separate OLS regressions (u_hat against v1_hat; u_hat against v2_hat; u_hat against v3_hat).

Not really. It means that e.g. when estimating E(Y|X, T_2, T_3) for the partialling out score we just fit the specified ML model to learn the effect of (X, T_2, T_3) on Y. However, one may be able to improve the prediction by explicitly using that T_2 and T_3 enter linearly as stated above.

In your second code snippet, we can impose the linear functional. However, I'm not sure about the validity of each coefficient standard errors and confidence interval in the regression summary - any suggestions welcomed.

The second code snippet just shows how the procedure suggested by you, i.e., "regress the residuals u against v1, v2 and v3 using linear regression." can be technically implemented. I am pretty sure that the standard errors obtained should better not be interpreted. When estimating them, the first-stage is ignored which usually results in too small standard errors. I further don't really understand how the Neyman orthogonal moment conditions corresponding to this procedure would look like. Therefore, it is also difficult to give you an advice how to obtain valid standard errors for this procedure.

In summary, the theta estimated using the first code snippet might be different compared to the second snipped since in the first approach we are running a regression without 'controlling' for other residuals.

I don't really agree with this summary. As stated above the first code snippet just shows how to implement the procedure suggested by you. What you exactly mean with "controlling for other residuals" is not clear to me. The model behind the first code snippet is briefly summarized here: https://docs.doubleml.org/stable/guide/sim_inf.html. It is closely related to Belloni et al. (2018). For each treatment effect that should be estimated there is one identifying orthogonal moment condition. Stacking them together then also allows to obtain joint confidence intervals etc.

cc @MartinSpindler