We should implement get_params and set_params across all our cate estimators so that they can be used many times interchangeably with base estimators form sklearn. This will also help prettier copying and printing of the estimators.

When X is None, then const_marginal_effect() and effect() and marginal_effect() return a result that is of shape (1, d_y, d_t). The leading 1 seems weird to me. In this extreme case where the user wants average treatment effects, then we should not me adding the leading 1 in these results.

The specification indicates that controls should be supported by our 2SLS implementation (see here), but this isn't actually implemented.

As @vasilismsr pointed out in his review of #75, we are currently using percentile rather than pivot bootstrap intervals. We should implement pivot intervals as the default instead, and consider adding additional methods (e.g. bootstrap t-intervals as described in https://web.as.uky.edu/statistics/users/pbreheny/621/F12/notes/9-18.pdf)

Our fit and effect method signatures suggest that X is optional, but not all of our estimators work when X is not supplied (e.g. our double ML estimator).

Currently we only support an interface for getting intervals for target quantities. We might want to extend the api so as to return other quantities related to distributional properties of the estimated quantity, that might help in hypothesis testing. Examples include, p-values, t-statistics and standard errors.

In this line:
https://github.com/microsoft/EconML/blob/50baec376bd3b906b0f16720d3b97711acdcd98c/econml/dml.py#L170

Why don't we simply call effect(X, T0=0, T1=T_test_pred). This ought to return Y_test_res_pred directly.

SHAP is a popular new model interpretability method. Seems to be gaining a lot of traction. Are there plans to incorporate these causal estimates with model agnostic interpretations? Purely out of curiosity, with only having read a little of your work

Since ortho forest is asymptotically normal, we can perform inference using the asymptotic normal confidence intervals and variance estimation. This would be computationally much more viable than bootstrapping the whole estimation. Also variance can be estimated via the bootstrap of little bags approach proposed in the GRF paper: https://arxiv.org/abs/1610.01271

When using sample weights, the current bootstrap is no longer reasonable. Consider the case of one binary feature and were data are summarized and consist of means conditional on each feature value, with sample_weight corresponding to the number of original observations used to calculate the mean. Then our summarized data set consists of two samples each with weights. The bootstrap will be non-sense here as it will be drawing sometimes one sample twice or both samples. That’s not the intention. We somehow need to simulate what bootstrap would have done in the original sample.

This is not possible simply with sample weights. Ideally we would need access to all original samples. One idea: also solicit sample variances (i.e. variances within segment); which will be stored in the sample_vars fit argument. Then what bootstrap can do is treat each segment as distributed from a normal with the given mean and variance. Then draw fresh random samples from that normal. Moreover the number of samples that we draw from each segment should be determined as follows: we draw n segment ids each with probability proportional to the original number of samples in that segment. Then for each of these segment samples we draw a y value drawn from the normal approximation of that segment.

For now, the package correctly throws a warning that “ the bootstrap is used with sample weights and we advise against as the current implementation of bootstrap resampling does not take into account sample weights.”

I'm struggling to save models using pickle.
The coding example I'm running

import cloudpickle
pickled = cp.dumps(self.deepIvEst)
cp.loads(pickled)

I got the error: NameError: name 'tf' is not defined

I get an error that finishes with the exception: NameError: name 'tf' is not defined. I added tf as a custom object then

File "/home/***/.local/share/virtualenvs/RayServices-dMeB-KyP/lib/python3.6/site-packages/keras/layers/__init__.py", line 168, in deserialize
    printable_module_name='layer')
  File "/home/***/.local/share/virtualenvs/RayServices-dMeB-KyP/lib/python3.6/site-packages/keras/utils/generic_utils.py", line 147, in deserialize_keras_object
    list(custom_objects.items())))
  File "/home/***/.local/share/virtualenvs/RayServices-dMeB-KyP/lib/python3.6/site-packages/keras/engine/network.py", line 1062, in from_config
    process_node(layer, node_data)
  File "/home/***/.local/share/virtualenvs/RayServices-dMeB-KyP/lib/python3.6/site-packages/keras/engine/network.py", line 1012, in process_node
    layer(unpack_singleton(input_tensors), **kwargs)
  File "/home/***/.local/share/virtualenvs/RayServices-dMeB-KyP/lib/python3.6/site-packages/keras/engine/base_layer.py", line 451, in __call__
    output = self.call(inputs, **kwargs)
  File "/home/***/.local/share/virtualenvs/RayServices-dMeB-KyP/lib/python3.6/site-packages/keras/layers/core.py", line 716, in call
    return self.function(inputs, **arguments)
  File "/home/***/.local/share/virtualenvs/RayServices-dMeB-KyP/lib/python3.6/site-packages/econml/deepiv.py", line 159, in <lambda>
    samp = L.Lambda(lambda pms: _zero_grad(sample(*pms), pms), output_shape=(d_t,))
NameError: name '_zero_grad' is not defined

You may be wondering why I decided to start with the error: NameError: name 'tf' is not defined.

Well, it's because I wanted to highlight what I think the nature of this problem is. It's a scope issue. Everything that's imported to make the model works isn't used when unpickling a model. The big question is how to make this work?

One thing I'm thinking of is creating a separate function directly in each estimator that would allow us to reconstruct an entire model by using the hyperparameters and computed gradients, then re-embed them into the rest of the setup.

Hello,

I am working through the ORF paper (arXiv:1806.03467v2 [cs.LG] 12 Jul 2018) and I have a small question that I couldn't figure out on my own.

What are the results from the simulation when Y and T are residualized on W only and not on x? Is there a reason these aren't included in the paper but are in the R code GRF_treatment_effects.R?

Here is the snippet from the paper:

(3) GRF-Res: We perform a naive combination of double ML and GRF by first residualizing the treatments and outcomes on the both the features x and controls W , then running GRF R package on the residualized treatments Tˆ , residualized outcomes Yˆ , and features x.

