[R] Bootstrapped CIs of MSE for (G)AM model

[Rui Barradas]

Hello,

There were several errors with your code. The following works but with 
the other CI types.


library(ISLR)
library(mgcv)
library(boot)

# function to obtain MSE
MSE <- function(data, indices, formula) {
   d <- data[indices, ] # allows boot to select sample
   fit <- gam(formula, data = d)
   ypred <- predict(fit)
   mean((d[["wage"]] - ypred)^2)
}

data(Wage)

# Make the results reproducible
set.seed(1234)

# bootstrapping with 1000 replications
results <- boot(data = Wage, statistic = MSE,
                 R = 1000, formula = wage ~ education + s(age, bs = 
"ps") + year))

# get 95% confidence intervals
# type = "bca" is throwing an error
ci.type <- c("norm","basic", "stud", "perc")
boot.ci(results, type = ci.type)



Hope this helps,

Rui Barradas

Às 20:36 de 22/11/2018, varin sacha via R-help escreveu:
> Dear R-experts,
> 
> I am trying to get the bootstrapped confidence intervals of Mean squared error (MSE) for a (G)AM model. I get an error message.
> Here below the reproducible R code. Many thanks for your response.
> 
> ####################
> 
> 
> install.packages("ISLR")
> 
> library(ISLR)
> 
> install.packages("mgcv")
> 
> library(mgcv)
> 
> install.packages("boot")
> 
> library(boot)
> 
> #MSE calculation
> 
> n=dim(Wage)[1]
> 
> p=0.667
> 
> GAM1<-gam(wage ~education+s(age,bs="ps")+year,data=Wage)
> 
> sam=sample(1 :n,floor(p*n),replace=FALSE)
> 
>   
> 
> Training =Wage [sam,]
> 
> Testing = Wage [-sam,]
> 
>   
> 
> ypred=predict(GAM1,newdata=Testing)
> 
> y=Testing$wage
> 
> MSE = mean((y-ypred)^2)
> 
> 
> # Bootstrap 95% CI for MSE
> 
> # function to obtain MSE
> MSE <- function(formula, data, indices) {
>    d <- data[indices,] # allows boot to select sample
>    fit <- gam(formula, data=d)
>    return(MSE(fit))
> 
> } # bootstrapping with 1000 replications
> results <- boot(data=Wage, statistic=MSE,
>     R=1000, formula=gam(wage ~education+s(age,bs="ps")+year,data=Wage)
> 
> )
> # get 95% confidence intervals
> boot.ci(results, type="bca")
> 
> ##########################
> 
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.
>