[R] mgcv: I can't manually reconstruct a P-spline from a GAM's coefficients
Alexander Engelhardt
alex at chaotic-neutral.de
Wed Jun 18 19:36:58 CEST 2014
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.
>
>
More information about the R-help
mailing list