We study the large-order asymptotic behavior of the perturbation series for short-distance coefficients in the NRQCD factorization formulas for the decays J/ψ → e + e − and η c → γγ. The short-distance coefficients of the leading matrix elements are calculated to all orders in the large-N f limit. We find that there is a universal Borel resummable renormalon associated with the cancellation of the Coulomb singularity in the short-distance coefficients. We verify that the ambiguities in the short-distance coefficients from the first infrared renormalon are canceled by ambiguities in the nonperturbative NRQCD matrix elements that contribute through relative order v 2 . Our results are used to estimate the coefficients of higher order radiative corrections in the decay rates for J/ψ → e + e − , η c → γγ, and J/ψ → γγγ.