and from the R code:
### Comparison #3
# Residualize on W, fit on x
forest <- causal_forest(input, data$res_Y_W, data$res_T_W)
tau_hat <- predict(forest, test_data)
results3[te_col] <- tau_hat$predictions

Thanks!

The DoublyRobustLearner currently does not do crossfitting for nuisance estimation and estimates propensities and direct models on the whole dataset.

Problem:
I think there should be some ways to identify subgroups suffer heterogeneous treatment effects. e.g tree structure can help us to tell that. It seems that SHAP targets only on 1-order (or maybe 2-order interaction), but what if we want more information, say higer orders?

I am facing an error on estimating conditional average treatment using the LinearCateEstimator(). I am following the example given in the documentation: https://econml.azurewebsites.net/spec/api.html#api-of-conditional-average-treatment-effect-package

Code:
est = LinearCateEstimator()
est.fit(y, T, X, W,Z)
print( est.effect(T0_test , T1_test, X))

est.effect(T0_test , T1_test, X) results in the following error:

IndexError Traceback (most recent call last)
in
1 est = LinearCateEstimator()
2 est.fit(y, T, X, W)
----> 3 print( est.effect(T0_test , T1_test, X))

~/.local/lib/python3.6/site-packages/econml/cate_estimator.py in effect(self, X, T0, T1)
148 # should it be possible to promote them to 2D arrays if that's what we saw during training?
149 eff = self.const_marginal_effect(X)
--> 150 m = shape(eff)[0]
151 if ndim(T0) == 0:
152 T0 = np.repeat(T0, m)

IndexError: tuple index out of range

Each of the T0_test, T1_test, X have the shape: (6500, 1)

I am not able to understand what leads to this error.

Also, I have faced issues in predicting CATE using DMLCateEstimator

est.effect(T0_test , T1_test, X) results in the following error for DMLCateEstimator:

python-input-30-ee269e75ba5f> in
1 est1 = DMLCateEstimator(model_y=LassoCV(), model_t=LassoCV())
2 est1.fit(Y, T, X, W)
----> 3 pred1= est1.effect(X, T0_test, T1_test)
4 print(pred1)
5 # est2 = DMLCateEstimator(model_y=RandomForestRegressor(),model_t=RandomForestRegressor())

~/.local/lib/python3.6/site-packages/econml/dml.py in effect(self, X, T0, T1)
137 T0 = self._one_hot_encoder.transform(reshape(self._label_encoder.transform(T0), (-1, 1)))[:, 1:]
138 T1 = self._one_hot_encoder.transform(reshape(self._label_encoder.transform(T1), (-1, 1)))[:, 1:]
--> 139 return super().effect(X, T0, T1)
140
141 def score(self, Y, T, X=None, W=None):

~/.local/lib/python3.6/site-packages/econml/cate_estimator.py in effect(self, X, T0, T1)
159 if ndim(eff) == ndim(dT): # y is a vector, rather than a 2D array
160 einsum_str = einsum_str.replace('y', '')
--> 161 return np.einsum(einsum_str, eff, dT)
162
163 def marginal_effect(self, T, X=None):

<array_function internals> in einsum(*args, **kwargs)

~/.local/lib/python3.6/site-packages/numpy/core/einsumfunc.py in einsum(*operands, **kwargs)
1354 # If no optimization, run pure einsum
1355 if optimize_arg is False:
-> 1356 return c_einsum(*operands, **kwargs)
1357
1358 valid_einsum_kwargs = ['out', 'dtype', 'order', 'casting']

ValueError: einstein sum subscripts string contains too many subscripts for operand 0

Documentation is outdated for most of our estimators due to the several updates we have done on the code. We need to an overall update to our docs to sync with the current state.

All our meta learner classes and the DoublyRobustLearner can be extended in a straightforward manner to accept multiple discrete treatments and not only binary. We should enable such functionality.

Dear all,

I have a problem implementing the cross price elasticity example. In particular I have a dataset consisting of an unbalanced panel of quantities sold of products and prices (so quantities sold and prices observed in time for different time slots each).

e.g. id_prod q_prod price Year

      1                  90             0.5        2001
      2                 10             3            2000
      2                 3              0.8          2001
      3.               23             2.9         2000   
      ...

I would like to recover the cross price elasticities: now, if I just use my "reshaped long" column of quantities as Y and the reshaped long column of prices as T, what I obtain is a single coefficient and NOT the matrix alpha whose I,j element is the elasticity of I when the price of j changes.

What I have done to obtain such a matrix has been to reshape the dataset in wide form and transpose it in order to obtain n (n = number of products) columns, one for each product and 14 rows, one for each year:
e.g. q_prod.1 q_prod.2 q_prodpr.3. price_1. price_2 price_3
2000 NaN. 10 23 NaN 3 2.9
2001. 90 3 7 0.5 0.8 0.6
2002. NaN NaN 110 NaN NaN 1.8

The code is below: y and t are reshaped data of quantities and prices; y_prova and t_prova represent the variables (I have taken a subset of products for computational speed) quantities and prices that I used as Y and T of the algorithm

y = pd.DataFrame(quantities_unstacked.T)
t = pd.DataFrame(prices_unstacked.T)
y_prova = y.iloc[0:13, 0:42648]
y_prova= np.array(y_prova, dtype=pd.DataFrame)
y_prova[pd.isnull(y_prova)] = 0

t_prova = t.iloc[0:13, 0:42648]
t_prova= np.array(t_prova, dtype=pd.DataFrame)
t_prova[pd.isnull(t_prova)] = 0

est = DMLCateEstimator(model_y=MultiTaskElasticNet(alpha=0.1),
model_t=MultiTaskElasticNet(alpha=0.1))
est.fit(y_prova, t_prova)

