I am running a regression with multiple interaction terms and multiple fixed effects. I am looking to plot the slope (derivative) of a variable, conditional upon it being interacted or not. I came across the marginaleffects package, which seems to be able to do this. However, it did not accept the felm model I was using from the lfe package, so I switched to using feols from the fixest package instead. My regression looks like the following:
reg <- feols(data = df, outcome ~ interact_dummy * (IV_Dummy1 + control1 + control2) + interact_dummy * (IV_Dummy2 + control1 + control2)| country_fe + year_fe)
I then plot attempt to plot the slope using the following
plot_slopes(reg, variables = "IV_Dummy1", condition = "interact_dummy")
But the plot is the following:

I would like instead to plot these points connected, showing the increase in the slope caused by the interaction effect. Ideally, this could also have shaded regions for the confidence interval (as the current plot has bounds). Is this issue specific to using a fixed effects model, or the fixest package? Is there an alternative package to use, or have I specified the plot_slopes() incorrectly
Additionally, the estimates from the feols regression are the same as when I ran the regression in the felm model, however their significance levels are slightly different. What accounts for this?
EDIT:
(Sample code)
reg2 <- feols(data = mtcars, mpg ~ vs * (am + hp + carb) + vs * (disp + hp + carb) | gear + carb)
plot_slopes(reg2, variables = "am", condition = "vs")
