We discuss factorization in heavy quarkonium production in high energy collisions using NRQCD. Infrared divergences at NNLO are not matched by conventional NRQCD matrix elements. However, we show that gauge invariance and factorization require that conventional NRQCD production matrix elements be modified to include Wilson lines or non-abelian gauge links. With this modification NRQCD factorization for heavy quarkonium production is restored at NNLO.