R help - Jun 2014 - mgcv: I can't manually reconstruct a P-spline from a GAM's coefficients

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

mgcv:gam automatically imposes sum-to-zero constraints on the smooths in 
a model (even if there is only one smooth). This is to avoid lack of 
identifiability with the intercept. The constraint removes one 
coefficient and shifts the curve...

best,
Simon

On 17/06/14 15:40, Alexander Engelhardt wrote:
> 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
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
Simon Wood, Mathematical Science, University of Bath BA2 7AY UK
+44 (0)1225 386603               http://people.bath.ac.uk/sw283