Incorporating multiple variables in Hierarchial model

Questions v5
1

Hi All, I am trying to build a Hierarchial model. I am new to Pymc and baysian so lot of learnings. All the examples I see Hierarchial models are using just single variable and one group. I have 4 different variable. x1 to x4, two groups group1, group2. I am expecting that variable x1 and x2 will have group specific interactions, while rest of variable x3 and x4 should have one common slope.

So in terms of building model. I have this

df.columns = ['group1', 'group2 , 'x1', 'x2', 'x3', 'x4', y]
n_unique_group1 = 10
n_unique_group2= 5
coords = {"group1": group1}

with pm.Model(coords=coords) as model:

    group1 = pm.MutableData('sku_cd', df.group1, dims='obs_id')
    x1 = pm.MutableData('x1', df.x1, dims= 'obs_id')
    x2 = pm.MutableData('x2', df.x2, dims='obs_id')
    x3 = pm.MutableData('x3', df.x3, dims= 'obs_id')
    x4 = pm.MutableData('x4', df.x4, dims = 'obs_id')
    
    # Priors
    alpha_mu = pm.Normal('alpha_mu', mu = 0.0, sigma = 10)
    alpha_sigma = pm.Exponential('alpha_sigma', 10)

    beta1_mu = pm.Normal('beta1_mu', mu = 0, sigma = 10)
    beta1_sigma = pm.Exponential('beta1_sigma', 1)

    beta2_mu = pm.Normal('beta2_mu', mu = 0, sigma = 10)
    beta2_sigma = pm.Exponential('beta2_sigma', 1)
    
    beta3_mu = pm.Normal('beta3_mu', mu = 0, sigma = 10)
    beta3_sigma = pm.Exponential('beta3_sigma', 1)  
    
    beta4_mu = pm.Normal('beta4_mu', mu = 0, sigma = 10)
    beta4_sigma = pm.Exponential('beta4_sigma', 1)

    # Random Intercept
    z_a = pm.Normal("z_a", mu=0, sigma=1, dims="group1")
    alpha = pm.Deterministic("alpha", alpha_mu + z_a * alpha_sigma, dims="group1")

    z_b1 = pm.Normal("z_b1", mu=0, sigma=1, dims="code")
    beta1 = pm.Deterministic("beta1", beta1_mu + z_b1 * beta1_sigma, dims="group1")

    z_b2 = pm.Normal("z_b2", mu=0, sigma=1, dims="code")
    beta2 = pm.Deterministic("beta2", beta2_mu + z_b2 * beta2_sigma, dims="group1")

    # Slopes for population level variables
    beta3 = pm.Normal("beta3", mu=beta3_mu, sigma=beta3_sigma)
    beta4 = pm.Normal("beta4", mu=beta4_mu, sigma=beta3_sigma)

    # Model Error
    sigma_y = pm.Exponential('sigma_y', 1)


    y_hat = alpha[group1] \
        + beta1[group1] * x1 \
            + beta2[group1] * x2 \
                + beta3* x3 \
                    + beta4 * x4 \

    # Likelihood

    y_like = pm.Normal('y_like', mu = y_hat, sigma = sigma_y, observed= df.y, dims='obs_id')

My Question is is that above formulation looks good? If not please do direct me in a right direction.
thank you.

I also get ton of convergence warnings, which I solved by non-centered approaches. Still a not of convergence issues, funnel shape for beta3 and beta4 the population specific slopes. One thing to note is my have like intotal around 3K data points, it is smaller size.