A Bayesian sampling algorithm for parameter estimation in a discrete-response model is presented, where the dependent variables contain two layers of binary choices and one ordered response. The investigation is motivated by an empirical study using such a double-selection rule for three labour-market outcomes, namely labour-force participation, employment and occupational skill level. It is of particular interest to measure the marginal effects of some mental health factors on these labour-market outcomes. The contribution is to present a sampling algorithm, which is a hybrid of Gibbs and Metropolis-Hastings algorithms. In Monte Carlo simulations, numerical maximization of likelihood fails to converge for more than half of the simulated samples. The proposed Bayesian method represents a substantial improvement: it converges in every sample, and performs with similar or better precision than maximum likelihood. The proposed sampling algorithm is applied to the double-selection model of labour-force participation, employment and occupational skill level, where marginal effects of explanatory variables, in particular the mental health factors, on the three labour-force outcomes are assessed through 95 Bayesian credible intervals. The proposed sampling algorithm can easily be modified for other multivariate nonlinear models that involve selectivity and are difficult to estimate by other means.