GLMM when the dependent variable has three levels I'm trying to fit mixed models in R for data where the dependent variable is categorical with three levels (choices of True, False, Unsure). My plan was to fit three separate binomial mixed models (for True vs. False-or-Unsure, False vs. True-or-Unsure, Unsure vs. True-or-False) using glmer in lme4, and then draw conclusions about the impact of the fixed effects on each choice by looking at the separate models. But I have some questions about whether this is the best approach:
Question 1: Is the correct way to implement this method to fit two separate binomial  models (e.g. True vs. False, True vs. Unsure) instead?
Question 2: After fitting separate binomial models, does one need to somehow combine them into a single model in order to draw conclusions about the influence of different effects on the dependent variable? Or is it fine to (say) use emmeans for post-hoc comparisons on each of the binomial models separately?
Question 3: Is it important for the random and fixed effects to be the same for each binomial model? Having applied likelihood ratio tests to various candidates for each of the binomial models, it seems that the random effects that produce the best models differ in each case. E.g. the addition of random slopes for workers produces a significant improvement in the model for False vs. True-or-Unsure but not in the one for Unsure vs. True-or-False. So is deciding on a single set of effects part of the model-selection process, or can there be methodological grounds to expect different random effects for each of the binomial models (e.g. if workers' strategies differ when choosing False but not for other choices)?
Question 4: Are there reasons to instead fit a multinomial mixed model via other packages in R? I've heard that the ordinal and MCMCglmm packages allow this. My concern was that interpreting and analysing such a model would be more challenging for someone who's new to statistics.
Update 1:
In light of the comments and answer below, it looks like my original plan of fitting separate binomial models isn't great. I've tried using ordinal::clmm() and mclogit::mblogit(), but have had some problems with both options, which I summarise below.
ordinal:

*

*The model converges using the code clmm(choice ~ A * B + (1|item) + (B|participant), data = targets). But I'm struggling to understand the coefficients in the summary. E.g., the estimate for level a2 of fixed effect A reads -4.5560; does this mean that a2 decreases the likelihood of a 'higher' choice? Even if this is right, I don't know how to interpret this with respect to my hypotheses. Some of my hypotheses concern the likelihood of Unsure (coded as 2) relative to False (1) or True (3); so what is the impact of a decreased likelihood of a 'higher' choice?


*I'm also confused about the threshold coefficients. The one for 1|2 reads -4.2262; so does this mean that there is a decreased likelihood of choice 1 relative to 2, calculated with respect to the model's reference level (which is presumably level a1 of effect A with level b1 of effect B)?


*The application of emmeans works fine. Presumably, the estimate for a1 b1-a2 b1 (= 4.556) indicates that there is an increased likelihood of a 'higher' choice for a1 b1 relative to a2 b1. But I again struggle to relate this to my hypotheses concerning the likelihood of Unsure.
mclogit:

*

*There are some problems with getting models that converge, perhaps due to incorrect syntax. The code mblogit(choice ~ A * B, data = targets, random = ~1|participants) works fine, but attempting to add a random slope or multiple random effects produces error messages. E.g., mblogit(choice ~ A * B, data = targets, random = ~1|item + ~1|participants) gives the error 'subscript out of bounds', and mblogit(choice ~ A * B, data = targets, random = ~B|participants) gives the error 'variable 'B' is absent, its contrast will be ignored' plus a failure to converge.


*The summary of the model that converges is easier to interpret than for clmm, since there are separate parts for 2 vs 1 and for 3 vs 1. E.g. for 2 vs 1, the estimate for a2 reads  -1.4687; this presumably means that a2 decreases the likelihood of choice Unsure (relative to False), which can be directly linked to my hypotheses about the likelihood of Unsure.


*Something goes wrong with the application of emmeans to the model that converges: all of the estimates are very small and all of the p-values = 1.
Main question: Should I use ordinal or mclogit, and how do I resolve the above queries for the advised package?
Update 2:
I've been pursuing the strategy of using mclogit::mblogit(). Most models now converge without errors when I use the code  mblogit(choice ~ A * B, data = targets, random = list(~1|item,  ~B|participants), method = "MQL". I have two remaining problems: 1. attempts at comparing models by using likelihood ratio tests via anova() don't go smoothly (unexpected degrees of freedom are often reported and the p-value is sometimes absent); 2. applying emmeans to models continues to produce strange results (the p-values are all listed as either 1 or NaN). If anyone is familiar with mclogit, then please tell me how to fix these problems (I can give more details) or whether there is a better way to compare models and carry out post-hoc comparisons.

I would not recommend fitting 3 seperate models. This will result in a loss of statistical power, and the the models themselves will not be directly comparable since the random structures will not be the same. In  generalised linear mixed models, the fixed effects are condtional on the random effects.
The mclogit can fit multinomial logit models in their conditional logit and baseline logit variants with random effects:
https://cran.r-project.org/web/packages/mclogit/index.html
