The use of the proportional odds (PO) model for ordinal regression is ubiquitous in the literature. If the assumption of parallel lines does not hold for the data, then an alternative is to specify a non-proportional odds (NPO) model, where the regression parameters are allowed to vary depending on the level of the response. However, it is often difficult to fit these models, and challenges regarding model choice and fitting are further compounded if there are a large number of explanatory variables. We make two contributions towards tackling these issues: firstly, we develop a Bayesian method for fitting these models, that ensures the stochastic ordering conditions hold for an arbitrary finite range of the explanatory variables, allowing NPO models to be fitted to any observed data set. Secondly, we use reversible-jump Markov chain Monte Carlo to allow the model to choose between PO and NPO structures for each explanatory variable, and show how variable selection can be incorporated. These methods can be adapted for any monotonic increasing link functions. We illustrate the utility of these approaches on novel data from a longitudinal study of individual-level risk factors affecting body condition score in a dog population in Zenzele, South Africa.