A new technique for Monte Carlo sampling of the hard-sphere collision force has been applied to study the interaction between a pair of spherical macroions in primitive-model electrolyte solutions with valences 1:2, 2:1, and 2:2. Macroions of the same charge can attract each other in the presence of divalent counterions, in analogy with earlier observations for planar and cylindrical geometries. The attraction is most significant at intermediate counterion concentrations. In contrast to the entropic depletion force between neutral particles, attraction between macroions is of energetic origin. The entropic contribution to the potential of mean force is generally repulsive at conditions corresponding to aqueous colloids with or without salt. For systems with divalent counterions, the potentials of mean force predicted by mean-field approximations like the Derjaguin-Landau-Verwey-Overbeek (DLVO) theory or the Sogami-Ise (SI) theory are qualitatively different from those observed in the simulations. However, for systems with monovalent counterions, predictions of DLVO theory are in fair agreement with simulation results.