In this paper, a computationally efficient approach is proposed for the determination of the nonstationary response statistics of hysteretic oscillators endowed with fractional derivative elements. This problem is of particular practical significance since many important engineering systems exhibit hysteretic/inelastic behavior optimally captured only through the concept of fractional derivative, and many natural excitations as seismic waves and atmospheric turbulence are both stochastic and nonstationary in time. Specifically, the approach is based on a statistical linearization scheme involving an equivalent system of augmented dimension. First, relying on a transformation scheme, the fractional derivative term is represented by a set of coupled linear ordinary differential equations. Next, the evolution of the system response statistics is captured by incorporating the statistical linearization technique in a nonstationary sense. This involves integrating in time a set of ordinary differential equations. Several numerical applications pertaining to classical hysteretic oscillators are considered, and the versatility of the proposed method is assessed via comparison with pertinent Monte Carlo simulations.