# Doubly robust learning with binary treatment and outcome

I'm trying to use doubly robust learning to estimate heterogenous treatment effects. My treatments T and outcomes y are both binary. I'm following the example listed under "How do I select the hyperparameters of the first stage models or the final model?" in the econML documentation.

For binary outcomes, the recommendation seems to be to simply wrap classifiers to make them behave like regressors, see for instance this github comment. A working example is provided below. However, the model outputs treatment effects which exceed 1, which doesn't make any sense for binary outcomes:

from econml.dr import DRLearner
import numpy as np
from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import GridSearchCV

# Wrapper to make classifier behave like a regressor
class RegWrapper:
def __init__(self, clf):
self._clf = clf

def fit(self, X, y):
self._clf.fit(X, y)
return self

def predict(self, X):
return self._clf.predict_proba(X)[:, 1]

### Simulate data
N = 10000
np.random.seed(0)
X = np.random.uniform(size=(N, 3))  # data affecting outcome and treatment
y = [int(np.random.uniform() < 0.1) for _ in range(N)]  # binary outcomes
T = [int(np.random.uniform() < 0.3) for _ in range(N)]  # binary treatment variables

# Propensity model - predicts T from X (classification)
model_clf = LogisticRegression()

# 'Regression' model - predicts y from X
model_reg = RegWrapper(LogisticRegression())

# Final model - predicts (I think) y_t=1 - y_t=0 from X
model_fin = GridSearchCV(
estimator=RandomForestRegressor(n_jobs=1),
param_grid={
'max_depth': [3, None],
'n_estimators': (10, 50)
}, cv=3, n_jobs=1, scoring='neg_mean_squared_error'
)

# Create and fit doubly robust learner
est = DRLearner(
model_regression=model_reg,
model_propensity=model_clf,
model_final=model_fin,
cv=2)
est.fit(y, T, X=X, W=None)

# Estimate heterogenous treatment effect for all points
point = est.effect(X, T0=0, T1=1)

print(max(point))  # prints 1.456


Can anyone point me in the right direction to get correct treatment effect estimates?

• The theory has not been worked out for the nonlinear case, and the code reflects that. The discussion you linked says that as well. Can you clarify why this is a problem? The ATE is unlikely to be outside the bounds, and the individual responses can be recoded to the appropriate boundary if they fall outside. Feb 16 at 22:23