Creating a Hierarchical Ordinal Logistic Model

Hi all,

I’m working on a problem where I’m trying to create ordinal regression models to predict ratings in multiple different contexts. However, I have reason to believe that these various contexts are all related, and so I’ve been trying to create a hierarchical model. The code I’ve created thus far is,

n_contexts = len(np.unique(context_codes))

with pm.Model() as model:
    
    # for weights
    mu_w = pm.Normal("mu_w", mu=0, sd=1)
    sigma_w = pm.HalfCauchy("sigma_w", beta=1)
    
    # for thresholds
    mu_thresh = pm.distributions.MvNormal(
        "mu_thresh",
        mu=np.array([-2, -1, 0, 1]),
        cov=np.eye(4),
        shape=4
    )
    sigma_thresh = pm.HalfCauchy("sigma_thresh", beta=2)
    
    # weight distribution
    w = pm.Normal("w", mu=mu_w, sd=sigma_w, shape=n_contexts) 
    
    # thresholds distribution
    thresholds = []
    for i in range(n_contexts):
        thresholds.append(
            pm.Normal(
                "thresholds_{}".format(i),
                mu=mu_thresh,
                sd=sigma_thresh,
                shape=4,
                transform=pm.distributions.transforms.ordered
            )
        )
    
    wx = w[context_codes] * x
    y_ = pm.OrderedLogistic(
        "y", 
        cutpoints=tt.stack(thresholds)[context_codes],
        eta=wx,
        observed=y
    )

When I try to create this model though, I get the error, ValueError: all the input array dimensions except for the concatenation axis must match exactly.

The model works when I assume fixed thresholds for all contexts, so I suspect that this has something to do with the way the cutpoints are being set, but I don’t know exactly what’s causing the error. Any help would be really appreciated.