Decomposition-order dependence of time development integrator on ensemble averages for the Nosé-Hoover dynamics is discussed. Six integrators were employed for comparison, which were extensions of the velocity-Verlet or position-Verlet algorithm. Molecular dynamics simulations by these integrators were performed for liquid-argon systems with several different time steps and system sizes. The obtained ensemble averages of temperature and potential energy were shifted from correct values depending on the integrators. These shifts increased in proportion to the square of the time step. Furthermore, the shifts could not be removed by increasing the number of argon atoms. We show the origin of these ensemble-average shifts analytically. Our discussion can be applied not only to the liquid-argon system but also to all MD simulations with the Nosé-Hoover thermostat. Our recommended integrators among the six integrators are presented to obtain correct ensemble averages.