Disclaimer: My expertise in hierarchical Bayesian models is also intermediate at best.
It looks to me like the response categories do not have an ordering. This is unfortunate from a methodological perspective, since an ordered categorical response is much simpler to model (and the results simpler to interpret) than if no ordering exists.
Your problem seems to be exactly what multinomial regression is for. Confusingly, these models are called by different names in different quarters – sometimes multinomial, sometimes polytomous, sometimes categorical.
Anyway, I think a hierarchical (mixed in frequentist terms) multinomial (categorical in brms terminology) logit or probit model in the correct approach. In the simplest approach, the responses would be a single vector (‘factor’ in R-speak) containing all the nominal-scale responses occurring in the data. Subject id would indeed constitute the group-level (random in frequentist terms) effect(s), and the uh…test conditions/stimuli would be the population-level (fixed in frequentist terms) effects. At least the brms package fits these models just fine. The results of such models are a little complex to interpret due to what you state above, i.e. the dependence of every category probability on every other category probability. The Nobel prize will surely go to whoever can come up with a solution to this. I’m not holding my breath. Meanwhile, third-party libraries exist to aid in the interpretation of these models.
Computational costs decrease if the data are grouped such that every distinct combination of explanatory variable values (including Subject) are grouped so that the multinomial response is a C-column matrix whose C columns represent the counts in each response category for each distinct combination of explanatory variable values. brms can do this for binomial outcomes, but I don’t know if this format is supported for categorical responses yet. Either way, the results for group and population-level effects should be identical in both approaches.