#a_hat[i,j] contains the elasticity of the demand of product i on the price of product j
a_hat = est.const_marginal_effect()
a_hat

The question is: is the procedure correct if I would like to obtain the cross elasticities price matrix (also taking into account of time dimension)? What if I would like to ignore NaN?

Sorry for the length of the question,

Thank you,

Federico

When the final stage of a DMLCateEstimator is a low dimensional linear model, then we can perform easier inference via asymptotic normality.

When samples represent summary statistics of many observations, then we need to pass in weights and variance properties of each sample. This will affect model fitting and inference.

Dear all,

sorry for rising another issue.
I cannot see any documentation for computing the standard errors. Is there a way to do it?

I would like to obtain the significance level of the following:

est.const_marginal_effect(X_test_first_atc1.T)

Thank you all!

The following is a slightly modification to the DeepIV notebook, but it no longer works. I guess it is because of the fact that I increased the max value of the X and T, which caused the overflow somewhere:

from econml.deepiv import DeepIVEstimator
from econml.bootstrap import BootstrapEstimator
import keras
import numpy as np
import matplotlib.pyplot as plt

get_ipython().run_line_magic('matplotlib', 'inline')


# ## Synthetic data
# 
# To demonstrate the Deep IV approach, we'll construct a synthetic dataset obeying the requirements set out above.  In this case, we'll take `X`, `Z`, `T` to come from the following distribution: 

# In[ ]:


# Initialize exogenous variables; normal errors, uniformly distributed covariates and instruments
e = np.random.normal(size=(445,1))
x = np.random.uniform(low=0.0, high=1000, size=(445,1))
z = np.random.uniform(low=0.0, high=1000, size=(445,1))

# Initialize treatment variable
t = np.sqrt((x+2) * z) + e

# Show the marginal distribution of t
plt.hist(t)
plt.xlabel("t")
plt.show()

plt.scatter(z[x < 1], t[x < 1], label='low X')
plt.scatter(z[(x > 4.5) * (x < 5.5)], t[(x > 4.5) * (x < 5.5)], label='moderate X')
plt.scatter(z[x > 9], t[x > 9], label='high X')
plt.legend()
plt.xlabel("z")
plt.ylabel("t")
plt.show()


# Here, we'll imagine that `Z` and `X` are causally affecting `T`; as you can see in the plot above, low or high values of `Z` drive moderate values of `T` and moderate values of `Z` cause `T` to have a bi-modal distribution when `X` is high, but a unimodal distribution centered on 0 when `X` is low.  The instrument is positively correlated with the treatment and treatments tend to be bigger at high values of x. The instrument has higher power at higher values of x 

# In[ ]:


# Outcome equation 
y = t*t / 10 - x*t / 10 + e

# The endogeneity problem is clear, the latent error enters both treatment and outcome equally
plt.scatter(t,z, label ='raw data')
tticks = np.arange(-2,12)
yticks2 = tticks*tticks/10 - 0.2 * tticks
yticks5 = tticks*tticks/10 - 0.5 * tticks
yticks8 = tticks*tticks/10 - 0.8 * tticks
plt.plot(tticks,yticks2, 'r--', label = 'truth, x=2')
plt.plot(tticks,yticks5, 'g--', label = 'truth, x=5')
plt.plot(tticks,yticks8, 'y--', label = 'truth, x=8')
plt.xlabel("t")
plt.ylabel("y")
plt.legend()
plt.show()


# `Y` is a non-linear function of `T` and `X` with no direct dependence on `Z` plus additive noise (as required).  We want to estimate the effect of particular `T` and `X` values on `Y`.
# 
# The plot makes it clear that looking at the raw data is highly misleading as to the treatment effect. Moreover the treatment effects are both non-linear and heterogeneous in x, so this is a hard problem!
# 
# ## Defining the neural network models
# 
# Now we'll define simple treatment and response models using the Keras `Sequential` model built up of a series of layers.  Each model will have an `input_shape` of 2 (to match the sums of the dimensions of `X` plus `Z` in the treatment case and `T` plus `X` in the response case).

# In[ ]:


treatment_model = keras.Sequential([keras.layers.Dense(128, activation='relu', input_shape=(2,)),
                                    keras.layers.Dropout(0.17),
                                    keras.layers.Dense(64, activation='relu'),
                                    keras.layers.Dropout(0.17),
                                    keras.layers.Dense(32, activation='relu'),
                                    keras.layers.Dropout(0.17)])

response_model = keras.Sequential([keras.layers.Dense(128, activation='relu', input_shape=(2,)),
                                   keras.layers.Dropout(0.17),
                                   keras.layers.Dense(64, activation='relu'),
                                   keras.layers.Dropout(0.17),
                                   keras.layers.Dense(32, activation='relu'),
                                   keras.layers.Dropout(0.17),
                                   keras.layers.Dense(1)])


# Now we'll instantiate the `DeepIVEstimator` class using these models.  Defining the response model *outside* of the lambda passed into constructor is important, because (depending on the settings for the loss) it can be used multiple times in the second stage and we want the same weights to be used every time.

# In[ ]:


keras_fit_options = { "epochs": 30,
                      "validation_split": 0.1,
                      "callbacks": [keras.callbacks.EarlyStopping(patience=2, restore_best_weights=True)]}

deepIvEst = DeepIVEstimator(n_components = 10, # number of gaussians in our mixture density network
                            m = lambda z, x : treatment_model(keras.layers.concatenate([z,x])), # treatment model
                            h = lambda t, x : response_model(keras.layers.concatenate([t,x])),  # response model
                            n_samples = 1, # number of samples to use to estimate the response
                            use_upper_bound_loss = False, # whether to use an approximation to the true loss
                            n_gradient_samples = 1, # number of samples to use in second estimate of the response (to make loss estimate unbiased)
                            optimizer='adam', # Keras optimizer to use for training - see https://keras.io/optimizers/ 
                            first_stage_options = keras_fit_options, # options for training treatment model
                            second_stage_options = keras_fit_options) # options for training response model


