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]]
