So, you can still use indexing to do this, but it will be quite inefficient as you need to create a index of direction + block, for example, index 0 would be parallel block 1.
A better approach is to formulate a design matrix:
from patsy import dmatrices
outcome, predictors = dmatrices("projection ~ block + direction + block * direction", data)
which is what usually call mixed effect model or hierarchical model. The pymc3 doc have quite an extensive session on this (e.g., https://docs.pymc.io/notebooks/GLM-hierarchical.html)
The differences is that, here you your case the linear prediction will goes into sigma instead of mu