Why is there "residual dots" in plot() of GAM? After building a generalized additive model (GAM) using mgcv package, we can use the plot function to visualize the smoother, like:
plot(M1, resid = TRUE)

The resid = TRUE argument ensures that residuals are added to the figure.
My question is: why this plot (to show the shape of smoother) is related to RESIDUAL? I mean, what we are interested should be the relationship between smoother and response variable, not residual. Then why residual occur in this smoother's plot?
In my immature understanding, "residual" is used in plots of model validation, e.g., residuals versus fitted values, residuals versus covariates.
 A: If you look at the help for the plot.gam function you'll see that the residuals argument indicates that what is added to the plots of the smooth functions are partial residuals, not the model residuals.
The estimated partial residual of the $i$th observation for the $j$th smooth function in the model, $\hat{\varepsilon}^{\text{partial}}_{ji}$, are given by (following Wood (2017, p.184))
$$\hat{\varepsilon}^{\text{partial}}_{ji} = \hat{f_{j}}(x_{ij}) + \hat{\varepsilon}^{\text{p}}_i ,$$
where $\hat{\varepsilon}^{\text{p}}_i$ is the estimated Pearson residual for the $i$th observation from the full model.
The reason we might want to add partial residuals to the partial effects plot of the smooth function is that these partial residuals should be evenly scattered about the smooth function if the model is well fitted to the data. They act as another diagnostic plot to assess the model fit.
References
Wood, S.N., 2017. Generalized Additive Models: An Introduction with R, Second Edition. CRC Press.
