I’m working with a dataset that has missing values, and I am using multiple imputation with the mice package in R. I then fit a Generalized Additive Model (GAM) using the mgcv package on the imputed datasets. I would like to plot the overall fit of the model across all imputations, but I’m not sure how to combine the plots from each imputed dataset.

library(mgcv)
library(mice)
library(ggplot2)

# Set seed for reproducibility
set.seed(123)

# Generate synthetic dataset
n <- 200
data <- data.frame(
  outcome = rnorm(n),
  exposure = rnorm(n),
  covariate1 = sample(0:1, n, replace = TRUE),
  covariate2 = rnorm(n, 40, 2),
  covariate3 = rnorm(n, 30, 1),
  covariate4 = sample(0:1, n, replace = TRUE),
  covariate5 = sample(0:1, n, replace = TRUE)
)
data[sample(1:n, 50), "exposure"] <- NA
data[sample(1:n, 50), "covariate1"] <- NA

# Perform multiple imputation
imp <- mice(data, m = 5, method = 'pmm', seed = 123)
model <- with(imp, gam(outcome ~ s(exposure) + covariate1 + covariate2 + covariate3 + covariate4 + covariate5, method = 'REML'))
# Plot
fitted_models <- model$analyses
par(mfrow = c(2, 3))  # Adjust as needed to fit all plots
for (i in 1:length(fitted_models)) {
  plot(fitted_models[[i]], shade = TRUE, rug = TRUE, residuals = TRUE, pch = 1, cex = 1, main = paste("Imputation", i))
}

This does not appear to be a solution

# pooled.fit <- pool(model)
# plot.gam(pooled.fit)