Implement DRScorer

[kgao]

Create initial implementation for drscorer for dr-learner based on dr-loss.

[kbattocchi]

This is a great start, but please address my comments, as well as the linting and testing failures in the automated checks.

[kbattocchi]

This doesn't look right to me - this should use the equation from the "Doubly Robust" section of https://github.com/py-why/EconML/blob/main/doc/spec/estimation/dr.rst

[kgao]

Mmm, if you mean the right equation should be the one right below "The Doubly Robust approach": Y_{i, t}^{DR} = g_t(X_i, W_i) + \frac{Y_i -g_t(X_i, W_i)}{p_t(X_i, W_i)} \cdot 1{T_i=t}
It's:
Y_DR[i,t] <- g_t(X[i], W[i]) + (Y[i] - g_t(X[i], W[i])) / p_t(X[i], W[i]) * (T[i] == t)
What I put here should be a short format combine with line 16 and line 18 (Where I put the input, weights)?

[kbattocchi]

Yeah, it's a bit tricky to write out - the key is that you're not multiplying by T at the end, you're multiplying by the indicator function selecting the specific case of T, and likewise you're dividing by the probability of that specific treatment, which is a bit awkward to express in pseudo-notation.

I think something like

Ydr(g,p) = g(X,W,T) + (Y - g(X,W,T)) / p_T(X,W)

would make it more obvious that there's only one term being included, it's not really being multiplied by T in a meaningful way, and it expresses the more-than-binary treatment outcome correctly. This does require writing out the arguments to g and p, but I think that's okay since otherwise it's hard to be precise about what's being computed.

[kgao]

Ok updated acoordingly.

[kbattocchi]

This code should basically reflect the code in

EconML/econml/dr/_drlearner.py

Lines 131 to 149 in 8b7fe33

	Y_pred, propensities = nuisances
	self.d_y = Y_pred.shape[1:-1] # track whether there's a Y dimension (must be a singleton)
	self.d_t = Y_pred.shape[-1] - 1 # track # of treatment (exclude baseline treatment)
	if (X is not None) and (self._featurizer is not None):
	X = self._featurizer.fit_transform(X)
	if self._multitask_model_final:
	ys = Y_pred[..., 1:] - Y_pred[..., [0]] # subtract control results from each other arm
	if self.d_y: # need to squeeze out singleton so that we fit on 2D array
	ys = ys.squeeze(1)
	weighted_sample_var = np.tile((sample_var / propensities**2).reshape((-1, 1)),
	self.d_t) if sample_var is not None else None
	filtered_kwargs = filter_none_kwargs(sample_weight=sample_weight,
	freq_weight=freq_weight, sample_var=weighted_sample_var)
	self.model_cate = self._model_final.fit(X, ys, **filtered_kwargs)
	else:
	weighted_sample_var = sample_var / propensities**2 if sample_var is not None else None
	filtered_kwargs = filter_none_kwargs(sample_weight=sample_weight,
	freq_weight=freq_weight, sample_var=weighted_sample_var)

where we take Y_pred from the nuisances and we form the doubly-robust estimate by subtracting Y_pred[:,0] from the other Y_pred[:,t] values. Since in the model fitting code the propensities nuisance is used only for adjusting sample_var, which we don't support here, I think you can ignore all of that code. So, I think this should look something more like:

Suggested change

	g, p = self.drlearner_._cached_values.nuisances
	Y = self.drlearner_._cached_values.Y
	T = self.drlearner_._cached_values.T
	Ydr = g + (Y - g) / p * T
	Y_pred, _ = self.drlearner_._cached_values.nuisances
	Y_dr = Y_pred[..., 1:] - Y_pred[..., [0]]

[kgao]

Good suggestion, thanks!