I have generated 10 model summaries after fitting the same model on 10 different subsets of the dataset which is as follows

library(mice)
data("nhanes")
head(nhanes)
imp <- mice(nhanes, print = FALSE, m = 10, seed = 24415)
df <- complete(imp, "long")
model_fit <- lapply(1:10, function(i) {
  model = lm(bmi ~ age + hyp + chl,
             data = subset(df, `.imp`==i))
})

From this I get different ggpredict objects

ggpredict(model_fit[[1]], c("age", "hyp"))
ggpredict(model_fit[[2]], c("age", "hyp"))
ggpredict(model_fit[[3]], c("age", "hyp"))
ggpredict(model_fit[[4]], c("age", "hyp"))
ggpredict(model_fit[[5]], c("age", "hyp"))
ggpredict(model_fit[[6]], c("age", "hyp"))
ggpredict(model_fit[[7]], c("age", "hyp"))
ggpredict(model_fit[[8]], c("age", "hyp"))
ggpredict(model_fit[[9]], c("age", "hyp"))
ggpredict(model_fit[[10]], c("age", "hyp"))

I am looking for an efficient way to a) Estimate the average of all the ggpredict objects b) Plot this using ggplot function.

So far I tried storing the results from each ggpredict function as list object and

`Reduce(`+`, list_ggpred)/length(list_ggpred)`

I got warning,

 " In Ops.factor(left, right) : `+1 not meaningful for factors.

Any suggestions highly appreciated. Thanks.