Plotting Slopes of interaction effect from Fixed effects model with marginaleffects package (R)

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")

The difference in significance between lfe and fixest is probably due to the fact that fixest reports robust standard errors by default. See this vignette: https://lrberge.github.io/fixest/articles/standard_errors.html

The reason why the plot is point-range by default, instead of a continuous line is that the moderating variable is binary, and because it rarely makes sense to plot continuous lines between binary categories (it implies that we can interpolate).

If that is still what you want to do, it is trivial to set the draw argument to FALSE, extract the raw data, and plot it with ggplot2. For example:

library(fixest)
library(marginaleffects)
library(ggplot2)

mod <- feols(hp ~ vs * am | carb, data = mtcars)

# extract raw data for plotting
dat <- plot_slopes(mod, variables = "vs", condition = "am", draw = FALSE)

# convert categorical (binary/factor) variable to numeric
dat$am <- as.numeric(as.character(dat$am))

ggplot(dat, aes(x = am, y = estimate, ymin = conf.low, ymax = conf.high)) +
    geom_line() +
    geom_ribbon(alpha = .2) +
    theme_minimal() +
    labs(y = "Slopes")