Peak ground fractional order responses (PGRα) as a generalization of conventional seismic intensity measures (IMs) such as peak ground acceleration (PGA) and peak ground velocity (PGV) can better predict the seismic performance of structural systems. This paper proposes the first ground motion prediction equation (GMPE) for PGRα for active shallow crustal regions using a subset of the PEER NGA-West2 ground motion database. The model development database consists of 4491 accelerograms from 82 different earthquakes in California with the magnitude and rupture distance ranges of Mw 4.0–7.9 and RRUP 0–300 km, respectively. PGRα intensity measures are computed from the modified Oustaloup’s recursive approximation to Caputo’s definition of differintegral operator. The main functional form of the predictive model is decided by implementing statistical ground motion data-driven testing methods such as the likelihood approach and Euclidean distance concept. The final functional form of the predictive model accounts for magnitude, distance, style-of-faulting, linear and nonlinear site, hanging wall, basin response, and anelastic distance attenuation effects, and models the aleatory variability with respect to Mw and VS30. The final predictive model produces PGA (α = 0), PGV (α = −1), and peak ground fractional order responses at 19 different α values ranging from −0.05 to −0.95 for the average horizontal component. The proposed predictive model draws estimates of ground motion amplitudes that are consistent with those from the NGA-West2 models for PGA and PGV for sample earthquake scenarios. Moreover, it can offer a basis for predictive modeling of peak ground fractional order response quantities for performance assessment of structures and infrastructures across a region.