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