Conversation
|
Check out this pull request on See visual diffs & provide feedback on Jupyter Notebooks. Powered by ReviewNB |
PR SummaryIntroduce maketables plug-in dunder methods on BaseExperiment (coef table, stats, depvar, vcov) with Bayesian/OLS handling, plus a demo notebook.
Written by Cursor Bugbot for commit c5bb622. This will update automatically on new commits. Configure here. |
| widths = sorted_samples[interval_size:] - sorted_samples[:n_intervals] | ||
| min_idx = np.argmin(widths) | ||
| ci_lower = float(sorted_samples[min_idx]) | ||
| ci_upper = float(sorted_samples[min_idx + interval_size]) |
There was a problem hiding this comment.
HDI calculation has off-by-one error in interval bounds
The HDI (Highest Density Interval) calculation in _maketables_coef_table_bayesian has an off-by-one error. For a 94% interval containing ceil(0.94 * n) samples starting at index i, the interval spans indices i to i + interval_size - 1. However, the code computes widths = sorted_samples[interval_size:] - sorted_samples[:n_intervals], which calculates sorted_samples[i + interval_size] - sorted_samples[i] instead of sorted_samples[i + interval_size - 1] - sorted_samples[i]. Similarly, ci_upper uses sorted_samples[min_idx + interval_size] when it needs sorted_samples[min_idx + interval_size - 1]. This results in intervals containing one extra sample and potentially selecting a suboptimal HDI.
| } | ||
| # Add R² if available | ||
| if hasattr(self, "score") and self.score is not None: | ||
| stat_map["r2"] = self.score |
There was a problem hiding this comment.
R² statistic returns Series instead of scalar value
The __maketables_stat__ method assigns self.score directly to stat_map["r2"] without handling the case where self.score is a pandas Series (which occurs for PyMC models). As shown in the notebook output, this causes the R² row in maketables to display the entire Series representation (unit_0_r2 0.836121 unit_0_r2_std 0.012656 dtype: float64) instead of a clean scalar value. Other code in the codebase (e.g., _bayesian_plot) properly extracts self.score["unit_0_r2"] when self.score is a Series.
| interval_size = int(np.ceil(0.94 * n)) | ||
| n_intervals = n - interval_size | ||
| widths = sorted_samples[interval_size:] - sorted_samples[:n_intervals] | ||
| min_idx = np.argmin(widths) |
There was a problem hiding this comment.
HDI calculation crashes with small MCMC sample sizes
The HDI calculation in _maketables_coef_table_bayesian crashes with ValueError: attempt to get argmin of an empty sequence when the number of MCMC samples is 16 or fewer. When n <= 16, interval_size = ceil(0.94 * n) equals n, making n_intervals = 0, which results in an empty widths array. Calling np.argmin on an empty array raises a ValueError. While unlikely in normal usage where MCMC runs with hundreds of samples, this could occur in quick testing scenarios with minimal sampling parameters.
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #600 +/- ##
==========================================
- Coverage 93.27% 92.60% -0.67%
==========================================
Files 37 37
Lines 5632 5682 +50
Branches 367 373 +6
==========================================
+ Hits 5253 5262 +9
- Misses 248 289 +41
Partials 131 131 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
Hi @drbenvincent ,
Attached a very drafty PR (with lots of support from Claude) to enable maketables support for CausalPy via the new Plug-In solution @dsliwka implemented.
Attached is also an example notebook.
For the RDD class, you e.g. get something like this:

Would you generally be interested in merging this PR?
Where I'd need your help - what coefficients should be reported by default? Which statistics might users be interested in beyond the defaults? Is my choice of 94% probability intervals a good one? Etc.
Note that this PR is more of a proof of concept and might need some cleaning up and unit tests. If this makes it into CausalPy, we should likely also implement tests over there to ensure we never break CausalPy flows.
Best, Alex
📚 Documentation preview 📚: https://causalpy--600.org.readthedocs.build/en/600/