Our objective was to evaluate the efficacy of robenacoxib in osteoarthritic dogs using four ordinal responses measured repeatedly over time. We propose a multivariate probit mixed effects model to describe the joint evolution of endpoints and to evidence the intrinsic correlations between responses that are not due to treatment effect. Maximum likelihood computation is intractable within reasonable time frames. We therefore use a pairwise modeling approach in combination with a stochastic EM algorithm. Multidimensional ordinal responses with longitudinal measurements are a common feature in clinical trials. However, the standard methods for data analysis use unidimensional models, resulting in a loss of information. Our methodology provides substantially greater insight than these methods for the evaluation of treatment effects and shows a good performance at low computational cost. We thus believe that it could be used in routine practice to optimize the evaluation of treatment efficacy.