# ## Fitting and predicting using the model
# Now we can fit our model to the data:

# In[ ]:


boot_est = BootstrapEstimator(deepIvEst, n_bootstrap_samples=2, n_jobs=1)
boot_est.fit(Y=y,T=t,X=x,Z=z)


# And now we can create a new set of data and see whether our predicted effect matches the true effect `T*T-X*X`:

# In[ ]:


n_test = 500
for i, x_enum in enumerate([2, 5, 8]):
    t = np.linspace(0,10,num = 100)
    y_true = t*t / 10 - x_enum*t/10
    y_pred = boot_est.predict(t, np.full_like(t, x_enum))
    plt.plot(t, y_true, label='true y, x={0}'.format(x_enum),color='C'+str(i))
    plt.plot(t, y_pred, label='pred y, x={0}'.format(x_enum),color='C'+str(i),ls='--')
plt.xlabel('t')
plt.ylabel('y')
plt.legend()
plt.show()


# You can see that despite the fact that the response surface varies with x, our model was able to fit the data reasonably well. Where is does worst is where the instrument has the least power, which is in the low x case.  There it fits a straight line rather than a quadratic, which suggests that the regularization at least is perfoming well.  

# ## Estimating the Confidence Interval of the Model

# In[ ]:


treatment_effect_interval = boot_est.marginal_effect_interval(X=x, T=t, lower=1, upper=99)

Could you help and investigate why?

Before our next release we should update our documentation to make sure that it hasn't gotten out of sync with our implementation.

If we enable sample weights in DML then we need to have a weighted lasso implementation and a weighted lasso CV implementation, where both the fitting and the cross validation splitting take into account the sample weights.

Right now, some estimators require that Y and T are 2D arrays, even when they just have a single column. The library would be more ergonomic if we allowed them to be 1D vectors as well.

I tried adding the bootstrap wrapper to the file Deep IV Examples.ipynb with n_jobs more than 1, but it gives me some errors:

deepIvEst = DeepIVEstimator(n_components = 10, # number of gaussians in our mixture density network
                            m = lambda z, x : treatment_model(keras.layers.concatenate([z,x])), # treatment model
                            h = lambda t, x : response_model(keras.layers.concatenate([t,x])),  # response model
                            n_samples = 1, # number of samples to use to estimate the response
                            use_upper_bound_loss = False, # whether to use an approximation to the true loss
                            n_gradient_samples = 1, # number of samples to use in second estimate of the response (to make loss estimate unbiased)
                            optimizer='adam', # Keras optimizer to use for training - see https://keras.io/optimizers/ 
                            s1=100, # number of epochs to train treatment model
                            s2=100) # number of epochs to train response model
boot_est = BootstrapEstimator(deepIvEst, n_bootstrap_samples=100, n_jobs=4)
boot_est.fit(Y=y,T=t,X=x,Z=z)

[Parallel(n_jobs=4)]: Using backend ThreadingBackend with 4 concurrent workers.
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-14-2528bf95d3ca> in <module>
----> 1 boot_est.fit(Y=y,T=t,X=x,Z=z)

~/Dropbox (Cambridge University)/R301/EconML/econml/bootstrap.py in fit(self, *args, **named_args)
     56         Parallel(n_jobs=self._n_jobs, prefer='threads', verbose=3)(
     57             (obj.fit, [arg[inds] for arg in args], {arg: named_args[arg][inds] for arg in named_args})
---> 58             for obj, inds in zip(self._instances, indices)
     59         )
     60         return self

~/anaconda3/lib/python3.6/site-packages/joblib-0.13.2-py3.6.egg/joblib/parallel.py in __call__(self, iterable)
    932 
    933             with self._backend.retrieval_context():
--> 934                 self.retrieve()
    935             # Make sure that we get a last message telling us we are done
    936             elapsed_time = time.time() - self._start_time

~/anaconda3/lib/python3.6/site-packages/joblib-0.13.2-py3.6.egg/joblib/parallel.py in retrieve(self)
    831             try:
    832                 if getattr(self._backend, 'supports_timeout', False):
--> 833                     self._output.extend(job.get(timeout=self.timeout))
    834                 else:
    835                     self._output.extend(job.get())

~/anaconda3/lib/python3.6/multiprocessing/pool.py in get(self, timeout)
    642             return self._value
    643         else:
--> 644             raise self._value
    645 
    646     def _set(self, i, obj):

~/anaconda3/lib/python3.6/multiprocessing/pool.py in worker(inqueue, outqueue, initializer, initargs, maxtasks, wrap_exception)
    117         job, i, func, args, kwds = task
    118         try:
--> 119             result = (True, func(*args, **kwds))
    120         except Exception as e:
    121             if wrap_exception and func is not _helper_reraises_exception:

~/anaconda3/lib/python3.6/site-packages/joblib-0.13.2-py3.6.egg/joblib/_parallel_backends.py in __call__(self, *args, **kwargs)
    565     def __call__(self, *args, **kwargs):
    566         try:
--> 567             return self.func(*args, **kwargs)
    568         except KeyboardInterrupt:
    569             # We capture the KeyboardInterrupt and reraise it as

~/anaconda3/lib/python3.6/site-packages/joblib-0.13.2-py3.6.egg/joblib/parallel.py in __call__(self)
    223         with parallel_backend(self._backend, n_jobs=self._n_jobs):
    224             return [func(*args, **kwargs)
--> 225                     for func, args, kwargs in self.items]
    226 
    227     def __len__(self):

