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

[Simon Wood]

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