Is `calculate_impact` correct?

[drbenvincent]

When I was looking at the multi-cell geolift notebook I started wondering if the calculate_impact method of the PyMCModel class was correct. In the post intervention period in the top plot you can see that the posterior expectation has narrow HDI's and the data points are far away. Yet when we look in the causal impact plot which simply looks at the difference between the data and the posterior mu distribution, we see larger HDI's which sometimes overlap with zero.

This was confirmed, the current implementation calculates the causal impact as the difference between the data and the posterior predictive distribution.

CausalPy/causalpy/pymc_models.py

Lines 169 to 173 in 714c48b

	def calculate_impact(
	self, y_true: xr.DataArray, y_pred: az.InferenceData
	) -> xr.DataArray:
	impact = y_true - y_pred["posterior_predictive"]["y_hat"]
	return impact.transpose(..., "obs_ind")

I think this should instead be a comparison between the data and the posterior expectation. So instead of y_pred["posterior_predictive"]["y_hat"] we should have y_pred["posterior_predictive"]["mu"]?

If so, then the implications are that our estimates of the causal impact will increase in precision. If the current code is in error, then we are not getting biased estimates, we are just getting less precise estimates than we should be getting out.

Making this change results in this...

[juanitorduz]

When I used to do this “by hand” I always considered the whole posterior distribution instead of the expectation. So, what if we just not simply change the top plot to reflect this?

I think in the CausalImpact paper by Google they also use the whole posterior (but I could not find any formula, maybe I missed it). That being said, I have no theoretical argument to look to the whole posterior (the expectation can be easier to explain, but that is not an argument).

[cluhmann]

Is the issue just the agreement between the plots? If so, what is being plotted in the top panel?

[drbenvincent]

I also chatted with @NathanielF about this.

Issue 1 - consistency between the plots

As @cluhmann also said, we don't want inconsistency between the plots. Currently the top plot HDI shows the model expectation and the lower plot shows the full posterior predictive distribution. This fix is easy enough.

We can always allow the user to specify which distribution they want in a new kwarg, the question is...

Issue 2 - is there a correct, or good default, distribution to pick

It could be that there's an academic paper out there on causal Bayesian inference that makes a good case for one of the other.

I think initially the causal impact was calculated using the posterior predictive, but there was a PR which changed this to be based on model expectation https://github.com/pymc-labs/CausalPy/pull/145/files#diff-4178e87a1cb5ec9e2aa470845a710302e3d87ed3ca9311df0a571b89e4e49bd7

but somehow that has drifted back to being based on the posterior predictive (yhat)

CausalPy/causalpy/pymc_models.py

Lines 169 to 173 in 8b26d42

	def calculate_impact(
	self, y_true: xr.DataArray, y_pred: az.InferenceData
	) -> xr.DataArray:
	impact = y_true - y_pred["posterior_predictive"]["y_hat"]
	return impact.transpose(..., "obs_ind")

My rationale for thinking that the model expectation would be the best to focus on was pretty simple. Basically, the frequentist approach to causal inference focusses on model expectations. And in the Bayesian approach we can simply express our posterior beliefs about what the model expectation is.

Of course, that doesn't nullify the utility of the posterior predictive. It will just tell you something else. But if the focus is on estimating causal effects, then my thought was to focus on the model expectation.

Happy to get more input and views on this.