Hello,

We have a question about how to retrieve the nonparametric curve in gam() function. Now, we get the estimates and draw the fitted curve using the following code.  You can see that our fitted curve is parallel to the output from gam() function but differs by a constant shift. However, the magnitude of the shift is not the intercept term. Also, gam's estimate appears more close to the (true value - intercept). Can anyone tell us how to get the correct estimate?

x <- runif(200,0,1)
f<-x^2+1
y <- rnorm(length(f), f, 1)
z<-matrix(0,nrow=200,ncol=2)
z[,1]<-y
z[,2]<-x
tt<-data.frame(z)
## finished simulating data, now fit model...
a<-gam(X1~s(X2,bs="tr",m=2),data=tt,family=gaussian)
bb<-s(X2,bs="tr",m=2)
basi<-smooth.construct.tr.smooth.spec(bb,tt,knots=list(a$smooth[[1]]$knots))

# plot, true value in red, our retrieved curve in blue, gam() plot in black.

# gam()'s estimate
plot(a, ylim=c(-.7, .7))

# True value - intercept
x.2<-sort(x)
lines(x.2, x.2^2+1-coef(a)[1], lty=1, col="Red")

# Our estimate; without intercept
f<-basi$X[,-1]%*%coef(a)[-1]
f.1=f[order(x)]
lines(x.2, f.1, lty=1, col="Blue")

Similarly, in the case of two nonparametric functions, e.g., y=beta0 + f(x1) + g(x2) + e, how can we retrieve the estimates of f(x1) and g(x2)?

We also want to calculate the confidence interval. We use the code:

var.y<-(basi$X)%*%a$Vp%*%t(basi$X)
diag.y=diag(var.y)[order(x)]

plot(a, ylim=c(-1, 1))

lines(x.2, x.2^2+1-coef(a)[1], lty=1, col="Red")

lines(x.2, f.1, lty=1, col="Blue", ylim=c(-1, 1))
lines(x.2, f.1-1.96*sqrt(diag.y), lty=4, col="Blue")
lines(x.2, f.1+1.96*sqrt(diag.y), lty=4, col="Blue")

However, the CI obtained this way is different from gam's output. Can anyone show us how to get the correct CI? Especially when we have a GLM, e.g., Poission model?

Thanks a lot for your help.

Jinsong Chen
Division of Biostatistics and Epidemiology
Dept. of Public Health Sciences
School of Medicine
University of Virginia

jc8ww@virginia.edu


	[[alternative HTML version deleted]]

