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
jc...@virginia.edu
[[alternative HTML version deleted]]
______________________________________________
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.