~/anaconda3/lib/python3.6/site-packages/joblib-0.13.2-py3.6.egg/joblib/parallel.py in <listcomp>(.0)
    223         with parallel_backend(self._backend, n_jobs=self._n_jobs):
    224             return [func(*args, **kwargs)
--> 225                     for func, args, kwargs in self.items]
    226 
    227     def __len__(self):

~/Dropbox (Cambridge University)/R301/EconML/econml/deepiv.py in fit(self, Y, T, X, Z)
    313         n_components = self._n_components
    314 
--> 315         treatment_network = self._m(z_in, x_in)
    316 
    317         # the dimensionality of the output of the network

<ipython-input-10-df8c54d69989> in <lambda>(z, x)
      1 deepIvEst = DeepIVEstimator(n_components = 10, # number of gaussians in our mixture density network
----> 2                             m = lambda z, x : treatment_model(keras.layers.concatenate([z,x])), # treatment model
      3                             h = lambda t, x : response_model(keras.layers.concatenate([t,x])),  # response model
      4                             n_samples = 1, # number of samples to use to estimate the response
      5                             use_upper_bound_loss = False, # whether to use an approximation to the true loss

~/anaconda3/lib/python3.6/site-packages/keras/engine/base_layer.py in __call__(self, inputs, **kwargs)
    455             # Actually call the layer,
    456             # collecting output(s), mask(s), and shape(s).
--> 457             output = self.call(inputs, **kwargs)
    458             output_mask = self.compute_mask(inputs, previous_mask)
    459 

~/anaconda3/lib/python3.6/site-packages/keras/engine/network.py in call(self, inputs, mask)
    562             return self._output_tensor_cache[cache_key]
    563         else:
--> 564             output_tensors, _, _ = self.run_internal_graph(inputs, masks)
    565             return output_tensors
    566 

~/anaconda3/lib/python3.6/site-packages/keras/engine/network.py in run_internal_graph(self, inputs, masks)
    719                                     kwargs['mask'] = computed_mask
    720                             output_tensors = to_list(
--> 721                                 layer.call(computed_tensor, **kwargs))
    722                             output_masks = layer.compute_mask(computed_tensor,
    723                                                               computed_mask)

~/anaconda3/lib/python3.6/site-packages/keras/layers/core.py in call(self, inputs)
    877 
    878     def call(self, inputs):
--> 879         output = K.dot(inputs, self.kernel)
    880         if self.use_bias:
    881             output = K.bias_add(output, self.bias, data_format='channels_last')

~/anaconda3/lib/python3.6/site-packages/keras/backend/tensorflow_backend.py in dot(x, y)
   1083         out = tf.sparse_tensor_dense_matmul(x, y)
   1084     else:
-> 1085         out = tf.matmul(x, y)
   1086     return out
   1087 

~/anaconda3/lib/python3.6/site-packages/tensorflow/python/ops/math_ops.py in matmul(a, b, transpose_a, transpose_b, adjoint_a, adjoint_b, a_is_sparse, b_is_sparse, name)
   2385       are both set to True.
   2386   """
-> 2387   with ops.name_scope(name, "MatMul", [a, b]) as name:
   2388     if transpose_a and adjoint_a:
   2389       raise ValueError("Only one of transpose_a and adjoint_a can be True.")

~/anaconda3/lib/python3.6/site-packages/tensorflow/python/framework/ops.py in __enter__(self)
   6081       if self._values is None:
   6082         self._values = []
-> 6083       g = _get_graph_from_inputs(self._values)
   6084       self._g_manager = g.as_default()
   6085       self._g_manager.__enter__()

~/anaconda3/lib/python3.6/site-packages/tensorflow/python/framework/ops.py in _get_graph_from_inputs(op_input_list, graph)
   5711         graph = graph_element.graph
   5712       elif original_graph_element is not None:
-> 5713         _assert_same_graph(original_graph_element, graph_element)
   5714       elif graph_element.graph is not graph:
   5715         raise ValueError("%s is not from the passed-in graph." % graph_element)

~/anaconda3/lib/python3.6/site-packages/tensorflow/python/framework/ops.py in _assert_same_graph(original_item, item)
   5647   if original_item.graph is not item.graph:
   5648     raise ValueError("%s must be from the same graph as %s." % (item,
-> 5649                                                                 original_item))
   5650 
   5651 

ValueError: Tensor("dense_1/kernel:0", shape=(2, 128), dtype=float32_ref) must be from the same graph as Tensor("concatenate_1/concat:0", shape=(?, 2), dtype=float32).

However, the problem is solved when I change n_jobs=1:

boot_est = BootstrapEstimator(deepIvEst, n_bootstrap_samples=10, n_jobs=1)

Is this a currency issue (i.e. DeepIV can only be trained on one job)? But it shouldn't be because in Bootstrapping, you should be able to do parallel processing ...

We need to add a version of DoublyRobustLearner that has a linear final CATE model and supports inference via asymptotically normal statistics.

Will pytorch support be offered for this package?

I can think of a number of biomedical applications that could largely benefit from these approaches.

When using inference, it's somewhat ambiguous whether a call like

est.effect_interval(X, low, high)

uses low and high as the optional treatments or as the lower and upper bounds of the confidence interval (the first interpretation is correct).

This confusion could be avoided by making T0 and T1 into keyword-only arguments, so that the user can call

est.effect_interval(X)
est.effect_interval(X, T0=low, T1=high)
est.effect_interval(X, lower=low, upper=high)

but not

est.effect_interval(X, low, high)

This would be a breaking change, but the safety benefits probably make it worthwhile anyway.

Right now when the neural networks are trained we're outputting a lot of info for each training epoch. We should find a way to give users more control (e.g. specifying extra keyword arguments for the training and response models' fit methods, which would let us set verbosity as well as other infrequently used arguments).

For orthogonal estimators that involve first stage nuisance estimation, we might want to add a lighter version of bootstrap inference that bootstraps only the computationally lighter second stage, while at the same time also being asymptotically valid.

Several hyperparameters in ortho forest should have a default setting that is data dependent and potentially set via some sophisticated method. We should enable such 'auto' defaults. Most important hyperparam seems to be the subsampling rate, which should definitely be set as a function of number of samples and features at fit time.

I'm running the following code to test different input shapes for DeepIV Estimator. I've been able to change input shapes via the input dimensions so far. The problem is that I've only been able to have a single dimension for my covariants, instruments, and treatments. I'm certain that's not the only way. I'm pretty sure the answer could be pretty simple, though I'm kind of jumping onto the deep-end with this library. I don't quite understand how.

To clarify the question is: How do I match different input shapes?

EconML Version: 0.4
Operating System: Ubuntu 18.10

Code that works:

from econml.deepiv import DeepIVEstimator
import keras
import numpy as np
import matplotlib.pyplot as plt
from count_service.primatives.base_models import treatment_model, response_model, create_callback

treatment_model = keras.Sequential([keras.layers.Dense(128, activation='relu', input_dim=8),
                                    keras.layers.Dropout(0.17),
                                    keras.layers.Dense(64, activation='relu'),
                                    keras.layers.Dropout(0.17),
                                    keras.layers.Dense(32, activation='relu'),
                                    keras.layers.Dropout(0.17)])

response_model = keras.Sequential([keras.layers.Dense(128, activation='relu', input_dim=8),
                                keras.layers.Dropout(0.17),
                                keras.layers.Dense(64, activation='relu'),
                                keras.layers.Dropout(0.17),
                                keras.layers.Dense(32, activation='relu'),
                                keras.layers.Dropout(0.17),
                                keras.layers.Dense(1)])


treatment_model2 = keras.Sequential([keras.layers.Dense(128, activation='relu', input_dim=2),
                                    keras.layers.Dropout(0.17),
                                    keras.layers.Dense(64, activation='relu'),
                                    keras.layers.Dropout(0.17),
                                    keras.layers.Dense(32, activation='relu'),
                                    keras.layers.Dropout(0.17)])

response_model2 = keras.Sequential([keras.layers.Dense(128, activation='relu', input_dim=2),
                                keras.layers.Dropout(0.17),
                                keras.layers.Dense(64, activation='relu'),
                                keras.layers.Dropout(0.17),
                                keras.layers.Dense(32, activation='relu'),
                                keras.layers.Dropout(0.17),
                                keras.layers.Dense(1)])

n = 5000

# Initialize exogenous variables; normal errors, uniformly distributed covariates and instruments
e = np.random.normal(size=(n,1))
x = np.random.uniform(low=0.0, high=10.0, size=(n,1))
z = np.random.uniform(low=0.0, high=10.0, size=(n,1))

e_single = np.random.normal(size=(n,))
x_single = np.random.uniform(low=0.0, high=10.0, size=(n,))
z_single = np.random.uniform(low=0.0, high=10.0, size=(n,))

e_large = np.random.normal(size=(n,3))
x_large = np.random.uniform(low=0.0, high=10.0, size=(n,3))
z_large = np.random.uniform(low=0.0, high=10.0, size=(n,3))

t_single = np.sqrt((x_single+2) * z_single) + e_single

y_single = t_single*t_single / 10 - x_single*t_single / 10 + e_single


# Initialize treatment variable
t = np.sqrt((x+2) * z) + e

# Outcome equation 
y = t*t / 10 - x*t / 10 + e


t_large = np.sqrt((x_large+2) * z_large) + e_large
y_large = t_large*t_large / 10 - x_large*t_large / 10 + e_large

deepIvEst = DeepIVEstimator(n_components = 10, # number of gaussians in our mixture density network
                            m = lambda z, x : treatment_model2(keras.layers.concatenate([z,x])), # treatment model
                            h = lambda t, x : response_model2(keras.layers.concatenate([t,x])),  # response model
                            n_samples = 1, # number of samples to use to estimate the response
                            use_upper_bound_loss = False, # whether to use an approximation to the true loss
                            n_gradient_samples = 1, # number of samples to use in second estimate of the response (to make loss estimate unbiased)
                            optimizer='adam', # Keras optimizer to use for training - see https://keras.io/optimizers/ 
                            first_stage_options = create_callback(), # options for training treatment model
                            second_stage_options = create_callback()) # options for training response model




deepIvEst.fit(Y=y_single,T=t_single,X=x_single,Z=z_single)
for i, x in enumerate([2, 5, 8]):
    t = np.linspace((1, 4),(1, 20),num = 1000)
    y_true = t*t / 10 - x*t/10
    fuller = np.full_like(t, x)
    y_pred = deepIvEst.predict(fuller, np.full_like(fuller, x))

Code that doesn't work:

n = 5000

# Initialize exogenous variables; normal errors, uniformly distributed covariates and instruments
e = np.random.normal(size=(n,1))
x = np.random.uniform(low=0.0, high=10.0, size=(n,1))
z = np.random.uniform(low=0.0, high=10.0, size=(n,1))

e_single = np.random.normal(size=(n,))
x_single = np.random.uniform(low=0.0, high=10.0, size=(n,))
z_single = np.random.uniform(low=0.0, high=10.0, size=(n,))

e_large = np.random.normal(size=(n,3))
x_large = np.random.uniform(low=0.0, high=10.0, size=(n,3))
z_large = np.random.uniform(low=0.0, high=10.0, size=(n,3))

t_single = np.sqrt((x_single+2) * z_single) + e_single

y_single = t_single*t_single / 10 - x_single*t_single / 10 + e_single


# Initialize treatment variable
t = np.sqrt((x+2) * z) + e

# Outcome equation 
y = t*t / 10 - x*t / 10 + e


t_large = np.sqrt((x_large+2) * z_large) + e_large
y_large = t_large*t_large / 10 - x_large*t_large / 10 + e_large

deepIvEst = DeepIVEstimator(n_components = 10, # number of gaussians in our mixture density network
                            m = lambda z, x : treatment_model2(keras.layers.concatenate([z,x])), # treatment model
                            h = lambda t, x : response_model2(keras.layers.concatenate([t,x])),  # response model
                            n_samples = 1, # number of samples to use to estimate the response
                            use_upper_bound_loss = False, # whether to use an approximation to the true loss
                            n_gradient_samples = 1, # number of samples to use in second estimate of the response (to make loss estimate unbiased)
                            optimizer='adam', # Keras optimizer to use for training - see https://keras.io/optimizers/ 
                            first_stage_options = create_callback(), # options for training treatment model
                            second_stage_options = create_callback()) # options for training response model




deepIvEst.fit(Y=y_single,T=t_single,X=x_large,Z=z_large)
for i, x in enumerate([2, 5, 8]):
    t = np.linspace((1, 4),(1, 20),num = 1000)
    y_true = t*t / 10 - x*t/10
    fuller = np.full_like(t, x)
    y_pred = deepIvEst.predict(fuller, np.full_like(fuller, x))

I get this error (rightfully so)

Traceback (most recent call last):
  File "count_service/junk/direct_sample.py", line 81, in <module>
    deepIvEst.fit(Y=y_single,T=t_single,X=x_large,Z=z_large)
  File "/home/kevin/.local/share/virtualenvs/RayServices-dMeB-KyP/lib/python3.6/site-packages/econml/deepiv.py", line 323, in fit
    treatment_network = self._m(z_in, x_in)
  File "count_service/junk/direct_sample.py", line 69, in <lambda>
    m = lambda z, x : treatment_model2(keras.layers.concatenate([z,x])), # treatment model
  File "/home/kevin/.local/share/virtualenvs/RayServices-dMeB-KyP/lib/python3.6/site-packages/keras/engine/base_layer.py", line 451, in __call__
    output = self.call(inputs, **kwargs)
  File "/home/kevin/.local/share/virtualenvs/RayServices-dMeB-KyP/lib/python3.6/site-packages/keras/engine/network.py", line 570, in call
    output_tensors, _, _ = self.run_internal_graph(inputs, masks)
  File "/home/kevin/.local/share/virtualenvs/RayServices-dMeB-KyP/lib/python3.6/site-packages/keras/engine/network.py", line 727, in run_internal_graph
    layer.call(computed_tensor, **kwargs))
  File "/home/kevin/.local/share/virtualenvs/RayServices-dMeB-KyP/lib/python3.6/site-packages/keras/layers/core.py", line 908, in call
    output = K.dot(inputs, self.kernel)
  File "/home/kevin/.local/share/virtualenvs/RayServices-dMeB-KyP/lib/python3.6/site-packages/keras/backend/tensorflow_backend.py", line 1133, in dot
    out = tf.matmul(x, y)
  File "/home/kevin/.local/share/virtualenvs/RayServices-dMeB-KyP/lib/python3.6/site-packages/tensorflow/python/util/dispatch.py", line 180, in wrapper
    return target(*args, **kwargs)
  File "/home/kevin/.local/share/virtualenvs/RayServices-dMeB-KyP/lib/python3.6/site-packages/tensorflow/python/ops/math_ops.py", line 2647, in matmul
    a, b, transpose_a=transpose_a, transpose_b=transpose_b, name=name)
  File "/home/kevin/.local/share/virtualenvs/RayServices-dMeB-KyP/lib/python3.6/site-packages/tensorflow/python/ops/gen_math_ops.py", line 5925, in mat_mul
    name=name)
  File "/home/kevin/.local/share/virtualenvs/RayServices-dMeB-KyP/lib/python3.6/site-packages/tensorflow/python/framework/op_def_library.py", line 788, in _apply_op_helper
    op_def=op_def)
  File "/home/kevin/.local/share/virtualenvs/RayServices-dMeB-KyP/lib/python3.6/site-packages/tensorflow/python/util/deprecation.py", line 507, in new_func
    return func(*args, **kwargs)
  File "/home/kevin/.local/share/virtualenvs/RayServices-dMeB-KyP/lib/python3.6/site-packages/tensorflow/python/framework/ops.py", line 3616, in create_op
    op_def=op_def)
  File "/home/kevin/.local/share/virtualenvs/RayServices-dMeB-KyP/lib/python3.6/site-packages/tensorflow/python/framework/ops.py", line 2027, in __init__
    control_input_ops)
  File "/home/kevin/.local/share/virtualenvs/RayServices-dMeB-KyP/lib/python3.6/site-packages/tensorflow/python/framework/ops.py", line 1867, in _create_c_op
    raise ValueError(str(e))
ValueError: Dimensions must be equal, but are 6 and 2 for 'sequential_5/dense_15/MatMul' (op: 'MatMul') with input shapes: [?,6], [2,128].

The key difference here is what I'm training the program with.

I've been reading up on input shapes, and part of what I'm reading suggests that I might be approaching this problem wrong. I can get input shapes of equal size working. How should I approach this?

We should create a BinaryCateEstimator base class that all estimators that solely accept binary treatments can inherit from. Alternatively we can create a DiscreteTreatmentCateEstimator which could also accept multiple discrete treatments. We could then make all our meta learners and the DoublyRobustLearner inherit from these classes. There are several things that can simplify in our API if we know that the cate estimator will only accept binary/discrete treatments and hence this can simplify development of child estimators for such classes.

Several of our estimators have notebooks demonstrating proper usage, but we're missing a few.

The logic of the effect in the ortho forest is wrong as it returns the const_marginal_effects for all treatments as opposed to the effect between two treatments.

It might be that the weighted model wrapper has a small bug in the case of linear models that fit an intercept:

def _weighted_inputs(self, X, y, sample_weight):
    normalized_weights = X.shape[0] * sample_weight / np.sum(sample_weight)
    sqrt_weights = np.sqrt(normalized_weights)
    weight_mat = np.diag(sqrt_weights)
    return np.matmul(weight_mat, X), np.matmul(weight_mat, y)

If there is an intercept then we also want the intercept to be multiplied by the weight, so this transformation would be invalid. The user can fix it by augmenting X with an intercept and then calling a linear model without intercept. But this might be cumbersome.

We might want to add verbosity params to most of our cate estimators which then each estimator can choose in an adhoc manner to print progress messages. This might be good since some methods might take time to run and progress messages can be good.

The OrthoForest already supports discrete treatments, but double ML doesn't. We should provide an initializer option which does the following:

• One-hot-encodes the treatments
• Uses predict_proba on the T model instead of predict

Panel data have several intricacies that we might want to create a panel wrapper to our estimators that call our estimators in the right way when the data have a panel structure. We might want to add accompanying applications of the panel wrapper to pricing applications.

I was wondering if the deepiv pacakge has causal inference built in? For example, in the causal random forest package, a confidence interval can be estimated for the treatment effects. Is there similar functionality here? So far I am only able to see a predicted value for treatment effect, but no standard errors.

Dear all,

thanks again for this wonderful ML tool.
I would like to understand if it is possible, in my case to have t-statistics exploiting the bootstrap algorithm you provided. In particular I have employed both

est = DMLCateEstimator(model_y=MultiTaskElasticNet(alpha=0.1), model_t=MultiTaskElasticNet(alpha=0.1))and est = DMLCateEstimator(model_y=MultiTaskElasticNet(alpha=0.1), model_t=MultiTaskElasticNet(alpha=0.1),model_final=Lasso(alpha=0.1),featurizer=PolynomialFeatures(degree=1),random_state=123)
obtaining estimates for my Databases (15 in total). What I would like to do is computing SE or t-statistics directly by exploiting the bootstrap method you provided. My attempt was the following but I end up making a mess (basically I created a Y a T and an X suitable for the method and simply inputted in boot_est):

for i in range(1,16): df_part_q['df_part_q_'+str(i)] = pd.DataFrame() df_part_p['df_part_p_'+str(i)] = pd.DataFrame() Y['Y_'+ str(i)] = [] T['T_'+ str(i)] = [] X['X_'+ str(i)] = []

for j in range(1,16): df_part_q['df_part_q_'+str(j)] = pd.concat([quantities_unstacked['quantities_unstacked'+str(j)].T, quantities_unstacked_othertatc1['quantities_unstacked_otheratc1_'+str(j)].T], axis = 1)
df_part_p['df_part_p_'+str(j)] = pd.concat([prices_unstacked['prices_unstacked'+str(j)].T, prices_unstacked_othertatc1['prices_unstacked_otheratc1_'+str(j)].T], axis = 1) df_part_r['df_part_r_'+str(j)] = db['db_'+str(j)]['Recalls rest product']

for j in range(1,16): Y['Y_'+str(j)].append(df_part_q['df_part_q_'+str(j)].values) T['T_'+str(j)].append(df_part_p['df_part_p_'+str(j)].values) X['X_'+str(j)].append(df_part_r['df_part_r_'+str(j)].values)

from econml.bootstrap import BootstrapEstimator boot_est=BootstrapEstimator(DMLCateEstimator(model_y=MultiTaskElasticNetCV(cv=3),model_t=MultiTaskElasticNetCV(cv=3)) ,n_bootstrap_samples=20)

boot_estimate = {} for j in range(1,16): boot_estimate['boot_est'+str(j)] = [] for j in range(1,16): boot_estimate['boot_est'+str(j)].append(boot_est.fit(Y['Y_'+str(j)][0], T['T_'+str(j)][0].reshape(-1, 1)))

te_pred_interval = {} for j in range(1,16): te_pred_interval['te_pred_interval'+str(j)] = []

for j in range(1,16): te_pred_interval['te_pred_interval'+str(j)].append(boot_estimate['boot_est'+str(j)][0].const_marginal_effect())

Can you please help me in this? :)

Thank you,

Federico

When the final stage of a DMLCateEstimator is a high dimensional sparse linear model, then we can perform easier inference via the debiased lasso and using asymptotic normality. We need to implement a debiased lasso estimator that has inference capabilities (e.g. supports confidence intervals)

Problem

When we use the ordinary classification models provided by sklearn, we get poor results because metalearners call the predict methods that return 0 or 1.

It would be because metalearners in EconML now focus on the continuous outcome.

Solution

I modified the predict methods in two ways:

• create wrappers of sklearn models
• override methods after the initialization of metalearners

I uploaded a jupyter notebook in gist to explain the modification.

Discussion and Proposition

I would like to find better solutions to handle the binary outcome (from my point of view, it would not be reasonable to add an option to handle the binary outcome in EconML right now). Do you have any ideas?

And I would appreciate it if you could provide a supplementary explanation about handling the binary outcome anywhere in your document.

In many cases (single treatment, single response models with binary treatments or linear responses) setting T0 to 0 and T1 to 1 for all rows of the data would be a sensible default that would save users from having to explicitly provide arrays.

The following example throws an exception because Lasso's predict returns a vector even when the input was a one-column array (originally reported as part of #85):

import numpy as np
from econml.dml import DMLCateEstimator
from sklearn.linear_model import Lasso

Y, T, X = [np.random.normal(size=(10,d)) for d in [1,1,2]]
DMLCateEstimator(model_t=Lasso(), model_y=Lasso()).fit(Y,T,X)

We ought to be able to identify this and reshape if necessary, either in DML or higher up in the class hierarchy.

Our doubly-robust metalearner ought to be able to use controls W in addition to features X.