Here is a similar plot I adapted from Clementini et al. (2019). It shows mean absolute Gaia G-magnitude (derived from Gaia-DR2 parallaxes) versus period for 998 known Cepheid variables. I guess that corresponds reasonably well with what you have done.
The dashed line was added by Clementini et al. They believe stars below this line are either not Cepheids or there is a problem with the Gaia parallaxes or the published periods. Interestingly, you have few of this population in your plot? So either that's an improvement with the Gaia DR3 data (I assume this is what you are using) or a cleaner parent catalogue of Cepheids.
I have added the red line (by hand) - this appears to be the red line plotted in your diagrams. Clearly this is not a good fit, although I would say your gradient look approximately the same. The reason why this gradient does not agree with the literature you quote is that the plot here uses absolute G-magnitude (which is not the same as the V-band) and there has been no attempt to remove extinction from the sample, which will certainly increase scatter and decrease the intercept value by the average extinction of an object in the sample. It will also bias the gradient since the more luminous Cepheids tend to be further away and perhaps have more extinction - this will artificially flatten (reduce) the gradient. This effect is clearly seen in Fig.6 of Clementini et al. where they compare the P-L relationship of galactic Cepheids over the whole sky, that are affected by varying levels of extinction, with a P-L relation for Cepheids in the Large Magellanic Cloud that have a close-to-uniform level of extinction. The P-L relation for the latter is clearly steeper.
It is a mystery however as to why there seems to be an offset between the data you show and the data plotted in Clementini et al. (2019). Are you using the intensity-weighted mean G magnitude (known as int_average_g in the Gaia Cepheid tables) ?
