[R] mgcv: I can't manually reconstruct a P-spline from a GAM's coefficients

[Alexander Engelhardt]

Hello R-helpers,

I am working through Simon Wood's GAM book and want to specify my own 
knot locations (on even tens, i.e. 10, 20, 30, etc.). Then, I want to 
compute a GAM on that area, and given the coefficients, reconstruct the 
same P-spline that is drawn in plot(my_gam).

I'm failing.
Here is what I am doing. I think the two plot calls at the end should 
produce the same line, but they don't:


bspline <- function(x,k,i,m=2){
     ## "draws" one B-spline basis function
     ## order of splines is m+1, i.e. m=2 means cubic splines of order 3
     if(m==-1){
         res <- as.numeric(x<k[i+1] & x>=k[i])
     } else {
         z0 <- (x-k[i]) / (k[i+m+1]-k[i])
         z1 <- (k[i+m+2]-x)/(k[i+m+2]-k[i+1])
         res <- z0*bspline(x,k,i,m-1) + z1*bspline(x,k,i+1,m-1)
     }
     res
}

construct_bspline <- function(x, knots, coefs, m){
     ## Constructs a complete spline from a set of knots and coefficients
     y <- rep(0, times=length(x))
     for(i in seq(coefs)){
         coef <- coefs[i]
         y <- y + coefs[i] * bspline(x, knots, i, m)
     }
     y
}

span_knots <- function(start, end, h, m){
     ## create knots of distance 'h' from 'start' to 'end',
     ##  to use with B-splines of order m+1
     first_knot <- start - start%%h - (m+1)*h
     last_knot <- end-1 + (h - (end-1)%%h) + (m+1)*h

     knots <- seq(first_knot, last_knot, by=h)
     knots
}

## create data:
x <- 125:197
y <- cos(sqrt(40*x)/3) + 3  # mean=3, s(0)=4
plot(x,y)

m <- 2  # use cubic splines of order m+1=3

myknots <- span_knots(start=min(x), end=max(x), h=10, m=m)

## Why do I have to set 'k' to length(myknots)-m-2 here? According to 
the book, it should be length(myknots)-m-1.
my_gam <- gam(y~s(x, bs="ps", m=m, k=length(myknots)-m-2), 
knots=list(x=myknots))

## Here, I construct my own spline with the coefficients of the GAM, 
except coef()[1] which is the intercept
spline_y <- construct_bspline(x, my_gam$smooth[[1]]$knots, 
coef(my_gam)[-1], mm)

plot(my_gam)
lines(x,spline_y)  # why is this different? It seems shifted on the x- 
and y-axis