Hi sorry for note getting to you earlier,
I think such a model would be most easily expressed using the brms package, which supports both ordinal outcomes (using the cumulative family, see https://journals.sagepub.com/doi/10.1177/2515245918823199 for more background) and fixed correlation matrices (Estimating Phylogenetic Multilevel Models with brms • brms).
I’ll however note that with such a big dataset, fitting with brms (or any MCMC method) might be computationally very demanding. You are however likely to also be able to fit such a model with R-INLA, which doesn’t use MCMC and uses an approximation instead and is thus quite a bit faster, while in high-data regime the approximation tends to be very good.
p-values and likelihood-ratio tests don’t have a direct counterpart in the Bayesian paradigm, so I am not sure you will obtain what you need :-) For my current best thinking on how to answer questions about model/hypothesis choice see Hypothesis testing, model selection, model comparison - some thoughts which lists some ways to express such questions when fitting Bayesian models.
Best of luck with your model!