# Gam plots y-axis range is different for model that that from predict function

I got some idea about gam (generalized additive model) plot from a course of Noam Ross. However, I have several confusion and got stucked. I fit a logistic gam and can interpret by using trans=plogis. But when interpreting in terms of link function, I don't know how to interpret them. The plot shows y-axis range that includes zero even if when there are no zero in data.

Also the range in y-axis is different than that of a plot using predict function.

Again, predicted value from predict function using terms and link also are different. Thus plot using these two also have different y-axis value.

What could be the reason behind it? Is there any way to have the actual y-axis value in terms of link function?

This is my code:

library(mgcv)
library(gss)
data(wesdr)
attach(wesdr)

g = gam(ret~s(dur)+s(gly)+s(bmi),family=binomial,method="REML")

p= predict(g, type="terms")  # predicted value for three variables using "terms"

par(mfrow=c(1,3))

plot(g,select=1) #smooth plot for first variable from model
plot(dur,p[,1])  #smooth plot for first variable using type 'terms'
plot(dur,p.s)  #smooth plot for first variable using type 'link'



All these plots have different yaxis!

Note that you have the plots in the wrong order in the included figure; panel 2 should be panel 3 and vice versa.

The first plot is the result from evaluating the smooth function $$f_j(\mathtt{dur}_i)$$ at 100 evenly spaced values over the range of dur. But as it also include the confidence interval, which is very wide for large values of dur, the y axis range is much larger than the values taken by the smooth function itself.

That this axis includes 0 is simply due to the sum-to-zero constraints imposed on the smooths. The 0 on the axis is the expected value of the response on the link scale; in other words, 0 represents the model constant (intercept) term, and values above 0 on the plot are larger than the expected value, etc. We need to do things like this so we can have an intercept in the model.

The third plot, the result of the call predict(g, type="terms") is showing you essentially the same thing as the first plot; the different axis range is simply the result of you not showing the uncertainty in the estimated partial effects of $$f_j(\mathtt{dur}_i)$$ via the confidence band. As the default is to return these term-wise contributions on the link scale, and for the same reason as above, the 0 on the axis is the result of the sum-to-zero constraint applied to the smooth.

The second plot, the result of the predict() call while trying to exclude certain terms is not what you think it is. {mgcv} is a little bit naughty here; the variables you tried to exclude from the prediction aren't being excluded because there are no terms in the model with the names you supplied. Here you'd need to provide the names of the smooths as they are shown in summary(g). So, what your plot showed is simply the fitted values for wesdr data set observations, but on the link scale. As the link scale is the logit scale, it is perfectly fine for the values to be negative (technically they could range from minus- to plus-infinity, but because of the way the logit works, they are most likely to be in the range +/-4). If you were to apply the inverse of the link function to the predicted values, you'd get values between 0 and 1.

The reason the values are all over the place (and not on a nice smooth line) is because you didn't actually exclude the effects of $$f_j(\mathtt{gly}_i)$$ and $$f_j(\mathtt{bmi}_i)$$ from the predictions.

What you needed to do is to call:

p.s <- predict(g, exclude = c("s(gly)","s(bmi)"), type="link")


Note that both your second and third plots are not going to result in something like plot 1 because you are predicting for the observed values of dur, which are not evenly spaced or in ascending order in the observed data frame. Hence if you were to plot the fitted values using a line, you'd get a horrible mess too. If you really want to do this stuff by hand, you want to create a data slice along dur holding gly and bmi at fixed representative values, and then use predict(..., type = "terms") (and then apply the inverse of the link function to map them to the 0,1 scale if that's what you want) or predict(g, exclude = c("s(gly)","s(bmi)"), type="response") I find this easier with my gratia package:

ds <- data_slice(g, dur = evenly(dur))
fv <- fitted_values(g, data = ds, exclude = c("s(gly)","s(bmi)"),
scale = "response")
fv


This results in

> fv
# A tibble: 100 × 7
dur   gly   bmi fitted    se lower upper
<dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>
1  1.2   12.2  22.9  0.262 0.285 0.169 0.383
2  1.75  12.2  22.9  0.278 0.253 0.190 0.387
3  2.29  12.2  22.9  0.295 0.222 0.213 0.393
4  2.84  12.2  22.9  0.313 0.196 0.237 0.400
5  3.38  12.2  22.9  0.330 0.173 0.260 0.409
6  3.93  12.2  22.9  0.348 0.156 0.282 0.421
7  4.47  12.2  22.9  0.366 0.144 0.303 0.434
8  5.02  12.2  22.9  0.383 0.136 0.323 0.448
9  5.56  12.2  22.9  0.400 0.131 0.340 0.463
10  6.11  12.2  22.9  0.416 0.128 0.356 0.478
# ℹ 90 more rows
# ℹ Use print(n = ...) to see more rows


and note that fitted-values() did the right thing and computed the fitted values and their confidence band on the link scale and then backtransformed the fitted values and the upper and lower bounds of the band back to the response scale.

Now you can plot this:

library("ggplot2")

fv |>
ggplot(aes(x = dur, y = fitted)) +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2) +
geom_line()


which gives

In general, you were on the right track but just missed a couple of important details:

1. the plot produced by plot(g) (or gratia::draw(g)) for the $$f_j(\mathtt{dur}_i)$$ term uses 100 evenly spaced values over the range of dur and it is drawn on the link scale, and as the smooth is subjected to the sum-to-zero identifiability constraint 0 values are perfectly OK, and
2. you didn't provide the correct term labels to exclude the effects of the other smooths when you predicted.