Hierarchical model, categorical outcome

Questions
1

Hi there,

I am new and I am stuck. :slight_smile: Sorry for these probably trivial questions. Any tips are much appreciated! Below my model code but here are also the data and the full analysis (simple model, and hierarchical model).

I have 3 categorical outcome variables (participants’ choice), 2 predictors; one binary (S), one continuous, standardised (Es). I tested participants, and I have multiple observations by participants and items.

So far I was able to implement a hierarchical model that accounts for differences between subjects (varying intercepts and slopes). At least this is what I think I did :D, please correct me if I’m wrong.

I have two questions about the model:

  1. How can I use vectorised parameters given the hierarchical structure?
    I tried to make the code cleaner by using vectorised parameters (one beta, and then calculate mu with pm.math.dot(X, beta). But I always get the shapes wrong.

  2. How to account for repeated measures by item? I should also account for the repeated measures by item but I am also afraid to make my model too complex if I were to attempt partial pooling for items as well. I don’t know whether this is a valid concern. Sampling already takes a while. Is there any neat way?

Since this is my first attempt, there must be other issues. Please let me know if that is the case! Thank you!

Here is the model code:

with pm.Model() as hier_mod:

mu_a = pm.Normal('mu_a', mu=0.0, sd=1.0, shape=2)
sigma_a = pm.HalfNormal('sigma_a', 1, shape=2)
mu_b1 = pm.Normal('mu_b1', mu=0., sd=.5, shape=(1,2))
mu_b2 = pm.Normal('mu_b2', mu=0., sd=.5, shape=(1,2))

intercept = pm.Normal('intercept', mu=mu_a, sd=sigma_a, shape=(n_sub, 2))
intercept_f = tt.concatenate([np.zeros((n_sub, 1)), intercept], axis=1)
                              
beta_1 = pm.Normal('beta_1', mu=mu_b1, sd=.5, shape=(n_sub, 2))
beta_1f = tt.concatenate([np.zeros((n_sub, 1)), beta_1], axis=1)  

beta_2 = pm.Normal('beta_2', mu=mu_b2, sd=.5, shape=(n_sub, 2))
beta_2f = tt.concatenate([np.zeros((n_sub, 1)), beta_2], axis=1)                            

mu = intercept_f[sub_idx] + beta_1f[sub_idx] * X_Es + beta_2f[sub_idx] * X_S
likelihood = pm.Deterministic('likelihood', tt.nnet.softmax(mu))

y1 = pm.Categorical('y1', p=likelihood, observed=y_r)

trace_hier = pm.sample(1000, tune=2000, cores=1)

FYI: n_sub = number of subjects, sub_idx = subject ID

This is it for now, thank you again everyone!
Alex