ValueError: shapes not aligned in basic model

Questions
1

I have the following bayesian hierarchical model, where y is of dimension (1313x1), X (1313, 2) and Z (1313x2).

y = data["Score"]
y = y.values.reshape(-1, 1)

X = data[["Var1", "Var2"]].values

Z = data[["Control1", "Control2"]].values

with pm.Model() as model:

    gamma = pm.Normal("gamma", mu=0, sigma=1, shape=(2, 2))

    u = pm.Normal("u", mu=0, sigma=1, shape=2)

    beta = pm.math.dot(Z, gamma) + u

    # Model error
    r = pm.Gamma("r", alpha=9, beta=4, shape=())

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

    # Likelihood
    y_like = pm.Normal("y_like", mu=y_hat, sigma=r, observed=y)

with model:
    trace = pm.sample(draws=500, chains=2, tune=2000, discard_tuned_samples=True, random_seed=SEED, progressbar=True,
                      cores=1, init="adapt_diag", return_inferencedata=True)

If I run the model, the following error message comes up:

ValueError: shapes (1313,2) and (1313,2) not aligned: 2 (dim 1) != 1313 (dim 0)

I considered transposing beta from (1313x2) to (2, 1313) but I am not sure whether its shape is correct at all.
However this gave me the following error

ValueError: Mass matrix contains zeros on the diagonal.
The derivative of RV u.ravel()[1] is zero.
The derivative of RV gamma.ravel()[2] is zero.
The derivative of RV gamma.ravel()[3] is zero.

I use pymc3 3.11.2

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.

3

That is not exactly what I want to model.
Anyway, if I use Bambi, the model formula is quite cumbersome, as I need to write the model as a string, without using matrix algebra, somehow like this:

model = bmb.Model("y ~ var_1*covar_1 + var_2*cov_2 + ... + var_n*covar_k", data, family="gaussian")

Or is there also a possibility to use dot products, like

model = bmb.Model(y = np.dot(X, Z), data, family=“gaussian”)

?

4

I’m not sure why you are including terms like "var_1*covar_1" in the model, which adds a particular type of interaction effect in Bambi. Is that what you want?

On the other hand, if you have a response Y and three predictors X_1, X_2, and X_3, no matter their nature, you would do

bmb.Model("y ~ x1 + x2 + x3", data)

without having to multiply the predictors by their coefficients (because these are the ones estimated with the model).

If you have many more predictors, you could do something like

predictors = ["x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8", "x9", "x10"]
formula = "y" + " + ".join(predictors)

bmb.Model(formula, data)