The standard way is to use a linear regression like approach. What you do is fine, but usually in a GLM context, you modify your group indicator so that:
|GROUP| --> Design matrix
| A | | 1 | 0 | 0 |
| B | | 1 | 1 | 0 |
| A | | 1 | 0 | 0 |
| C | | 1 | 0 | 1 |
And the coefficient will naturally coded for the contrast between each combination of groups. Usually, doing so you get better conditioned matrix and inference is usually easier as well (as the intercept coded for the overall mean, all the rest of the coefficients are usually in the same scale). This is also what the glm component in pymc3 is doing (using patsy internally to generate the design matrix)