I’ve recently worked with stan again a little after mostly using PyMC for years, and the thing I missed most was (perhaps surprisingly) support for named dimensions. It might sound like a pretty minor detail, but to me it really makes a difference in practice.
And that is even though the pymc support is hardly perfect. The random variables remember the dimension names, but all intermediate variables don’t and just use numpy broadcasting rules.
If I were to redesign something, I think that would be pretty far up my list of things I want to support.