Double check correctness of causal impact calculation in scikit-learn flavour of `DifferenceInDifferences`

[drbenvincent]

At the moment we have a TODO comment which needs some attention

# This is the coefficient on the interaction term
# TODO: THIS IS NOT YET CORRECT ?????
self.causal_impact = (
    self.y_pred_treatment[1] - self.y_pred_counterfactual[0]
)[0]

From memory, this is not wrong, but achieved the right effect in a different manner to the pymc flavour of difference in differences.

Either way, we need to double check and resolve.

[roesta07]

Hi @drbenvincent I recently worked on a DID project using Pymc. I would be happy to contribute if possible. Could you confirm on this bug?
Did you mean to directly reference the interaction parameter(causal_estimate) from the sklearn model's estimate as we did it for pymc model?

[drbenvincent]

Hi @roesta07. Sorry for the very delayed reply on this! The notification slipped through the net.

This whole section feels like it could be tightened up to be honest:

CausalPy/causalpy/experiments/diff_in_diff.py

Lines 193 to 206 in fdce5b0

	if isinstance(self.model, PyMCModel):
	# This is the coefficient on the interaction term
	coeff_names = self.model.idata.posterior.coords["coeffs"].data
	for i, label in enumerate(coeff_names):
	if "post_treatment" in label and self.group_variable_name in label:
	self.causal_impact = self.model.idata.posterior["beta"].isel(
	{"coeffs": i}
	)
	elif isinstance(self.model, RegressorMixin):
	# This is the coefficient on the interaction term
	# TODO: CHECK FOR CORRECTNESS
	self.causal_impact = (
	self.y_pred_treatment[1] - self.y_pred_counterfactual[0]
	)

For the PyMC models...

It's clear what the code is doing, but there is no comment about why it's looking for a coefficient containing "post_treatment". There could be some explanation there. And the way it's implemented with a loop and if (with not break) is a bit lazy.

Actually, looking a bit deeper I'm not 100% sure it's robust. When we look at the result summary in the notebook there is the main effect and the interaction term which contains post_treatment. It is clear from looking at the figure that the right term (the interaction) is being taken as the causal impact, but it's not clear that the code will always robustly get that coefficient.

Model coefficients:
    Intercept                     1.1, 94% HDI [1, 1.1]
    **post_treatment**[T.True]        0.99, 94% HDI [0.92, 1.1]
    group                         0.16, 94% HDI [0.094, 0.23]
    group:**post_treatment**[T.True]  0.5, 94% HDI [0.4, 0.6]
    sigma                         0.082, 94% HDI [0.066, 0.1]

NOTE: Looking back at the code, it's clear that post_treatment is a required column in the input data and formula. But it would be nice to have a comment explaining that. On a separate note, it would probably make sense to either a) create an issue to explain that more in the docs, or b) allow the user to use whatever column name they want, but to define that in the construction of the DifferenceInDifferences object. But like I say, that's a separate or side issue.

For the frequentist models...

I think the current line calculating the causal impact should equal the coefficient on the interaction term, but it would just feel more secure if we double checked that. Maybe even work it out both ways and check. And the indexing (self.y_pred_treatment[1] - self.y_pred_counterfactual[0]) with no explanation makes me feel uneasy. When #494 is merged, these should be xarray objects (if they aren't already), so at that point we might be able to use named/labelled indexing for greater clarity and certainty.

---

Sorry again about the delay. If you are still willing and able to take a look, then that would be great. Otherwise I guess these notes will be useful for me or someone else to move forward with this issue :)

[roesta07]

No worries, @drbenvincent . I’ll be mostly free for the next few weeks, so I’ll check and get back to you.

[roesta07]

NOTE: Looking back at the code, it's clear that post_treatment is a required column in the input data and formula. But it would be nice to have a comment explaining that. On a separate note, it would probably make sense to either a) create an issue to explain that more in the docs, or b) allow the user to use whatever column name they want, but to define that in the construction of the DifferenceInDifferences object. But like I say, that's a separate or side issue.

Hi @drbenvincent ;looking at our api, I think we should go with b) allow the user to use whatever column name they want, but to define that in the construction of the DifferenceInDifferences object.
We can add a new parameter post_treatment_variable_name to the __init__ method and set its default value to 'post_treatment' so that the previous codes does not break.

something like this along with fixing inside logic

def __init__(
    self,
    data: pd.DataFrame,
    formula: str,
    time_variable_name: str,
    group_variable_name: str,
    post_treatment_variable_name: str = 'post_treatment',  <<<- add new parameter which defaults to post_treatment
    model=None,
    **kwargs,
):

let me know, i can open up a new issue and make necessary fixes

[roesta07]

@drbenvincent and once the above is sorted out then the original issue i.e to get the coefficient of sklearn model would not be that tricky.

        elif isinstance(self.model, RegressorMixin):
            # This is the coefficient on the interaction term
            # TODO: THIS IS NOT YET CORRECT ?????
            # self.causal_impact = (
            #     self.y_pred_treatment[1] - self.y_pred_counterfactual[0]
            # )[0]


            print("### Updated Code")
            ## store the coefficient into dictionary {intercet:value}
            coef_map = dict(zip(self.labels, self.model.get_coeffs()))
            ## create and find the interaction term based on the values user provided
            interaction_term = f"{self.group_variable_name}:post_treatment"  # change , once the above issue is sorted out 
            matched_key = next((k for k in coef_map if interaction_term in k), None)
            att = coef_map.get(matched_key)
            self.causal_impact = att

for some reason for the model.summary() to work i had to initiate the sklearn model with fit_intercept=False other wise it would output the value of intercept to 0

result2 = cp.DifferenceInDifferences(
    df_new,
    formula="y ~ 1 + group*post_treatment",
    time_variable_name="post_treatment",
    group_variable_name="group",
    model=sk.linear_model.LinearRegression(fit_intercept=False),
)
print(result2.causal_impact)
print(result2.summary())