library(marginaleffects)
dat = palmerpenguins::penguins
mod = lm(bill_length_mm ~ flipper_length_mm * species, data = dat)
pre = predictions(
    mod,
    newdata = datagrid(
        flipper_length_mm = fivenum,
        species = unique
    )
)
library(plot2)
par(las = 1)
 with(
     pre,
     plot2(
         x = flipper_length_mm, y = estimate,
         ymin = conf.low, ymax = conf.high,
         by = species,
         type = "ribbon",
         palette = "harmonic",
         frame = FALSE, grid = TRUE
     )
 )