Plotting GAM with Robust Clustered Std. Errors

Question

So, I have a GAM with fixed effects and I want to plot it with robust clustered standard errors.

Basically, I have a similar model:

gam1 <- gam(iris$Sepal.Length ~ s(Petal.Width) + factor(Species), data = iris)

where I have, in this case, fixed effects for Species. Now, I wanted a plot similar to:

plot(gam1, select = 1, seWithMean = T, shift = coef(gam_m1)[1],  shade = TRUE)

But with robust standard errors clustered by Species.

---

I already tried things similar to:

robust_se <- sqrt(diag(vcovCL(gam1, cluster = iris$Species)))

# Calculate predicted values
# predDF <- ggpredict(gam1, terms = "ideo")

predDF <- predict.gam(gam1, se = T)

predDF$u_ci <- predDF$fit + 1.96 * robust_se
predDF$l_ci <- predDF$fit - 1.96 * robust_se

predDF$Sepal.Length <- predDF$fit
predDF$Petal.Width <-  iris$Petal.Width

predDF$Species <- iris$Species

predDF <- data.frame(predDF)

ggplot(predDF, aes(x = Petal.Width, y = Sepal.Length)) +
  geom_smooth(method = "gam", formula = y ~ s(x),
              method.args = list(cluster = predDF$Species))
              se = T) +
  geom_ribbon(aes(ymin = l_ci, ymax = u_ci), fill = "blue", alpha = 0.2) 
                    

but the plot is all wrong (and tbh I'm not even entirely sure I'm doing it right).

Thanks in advance! :)