# Unexpected behavior with arviz.plot\_hdi() with categorical x

**URL:** <https://discourse.pymc.io/t/unexpected-behavior-with-arviz-plot-hdi-with-categorical-x/16403>\
**Category:** General\
**Tags:** bambi, arviz\
**Created:** [January 19, 2025, 12:07am UTC](https://discourse.pymc.io/t/unexpected-behavior-with-arviz-plot-hdi-with-categorical-x/16403 "2025-01-19T00:07:14Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![milesalanmoore](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/milesalanmoore/32/8921_2.png) [@milesalanmoore](https://discourse.pymc.io/u/milesalanmoore)\
**Post date:** [January 19, 2025, 12:07am UTC](https://discourse.pymc.io/t/unexpected-behavior-with-arviz-plot-hdi-with-categorical-x/16403/1 "2025-01-19T00:07:14Z")

</div>

Hi all!

As a big Stan and `brms` user, I am really enjoying learning to fit Bayesian models in Python with the PyMC ecosystem.

While attempting to refit models from my own work in `bambi`, I ran into some trouble plotting posterior expectations and 95% HDIs from an ANOVA model. Here’s a simple simulation to demonstrate the issue

```Python
import arviz as az
import bambi as bmb
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# Simulate data
np.random.seed(42)
x = ['A', 'B', 'C']
yA = np.random.normal(loc=5, scale=3, size=30)
yB = np.random.normal(loc=2, scale=4, size=30)
yC = np.random.normal(loc=7, scale=1.8, size=30)

# Create a DataFrame
data = pd.DataFrame({
    'y': np.concatenate([yA, yB, yC]),
    'group': np.repeat(x, repeats=30)
})
data['group'] = data['group'].astype('category')

# Plot the data
plt.figure(figsize=(8, 6))
sns.boxplot(x='group', y='y', data=data, palette='Set3')
sns.stripplot(x='group', y='y', data=data, color='black', alpha=0.5, jitter=True)
plt.title('Distribution of y across groups')
plt.xlabel('Group')
plt.ylabel('y')
plt.show()

```

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/5/58e31f692e57981d1a295a13965a0d8f2018ded4.png)

```Python
# Fit a Bayesian ANOVA model using Bambi
model = bmb.Model('y ~ group', data)
idata = model.fit()

preds = model.predict(idata, kind="response_params", inplace=False)
y_mu = az.extract(preds["posterior"])["mu"].values
group = data.loc[:, "group"].values

```

Now I want to use az.plot\_hdi() to plot the expectation of the posterior distribution that I extracted (in `preds` above).

Following the documention for `arviz`, I first tried:

```Python

az.plot_hdi(x=group, y=y_mu.T)

#UFuncTypeError: ufunc 'multiply' did not contain a loop with signature matching types (dtype('<U1'), dtype('float64')) -> None

```

It turns out that `az.plot_hdi` has a parameter called _smooth_ that is **True** by default. In this case, the function attempts to interpolate across all `x`, which is not possible for `numpy` to do when `x` is categorical. Setting _smooth_ to **False** works around the error, but does not produce the plot I’d expect it to (interpolating the HDI across groups instead of plotting, say, segments for the HDI region):

```Python
az.plot_hdi(x=group, y=y_mu.T, smooth=False)

```

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/4/4f1d518a47f9589c1f71ff8146173c7f949143be.png)

Am I missing something with respect to how `az.plot_hdi()` is intended to behave? The documentation doesn’t seem to suggest the function should only work if x is numeric, but that seems to be the behavior? Perhaps there’s a plot\_kwarg that I should be passing through to `matplotlib` to get a more fitted ANOVA like plot?

If this is not expected, I am happy to try and figure out a solution and contribute to the repo. I wanted to open up a discussion before a github issue though.

---

<div class="post-metadata">

**Author:** ![tcapretto](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/tcapretto/32/4902_2.png) [@tcapretto](https://discourse.pymc.io/u/tcapretto)\
**Post date:** [January 19, 2025, 1:04pm UTC](https://discourse.pymc.io/t/unexpected-behavior-with-arviz-plot-hdi-with-categorical-x/16403/2 "2025-01-19T13:04:19Z")

</div>

Hi and welcome!

I don’t think it’s possible to modify the behavior of `arviz.plot_hdi()` without modifying its internals. However, you can still do what you want with the `interpret` sub-package in Bambi.

```python
bmb.interpret.plot_predictions(
    model=model,
    idata=idata,
    conditional="group",
);

```

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/c/cf4e84517eabc9b2b531197c3efb77d2cadd95ef.png)

That displays the 94% HDI for the mean parameter. If you want it for the predictive distribution, you can do

```python
bmb.interpret.plot_predictions(
    model=model,
    idata=idata,
    conditional={"group": ["A", "B", "C"]},
    pps=True
);

```

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/e/ea092547173dec082f6fcb53e031423f2204bd21.png)

You can see more examples in [Tools to interpret model outputs](https://bambinos.github.io/bambi/notebooks/#tools-to-interpret-model-outputs).

* * *

It should be possible to map group levels to colors but there’s a bug right now. I’m opening an issue in the Bambi repository to fix it.

**Edit** issue here [`bambi.interpret_plot_predictions()` fails when we condition and color by the same categorical variable · Issue #870 · bambinos/bambi](https://github.com/bambinos/bambi/issues/870)

---

<div class="post-metadata">

**Author:** ![milesalanmoore](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/milesalanmoore/32/8921_2.png) [@milesalanmoore](https://discourse.pymc.io/u/milesalanmoore)\
**Post date:** [January 19, 2025, 2:36pm UTC](https://discourse.pymc.io/t/unexpected-behavior-with-arviz-plot-hdi-with-categorical-x/16403/3 "2025-01-19T14:36:34Z")

</div>

Thank you, @tcapretto! This makes sense to me. [Peeking at the function definition](https://github.com/arviz-devs/arviz/blob/main/arviz/plots/hdiplot.py), I see what you mean about needing to modify internals. Thank you for providing an alternative solution using the `interpret` sub module. I will play around some with this! I think I should do some reading about plotting models fit with `PyMC` and `bambi`. I am much less comfortable here than I am with `brms` and `bayesplot`’s add\_pred\_draws() and `ggplot2` / `ggdist` workflow in R. Though I am excited about how lively the `PyMC` community is!

If you’d like, I can open an issue on `arviz` regarding the `az.plot_hdi` function. I think a ValueError if categorical or str type x’s are supplied would be more informative than the error that is currently returned from `np.interpolate` when this is the case (when `smooth=True`). This is something I am more than happy to submit a PR for as well.

-Miles

* * *

**Edit** : Issue and PR created: [[WIP] `plot_hdi` raise exception when `x` is string (#2412) by milesalanmoore · Pull Request #2413 · arviz-devs/arviz · GitHub](https://github.com/arviz-devs/arviz/pull/2413)

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [January 22, 2025, 8:46am UTC](https://discourse.pymc.io/t/unexpected-behavior-with-arviz-plot-hdi-with-categorical-x/16403/4 "2025-01-22T08:46:32Z")

</div>

11 posts were split to a new topic: [User experience: Python vs R, PyMC vs Stan vs JAX](https://discourse.pymc.io/t/user-experience-python-vs-r-pymc-vs-stan-vs-jax/16426)