For multilayered or coated substrates in elastohydrodynamic-lubricated (EHL) contacts, the subsurface stress distributions under a normal load combined with shear traction have been analyzed in this article through computer simulations. The Papkovich-Neuber potentials and Fourier transform are adopted to deduce the pressure-displacement, pressure-stress, and shear traction-stress response functions in frequency domain for the coated substrates, and to calculate distributions of pressure and subsurface stress. The results from the analysis of EHL contacts on coated substrates are compared with those from dry contact model in which shear traction is assumed to obey Coulomb's law. Effects of the Young's modulus of coatings, the properties of lubricants, and the magnitude of traction are discussed. Similar to the results in dry contacts, hard coatings in lubricated cases tend to increase the von Mises stress, whereas soft coatings decrease the stress. Shear traction makes the max von Mises stress increasing and moving closer to surface. However, the changes in subsurface stress due to shear traction are less obvious in lubricated contacts. Comparison between EHL and dry contact models reveals that lubrication can reduce the von Mises stress in the coating layer due to smaller shear traction. The analyses show that pressure, film thickness, and subsurface stress distributions are influenced by surface coatings, sliding velocity, rheological models, and pressure-viscosity behaviors.