Thanks for the clarification. First off: that notebook rendering is really cool! Second, the xarray tricks you showcase are pretty handy. To the extent these notebooks can serve as something of a rosetta stone going from stan/R to pymc3 I imagine others might be interested in getting a better handle on all the coordinate/broadcasting magic. I ended doing something for option (1) by putting the data into pandas, but it feels a bit clunky.
pm_data = az.from_pymc3(
trace=trace, coords={"treatment": data["treatment"], "actor": data["actor"]},
)
samples = pd.DataFrame(np.vstack(pm_data.posterior["p"]).T)
samples["treatment"] = pm_data.constant_data.treat_id
samples["actor"] = pm_data.constant_data.actor_id
samples_long_form = samples.melt(id_vars=["treatment", "actor"])
samples_long_form.groupby(["treatment", "actor"]).mean()