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

[Alexander Engelhardt]

Hello Simon,

thanks for your response!

I incorporated the intercept, i.e. coef(my_gam)[1] into the spline. That 
is, my B-spline coefficients are the intercept for the first one, and 
for all others I added the intercept to them (see the code below if I'm 
not clear)

Now I *almost* reproduce the spline:

http://i.imgur.com/pDZuHev.png

The red line is what plot(gam) does, the black line is my manual 
reproduction.
Could you hint me towards where my mistake is?

Best wishes,
  Alex

New code:

my_gam <- gam(y~s(x, bs="ps", m=m, k=length(myknots)-m-2), 
knots=list(x=myknots))

intercept <- coef(my_gam)[1]
coefs <- rep(NA, length(coef(my_gam)))
coefs[1]  <- intercept
coefs[-1] <- intercept + coef(my_gam)[-1]

spline_y <- construct_bspline(x, my_gam$smooth[[1]]$knots, coefs, m)

plot(my_gam, lwd=2, col="red")
lines(x,spline_y - mean(spline_y))  # subtract mean(spline_y) to have it 
centered around 0


On 06/18/2014 11:03 AM, Simon Wood wrote:
> 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.
>
>