There are going to be several different ways to handle this, but it’s going to impact both your model and your hierarchical_normal() function. The simplest way would be to just do this kind of thing:
age = data_df.age_num.astype(int).values
age_idx,ages = pd.factorize(data_df["age_num"])
n_age = len(ages)
region_idx,regions = pd.factorize(data_df["region"])
n_region = len(regions)
and then index into your parameter arrays like this:
η = β0 + α_age[age_idx] + α_region[region_idx]
However, you can definitely utilize pm.MutableData objects if the flexibility would be of use to you. I would also recommend utilizing the dimension/coordinate functionality to have more meaningful index labels. Check here and here.