An implementation of time-correlator theory to model resonance Raman scattering in systems where thermal effects are important is presented. This treatment is applicable to strong coupling and high-temperature cases where the first-order approximation is inaccurate. To illustrate important effects of thermalization in molecules with Raman-active low-frequency vibrations or in molecules at high temperatures, calculations are presented for a four-mode model molecule with low-frequency modes and large electron-phonon coupling (total S = 2). For such large coupling, an increase in temperature causes a nearly symmetric broadening of the absorption spectrum, an overall decrease in Stokes scattering intensities for both low- and high-frequency modes and complex changes in anti-Stokes intensities. Experimental studies of the protein systems rhodopsin and bacteriorhodopsin (bR) are then modeled to illustrate that such effects are observed in real molecular systems. Analysis of the resonance Raman intensities of rhodopsin is presented for a fully harmonic model of the excited-state potential surface. Using the resulting parameters, the contribution of displaced low-frequency modes to a Gaussian homogeneous electronic absorption line shape is examined quantitatively, and both the first-order approximation and an effective zero-temperature approximation using scaled excited-state displacements are found to be inadequate for calculating the Raman intensities of these modes.