[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