4

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:

enter image description here enter image description here

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

kjetil b halvorsen
  • 77,844
Nonya
  • 337
  • 1
  • 10

1 Answers1

4

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.

Gavin Simpson
  • 47,626
  • 1
    Could you maybe explain in more detail how visreg produces the plot? Like, the maths behind it. I thought it was the same/similar to avPlot but I'm a bit confused now about the "predict(gam, newdata = foo)". Also thank you for the answer – Nonya Jul 24 '21 at 01:04