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@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
______________________________________________
R-help@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.
