I’m working with a hierarchical model and I have some conceptual doubts about how to declare hierarchies.

I have hierarchies that are product categories (e.g., food, beverages) and product families (e.g., within food → fruits, chocolates; within beverages → water, soft drinks). Therefore, the fruit family exists only within the food category, and has no relationship to the beverages category.

The way I’m stating the hierarchies today is as follows:

```
idx_category, category = df["category"].factorize()
idx_family, family = df["family"].factorize()
coords = {
"obs_id": idx_category,
"category": category,
"family": family,
}
with pm.Model(coords=coords) as model:
price = pm.MutableData("price", df["log_price"].values, dims="obs_id")
# price influence
a_global = pm.Normal("a_global", mu=5, sigma=10)
sigma_a_cat = pm.HalfNormal("sigma_a_cat", sigma=2)
a_cat = pm.Normal("a_cat", mu=a_global, sigma=sigma_a_cat, dims="category")
sigma_a_fam = pm.HalfNormal("sigma_a_fam", sigma=2, dims="category")
a_fam = pm.Normal("a_fam", mu=a_cat, sigma=sigma_a_fam, dims=("family", "category"))
# intercept
T0_global = pm.Normal("T0_global", mu=10, sigma=10)
sigma_T0_cat = pm.HalfNormal("sigma_T0_cat", sigma=2)
T0_cat = pm.Normal("T0_cat", mu=T0_global, sigma=sigma_T0_cat, dims="category")
sigma_T0_fam = pm.HalfNormal("sigma_T0_fam", sigma=2, dims="category")
T0_fam = pm.Normal("T0_fam", mu=T0_cat, sigma=sigma_T0_fam, dims=("family", "category"))
mu = T0_fam[idx_family, idx_category] + a_fam[idx_family, idx_category] * p
Y_obs = pm.Exponential("Y_obs", lam=1/mu, observed=df["y_obs"].values)
```

However, when I use this approach, because of `dims=("family", "category")`

, I end up with parameters like `a_fam[fruits, beverages]`

- which don’t make sense. Sampling is taking a long time to run, and my hypothesis here is that this explosion of unnecessary parameters is affecting performance.

What would be the best way to declare hierarchies in this situation?

Bonus question: In this model, I only have hierarchies on the `mu`

parameters for each `pm.Normal`

- i.e., `a_cat`

`mu`

’s comes from `a_global`

and `a_fam`

`mu`

’s comes from `a_cat`

- but `sigma_a_cat`

and `sigma_a_fam`

are independent of the hierarchical leveling. Is this the right approach or should I also create a hierarchy of sigmas? My fear here is that with independent sigmas, things deviate a lot from the higher hierarchy.