Using coefficients instead of the predict function to make predictions for a GAM in R

Question

I have a question in regards to using the coefficients from a GAM which uses smooth variables.

I want to know how to do so manually, i.e. not just plugging into the predict function.

I tried an example, but just by looking at the coefficients for the smoothed terms I immediately knew that plugging them into the equation would not yield the same values as the predict function.

Here is my example, which comes from the Hitters dataset in ISLR2:

library(ISLR2)
library(mgcv)

### One more model, for stackexchange
gam_stack=gam(Walks~s(RBI, k=3)+s(Hits, k=5), data=Hitters, method="REML")

### Our new data:
new_data=data.frame(RBI=c(25), Hits=c(110))

### The prediction using predict:
test_pred=predict(gam_stack, new_data)
test_pred
       1 
37.15598 

 ### getting the coefficients from the model:
 coefs_gam_stack=gam_stack$coefficients
 coefs_gam_stack
(Intercept)    s(RBI).1    s(RBI).2   s(Hits).1   s(Hits).2   s(Hits).3   s(Hits).4 
  38.742236    6.844722    5.544233    6.384631   -2.667766   15.174375    2.232861 

I want to use the coefficients extracted, coefs_gam_stack above, to predict walks, the way one could for a linear model. Simply looking at the coefficients it is easy to see that simply plugging in the values in the new_data vector above will not yield the same prediction as using the predict function. The predict function for the new data gives a predicted number of walks of 37.15598, which is test_pred in my code. Obviously looking at the coefficients, which are found in coefs_gam_stack above, simply plugging in the values in new_data won't yield the same result as the predict function.

Obviously I'm doing something wrong, and as I am somewhat new to GAM modelling it is clearly my incomplete understanding of both the methodology and the output of the gam function in R.

Any help would be greatly appreciated, obviously you can't just plug the coefficients and new_data values in the way one can with a linear model. I suspect what I am missing has something to do with knots, but I'm not quite sure.

I have seen similar stackexchange questions answered using the lpmatrix from the predict function, but this seems only to get a matrix of the data used to build the model, whereas I want to be able to manually predict with new data.

Thanks in advance! I'm hoping this will really solidify my understanding of GAM modelling. Right now I can still use it, and use the predict function, but I want to bolster my knowledge of just what is going on underneath the hood.

Answer 1

predict.gam is doing essentially three things:

1. evaluating the basis functions for the spline(s) in your model at the values of the covariates supplied to it (or the observed data if you don't provide anything to the newdata argument) to create what is called the linear predictor matrix, denoted by $\mathbf{X}_p$,
2. computing the result of the operation $\boldsymbol{\hat{\eta}}_p = \mathbf{X}_p \boldsymbol{\hat{\beta}}$, i.e. computing the linear predictor values for the supplied data, and
3. applying the inverse of the link function, $g^{-1}$, to yield predicted values on the response scale.

As @Roland mentions in their comment, you can work through these steps manually as:

library("mgcv")
set.seed(1)
df <- gamSim()
m <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = df, method = "REML")

# step 1.
Xp <- predict(m, newdata = df, type = "lpmatrix")

# step 2.
eta_p <- drop(Xp %*% coef(m))

# step 3.
y_hat <- family(m)$linkinv(eta_p)

Compare:

head(y_hat)
head(pred_gam <- predict(m, newdata = df, type = "response"))

all.equal(y_hat, pred_gam, check.attributes = FALSE)

yielding

> head(y_hat)
  head(pred_gam <- predict(m, newdata = df, type = "response"))

  all.equal(y_hat, pred_gam, check.attributes = FALSE)
        1         2         3         4         5         6
 5.895958  3.148708  8.281588  8.650655 15.683875  8.384466
        1         2         3         4         5         6
 5.895958  3.148708  8.281588  8.650655 15.683875  8.384466
[1] "target is numeric, current is array"

The only differences are because of how Simon returns predicted values from predict.gam(), as an array, not a vector (I suppress some of these differences with check.attributes = FALSE.)

If you want to get even deeper into what predict.gam(..., type = "lpmatrix") is doing, you can use mgcv::smoothCon() and mgcv::PredictMat(), which are used to create a smooth and the evaluate the basis function of a smooth at supplied covariate values. Note that this works on a single smooth at a time, so you'd need to combine the outputs from four calls to PredictMat() to get the lpmatrix outputs, plus add on a constant term.

If you want to go deeper than that, you'll need to learn the math behind the spline type you are using and study Simon's C code.