Partial residuals plot I am trying to understand how the gam package in R generates the partial residuals plots, so I tried to create one from scratch to compare to the one generated by plot.gam(), but they don't match. I wanted to try as simple an example as possible, so I am only using one variable, $lstat$, as a predictor of $medv$, in the Boston housing prices dataset. 
library(MASS)
library(gam)
lm.fit.Boston.1 <- lm(medv ~ lstat, data = Boston)
summary(lm.fit.Boston.1)
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.55384    0.56263   61.41   <2e-16 ***
## lstat       -0.95005    0.03873  -24.53   <2e-16 ***

Then plot the partial residuals:
plot.gam(lm.fit.Boston.1, se=TRUE)
grid()


If I am understanding Introduction to Statistical Learning, p. 285 (footnote at bottom), correctly, the partial residual plot is the best-fit line of $lstat$ as a predictor of the partial residuals, in this case $r_i=y_i-\beta_0$.
beta_0 <- coef(lm.fit.Boston.1)[1]
y <- Boston$medv
r <- y - beta_0
lm.residuals <- lm(r ~ Boston$lstat)
plot(Boston$lstat, lm.residuals$fitted.values)
grid()


What am I doing wrong, why aren't the plots identical? Thank you.
 A: You don't compute the partial residuals correctly.
As the footnote in Section 7.7.1 GAMs for Regression Problems of ISLR explains:

Partial residuals for $X_3$ have the form $r_i = y_i - f_1(x_{i1}) - f_2(x_{i2})$.

To compute the partial residuals for $X_3$ we estimate $f_1$ and $f_2$ by fitting a (GAM) model for $Y$ on $X_1$ and $X_2$. We don't use $X_3$ in this "partial" model. In your example there is only one predictor, lstat. So the intercept 1 plays the role of $X_1$ and $X_2$ and lstat — the role of $X_3$.
What you did wrong: You used the full model medv ~ 1 + lstat to compute partial residuals for lstat. 
What you should do instead: Use the partial model medv ~ 1.
[ISLR] G. James, D. Witten, T. Hastie, and R. Tibshirani. An Introduction to Statistical Learning with Applications in R. Springer, 2nd edition, 2021. Available online.
Here are the partial residuals plots. To make it easier to compare plots, I use the same x:y ratio.
lm.full    <- lm(medv ~ 1 + lstat)
lm.partial <- lm(medv ~ 1)

# This is the plot we want to reproduce by hand
make_gam_plot(lm.full)


# This is the wrong way to compute partial residuals for `lstat`
# as we use `lstat` in the model which produces the residuals
resid <- medv - coef(lm.full)["(Intercept)"]
lm.full.residuals <- lm(resid ~ lstat)

make_base_plot(lstat, fitted(lm.full.residuals))


# This is the right way to compute partial residuals for `lstat`
resid <- medv - coef(lm.partial)["(Intercept)"]
lm.partial.residuals <- lm(resid ~ lstat)

make_base_plot(lstat, fitted(lm.partial.residuals))


Created on 2022-07-10 by the reprex package (v2.0.1)
R code to reproduce the figures:
attach(MASS::Boston)

xlim <- c(0, 40)
rlim <- c(-30, 20)

make_gam_plot <- function(model) {
  gam::plot.Gam(
    model,
    se = TRUE,
    xlim = xlim,
    ylim = rlim
  )
  grid()
}

make_base_plot <- function(x, y) {
  plot(x, y,
    xlim = xlim,
    ylim = rlim
  )
  grid()
}

lm.full <- lm(medv ~ 1 + lstat)
lm.partial <- lm(medv ~ 1)

make_gam_plot(
  lm.full
)

resid <- medv - coef(lm.full)["(Intercept)"]
lm.full.residuals <- lm(resid ~ lstat)

make_base_plot(
  lstat, fitted(lm.full.residuals)
)

resid <- medv - coef(lm.partial)["(Intercept)"]
lm.partial.residuals <- lm(resid ~ lstat)

make_base_plot(
  lstat, fitted(lm.partial.residuals)
)

