Why does my visreg plot not match my partial effect plot for GAM? Here's my code (looking at only one variable):
gam <- gam(heart_data$DEATH_EVENT~ s(age) + anaemia + 
           s(creatinine_phosphokinase) + diabetes +
           s(ejection_fraction) + high_blood_pressure + 
           s(platelets) + s(serum_creatinine) + 
           s(serum_sodium) + sex + smoking, family = "binomial", 
           data = heart_data, method = "REML")

plot(gam, select = c(5), trans = plogis, rug = TRUE, 
     residuals = TRUE,
     pch = 1, cex = 1, shift = coef(gam)[1],  seWithMean = TRUE)
visreg(gam, "serum_creatinine", scale="response", rug=2, 
       xlab="serum creatinine",
       ylab="P(mortality)")

I get:


I assumed that after correcting for the intercept they should match. What am I missing?

The visreg plot is conditional upon some reference values for the other covariates, which take their median values by default.
The plot produced by plot.gam is a partial plot, showing only the change in $E(Y)$ as a function of the specified covariate. In effect, this plot is totally ignoring the effects of the other covariates on the response (the effect of the covariate of interest on $E(Y)$ was estimated conditional upon the other terms in the model of course). All you are really doing here is rescaling the y-axis to put the original partial effect plot on a more natural scale.
Another way to think about the difference is that the visreg() plot shows the results of predict(gam, newdata = foo) where foo is a data frame containing a vector of values for the covariate of interest and constant values (their median value) for all other terms in the model. It is the effects of — and the additional uncertainties thereof — these additional covariates that shift the curve up or down and alter the credible interval displayed.
