[eng] The analysis of heavy quark-antiquark bound states -also called quarkonium states- produced in high-energy collisions offers a challenging opportunity to test our understanding of the Quantum Chromodynamics (QCD), the theory of strong interactions. Production of quarkonium states involves low-energy non-perturbative effects that are inherent to the dynamics of the bound state. The creation of the heavy quarks, owing to their large mass, also implies an energy scale at which the dynamics of strong interactions can be analyzed as a perturbation by virtue of the asymptotic freedom of QCD. Non-Relativistic QCD (NRQCD) provides a theoretical framework that disentangles the high- and low-energy scales in quarkonium production. It can be used to factorize out the non-perturbative effects into a restricted number of process-independent parameters in the expression of the rate, leaving a process-dependent part that can be computed within perturbative QCD. In this thesis, we investigate the implications of the NRQCD factorization on the phenomenology of quarkonium production in electron-positron annihilation, in photoproduction and in hadronic collisions. In particular, we study the impact of the radiative and the relativistic corrections on specific observables that have been measured. By comparing our new more precise predictions to the experimental data we emphasize that some features of quarkonium production still need to be elucidated. To gain further insight on the mechanisms at work in quarkonium production, we propose an innovative approach that involves extending the set of observables analyzed so far. We present a new algorithm based on Monte-Carlo techniques aimed at analyzing the detailed kinematics of quarkonium events. We consider new signatures, such as the hadroproduction of a quarkonium state in association with a heavy-quark pair of the same flavor, showing how these measurements will be helpful in improving our understanding of quarkonium production.