This was very helpful and solved my issue. I ultimately realized that I had to include my ‘cluster_index’ within X and remove the dimension from \delta.
β = pm.Normal("β", mu=δ, sigma=3., dims= ("clusters", "num_cols"))
μ = (X[cluster_index, :] * β[cluster_index, :]).sum(axis=1)
I need to deepen my understanding of dimensions, but this is providing me with momentum. Thanks!