The ability of density functional theory to predict the relative energies of different spin states, especially for systems containing transition metal atoms, is of great importance for many applications. Here, in order to sort out the key factors determining accuracy, we compare the predictions of 60 density functional approximations of 10 different types [local spin density approximation, generalized gradient approximation (GGA), nonseparable gradient approximation (NGA), global-hybrid GGA, range-separated hybrid GGA, range-separated hybrid GGA plus molecular mechanics, meta-GGA, meta-NGA, global-hybrid meta-GGA, and range-separated hybrid meta-GGA] for their ability to represent the spin-flip transitions of all 4d transition metal atoms of groups 3-10 (Y through Pd) and their singly positive cations. We consider all 16 excitation energies connecting the ground states (of the neutral atoms and the cations) to their first excited states of different multiplicities, and we also consider all eight ionization potentials. We also test the Hartree-Fock method. All density functional and Hartree-Fock calculations are converged to a stable solution, in which the spatial symmetry is allowed to be completely broken to achieve the lowest possible energy solution. By analyzing the fractional subshell occupancies and spin contaminations, we are able to sort out the effects of s orbital vs d orbital bias and high-spin vs low-spin bias. A reliable functional should have little or no bias of either type rather than succeeding for a limited subset of cases by cancellation of errors. We find that the widely used correlations of spin splittings to percentage of Hartree-Fock exchange are not borne out by the data, and the correlation functionals also play a significant role. We eventually conclude that SOGGA11-X, B1LYP, B3V5LYP, and MPW3LYP are the most consistently reliable functionals for balanced treatments of 4d transition metal atoms and their cations.