In categorical data analysis, there is rich literature for modeling binary and polychotomous responses. However, existing methods are inadequate for handling combinatorial responses, where each response is an array of integers subject to additional constraints. Such data are increasingly common in modern applications, such as surveys collected under skip logic, event propagation on a network, and observed matching in ecology. Ignoring the combinatorial structure in the response data may lead to biased estimation and prediction. The fundamental challenge for modeling these integer-vector data is the lack of a link function that connects a linear or functional predictor with a probability respecting the combinatorial constraints. In this paper, we propose a novel augmented likelihood, in which a combinatorial response can be viewed as a deterministic transform of a continuous latent variable. We specify the transform as the maximizer of integer linear program, and characterize useful properties such as dual thresholding representation. When taking a Bayesian approach and considering a multivariate normal distribution for the latent variable, our method becomes a direct generalization to the celebrated probit data augmentation, and enjoys straightforward computation via Gibbs sampler. We provide theoretical justification for the proposed method at an interesting intersection between duality and probability distribution and develop useful sufficient conditions that guarantee the applicability of our method. We demonstrate the effectiveness of our method through simulation studies and a real data application on modeling the formation of seasonal matching between waterfowl.