Sorting and spillovers can create correlation in individual outcomes. In this situation, standard discrete choice estimators cannot consistently estimate the probability of joint and conditional events, and alternative estimators can yield incoherent statistical models or intractable estimators. I propose a random effects estimator that models the dependence among the unobserved heterogeneity of individuals in the same cluster using a parametric copula. This estimator makes it possible to compute joint and conditional probabilities of the outcome variable, and is statistically coherent. I describe its properties, establishing its efficiency relative to standard random effects estimators, and propose a specification test for the copula. The likelihood function for each cluster is an integral whose dimension equals the size of the cluster, which may require high-dimensional numerical integration. To overcome the curse of dimensionality from which methods like Monte Carlo integration suffer, I propose an algorithm that works for Archimedean copulas. I illustrate this approach by analysing labour supply in married couples.