[R] Bootstraping GAMs for confidence intervales calculation
Luis Reino
luisreino at isa.utl.pt
Tue Jul 29 18:11:19 CEST 2008
Dear R-Users,
I am resending this message just to reminder my question regarding the
calculation of a bootstrap confidence intervals for a GAM plot.
I am trying to apply a bootstrap to a GAM in order to calculate the 95%
confidence intervals for a smooth curve obtained by the “plot.gam”
function of the mgcv package. Nonetheless, I am getting some
difficulties in transposing the results for the graphs.
I used the following commands in R, “mgcv” and “boot” packages:
*> attach(bbvc_11Jul08)*
*> model.boot<-function(data,indices){*
*+ sub.data<-data[indices,]*
*+
model<-gam(asin(sqrt(Cov0_30))~s(age,k=4,bs="cr"),family=gaussian,data=sub.data)*
*+ coef(model)}*
*> gam.boot<-boot(bbvc_11Jul08,model.boot,R=1000)*
* *
*> gam.boot*
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = bbvc_11Jul08, statistic = model.boot, R = 1000)
Bootstrap Statistics :
original bias std. error
t1* 0.40370112 0.0002762674 0.02017378
t2* 0.04715501 -0.0144132280 0.07648348
t3* 0.14590936 -0.0135820501 0.05161244
t4* 0.13520098 -0.0096405733 0.04470793
I obtained 4 indexes in the “bootstrap” output. I know that for a lm
function, with an independent variable, t1 is the index of the intercept
and t2 is the index of the variable, but for GAMs I don't find what each
indexes really mean!
Moreover, to calculate the confidence intervals of 95% I proceed as follows:
* *
*>boot.ci(gam.boot,conf=0.95) *
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 1000 bootstrap replicates
CALL :
boot.ci(boot.out = gam.boot, conf = 0.95)
Intervals :
Level Normal Basic Studentized
95% ( 0.3639, 0.4430 ) ( 0.3635, 0.4434 ) ( 0.3583, 0.4507 )
Level Percentile BCa
95% ( 0.3640, 0.4439 ) ( 0.3623, 0.4434 )
Calculations and Intervals on Original Scale
Warning message:
In sqrt(tv[, 2]) : NaNs produced
My first doubt is when I make "boot.ci", which of the four should I
choose? By default, R calculates "index=1:min(2,length(boot.out$t0)" but
I don't know if this is the correct form to do it.
Finally, I don’t know how to introduce the 95% confidence intervals in
the "plot.gam".
I hope you can help me. Thanks in advance.
Luís Reino
--
Luis Reino, PhD Student
Centro de Estudos Florestais
Departamento de Engenharia Florestal
Instituto Superior de Agronomia
Universidade Técnica de Lisboa
1349-017 Lisboa
Portugal
More information about the R-help
mailing list