We discuss heavy quarkonium production through parton fragmentation, including a review of arguments for the factorization of high-p_T particles into fragmentation functions for hadronic initial states. We investigate the further factorization of fragmentation functions in the NRQCD formalism, and argue that this requires a modification of NRQCD octet production matrix elements to include nonabelian phases, which makes them gauge invariant. We describe the calculation of uncanceled infrared divergences in fragmentation functions that must be factorized at NNLO, and verify that they are absorbed into the new, gauge invariant matrix elements.