The development and use of statistical models for categorical outcomes is a common task in the social sciences, e.g. when analyzing large-scale survey data. These models allow the scientist to identify explanatory variables that are significant for the different choices encoded in the response. When the data are hierarchically structured, multilevel modeling through e.g. mixed-effects is customary to avoid structure-induced bias. In this project, we develop a Bayesian multilevel model to account for group-specific effects when the outcome of the experiment is assumed to be categorical, i.e. a one-hot-encoding vector that selects one category out of K possibilities. The Bayesian inference is done via Gibbs sampling with Polya-gamma data augmentation, and we further perform variable selection for all effects in the hierarchical model (both random and fixed-effects), using spike-and-slab priors. We perform simulation tests on synthetic data, and will present one real-world application using survey data.