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