The principle of maximum conformality (PMC) provides a way to eliminate the conventional renormalization scale ambiguity in a systematic way. By applying the PMC scale setting, all non-conformal terms in a perturbative series are summed into the running coupling, and one obtains a unique, scale-fixed prediction at any finite order. In this study, we make a detailed PMC analysis for the spin-singlet heavy quarkoniums decay (into light hadrons) at the next-to-leading order. After applying the PMC scale setting, the decay widths for all those cases are almost independent of the initial renormalization scales. The PMC scales for η_c and h_c decays are below 1 GeV; to achieve a confidential pQCD estimation, we adopt several low-energy running coupling models to carry out the estimation. By taking the MPT model, we obtain Γ(η_c →LH)=25.09^+5.52_−4.28 MeV,Γ(η_b →LH)=14.34^+0.92_−0.84 MeV, Γ(h_c →LH)=0.54^+0.06_−0.04 MeV and Γ(h_b →LH)=39.89^+0.28_−0.46 keV, where the errors are calculated by taking m_c∈[1.40 GeV, 1.60 GeV] and m_b∈[4.50 GeV, 4.70 GeV]. These decay widths agree with the principle of minimum sensitivity estimations, in which the decay widths of η_c,b are also consistent with the measured ones.