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

Simon Wood s.wood at bath.ac.uk
Wed Jun 18 11:03:31 CEST 2014

```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