ValueError: shapes not aligned in basic model

Questions
2

Considering that y is of dimension (1313x1), X is (1313, 2), and Z is (1313x2), we can see that

beta = Z * gamma + u results in an object of shape (1313, 2). Then, when you do the following:

y_hat = pm.math.dot(X, beta)

You’re trying to compute the dot product of a matrix of shape (1313x2) with a matrix of shape (1313x2), which is not possible.

Following classic hierarchical models notation, I think that what you want is something like

with pm.Model() as model:
  beta = pm.Normal("beta", mu=0, sigma=1, shape=(2, 1))
  u = pm.Normal("u", mu=0, sigma=1, shape=(2, 1))
  r = pm.Gamma("r", alpha=9, beta=4)

  y_hat = pm.math.dot(X, beta) + pm.math.dot(Z, u)
  y_like = pm.Normal("y_like", mu=y_hat, sigma=r, observed=y)

Edit: For models like this, consider using Bambi, which makes your hierarchical GLM modelling much easier.