95% bootstrap CIs

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

95% bootstrap CIs

R help mailing list-2
Dear R-Experts,

Here is my reproducible R code to get the Mean squared error of GAM and MARS after I = 50 iterations/replications.
If I want to get the 95% bootstrap CIs around the MSE of GAM and around the MSE of MARS, how can I complete/modify my R code ?

Many thanks for your precious help.

##################

library(mgcv)
 library(earth)  
my.experiment <- function() {
n<-500
x <-runif(n, 0, 5)
 z <- rnorm(n, 2, 3)
a <- runif(n, 0, 5)
y_model <- 0.1*x^3 - 0.5*z^2 - a + x*z + x*a + 3*x*a*z + 10
 y_obs <- y_model +c( rnorm(n*0.97, 0, 0.1), rnorm(n*0.03, 0, 0.5) )
 gam_model<- gam(y_obs~s(x)+s(z)+s(a))
mars_model<-earth(y_obs~x+z+a)
MSE_GAM<-mean((gam_model$fitted.values - y_model)^2)
 MSE_MARS<-mean((mars_model$fitted.values - y_model)^2)
 return( c(MSE_GAM, MSE_MARS) )
}  
my.data = t(replicate( 50, my.experiment() ))
colnames(my.data) <- c("MSE_GAM", "MSE_MARS")
summary(my.data) 

##################

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: 95% bootstrap CIs

Rui Barradas
Hello,

Is this what you are looking for?


ci95 <- apply(my.data, 2, quantile, probs = c(0.025, 0.975))


Hope this helps,

Rui Barradas

Às 20:42 de 23/09/19, varin sacha via R-help escreveu:

> Dear R-Experts,
>
> Here is my reproducible R code to get the Mean squared error of GAM and MARS after I = 50 iterations/replications.
> If I want to get the 95% bootstrap CIs around the MSE of GAM and around the MSE of MARS, how can I complete/modify my R code ?
>
> Many thanks for your precious help.
>
> ##################
>
> library(mgcv)
>   library(earth)
> my.experiment <- function() {
> n<-500
> x <-runif(n, 0, 5)
>   z <- rnorm(n, 2, 3)
> a <- runif(n, 0, 5)
> y_model <- 0.1*x^3 - 0.5*z^2 - a + x*z + x*a + 3*x*a*z + 10
>   y_obs <- y_model +c( rnorm(n*0.97, 0, 0.1), rnorm(n*0.03, 0, 0.5) )
>   gam_model<- gam(y_obs~s(x)+s(z)+s(a))
> mars_model<-earth(y_obs~x+z+a)
> MSE_GAM<-mean((gam_model$fitted.values - y_model)^2)
>   MSE_MARS<-mean((mars_model$fitted.values - y_model)^2)
>   return( c(MSE_GAM, MSE_MARS) )
> }
> my.data = t(replicate( 50, my.experiment() ))
> colnames(my.data) <- c("MSE_GAM", "MSE_MARS")
> summary(my.data)
>
> ##################
>
> ______________________________________________
> [hidden email] 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.
>

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: 95% bootstrap CIs

R help mailing list-2
In reply to this post by R help mailing list-2
Dear R-experts,

Below the reproducible example. I have tried to write a function that returns the statistic of interest (MSE in my case). I have run boot( ) where the function is included in the statistic argument. I have run boot.ci with the result from boot( ). I guess the error comes from the data : bootResults <- boot(data=?????,statistic=mse, R=1000)
Many thanks for your help.

##################################################
library(mgcv)

library(earth)

library(boot)

 
n<-2000

x <-runif(n, 0, 5)

z <- rnorm(n, 2, 3)

a <- runif(n, 0, 5)


y_model<- 0.1*x^3 - 0.5 * z^2 - a + 10

y_obs<-rnorm(n, y_model, 0.1)

gam_model<- gam(y_obs~s(x)+s(z)+s(a))

mars_model<-earth(y_obs~x+z+a)

 
MSE_GAM<-mean((gam_model$fitted.values - y_model)^2)

MSE_MARS<-mean((mars_model$fitted.values - y_model)^2)

 
MSE_GAM

MSE_MARS

 

mse <- function(data,i) {

boot.gam <- gam(y_obs~s(x)+s(z)+s(a),data=data[i,])

return(mean(boot.gam$residuals^2))

}

bootResults <-boot(data=data,statistic=mse,R=1000)

 

mse <- function(data,i) {

boot.earth <- earth((y_obs~x+z+a),data=data[i,])

return(mean(boot.earth$residuals^2))

}

bootResults <-boot(data=data,statistic=mse,R=1000)
##################################################

 









Le lundi 23 septembre 2019 à 21:42:56 UTC+2, varin sacha via R-help <[hidden email]> a écrit :





Dear R-Experts,

Here is my reproducible R code to get the Mean squared error of GAM and MARS after I = 50 iterations/replications.
If I want to get the 95% bootstrap CIs around the MSE of GAM and around the MSE of MARS, how can I complete/modify my R code ?

Many thanks for your precious help.

##################

library(mgcv)
library(earth)  
my.experiment <- function() {
n<-500
x <-runif(n, 0, 5)
z <- rnorm(n, 2, 3)
a <- runif(n, 0, 5)
y_model <- 0.1*x^3 - 0.5*z^2 - a + x*z + x*a + 3*x*a*z + 10
y_obs <- y_model +c( rnorm(n*0.97, 0, 0.1), rnorm(n*0.03, 0, 0.5) )
gam_model<- gam(y_obs~s(x)+s(z)+s(a))
mars_model<-earth(y_obs~x+z+a)
MSE_GAM<-mean((gam_model$fitted.values - y_model)^2)
MSE_MARS<-mean((mars_model$fitted.values - y_model)^2)
return( c(MSE_GAM, MSE_MARS) )
}  
my.data = t(replicate( 50, my.experiment() ))
colnames(my.data) <- c("MSE_GAM", "MSE_MARS")
summary(my.data) 

##################

______________________________________________
[hidden email] 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.

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: 95% bootstrap CIs

Rui Barradas
Hello,

In your reproducible example you forget to define 'data'.
You should also

set.seed(<some_int_number>)


The following works.


data <- data.frame(a, x, z, y_obs)
boot.ci.type <- c("norm","basic", "perc")

mse_gam <- function(data,i) {
   boot.gam <- gam(y_obs~s(x)+s(z)+s(a),data=data[i,])
   mean(boot.gam$residuals^2)
}

bootResults_gam <-boot(data=data, statistic=mse_gam, R=1000)
boot.ci(bootResults_gam, type = boot.ci.type)


mse <- function(data,i) {
   boot.earth <- earth((y_obs~x+z+a),data=data[i,])
   mean(boot.earth$residuals^2)
}

bootResults <- boot(data=data, statistic=mse, R=1000)
boot.ci(bootResults, type = boot.ci.type)



Hope this helps,

Rui Barradas

Às 13:43 de 25/09/19, varin sacha via R-help escreveu:

> Dear R-experts,
>
> Below the reproducible example. I have tried to write a function that returns the statistic of interest (MSE in my case). I have run boot( ) where the function is included in the statistic argument. I have run boot.ci with the result from boot( ). I guess the error comes from the data : bootResults <- boot(data=?????,statistic=mse, R=1000)
> Many thanks for your help.
>
> ##################################################
> library(mgcv)
>
> library(earth)
>
> library(boot)
>
>  
> n<-2000
>
> x <-runif(n, 0, 5)
>
> z <- rnorm(n, 2, 3)
>
> a <- runif(n, 0, 5)
>
>
> y_model<- 0.1*x^3 - 0.5 * z^2 - a + 10
>
> y_obs<-rnorm(n, y_model, 0.1)
>
> gam_model<- gam(y_obs~s(x)+s(z)+s(a))
>
> mars_model<-earth(y_obs~x+z+a)
>
>  
> MSE_GAM<-mean((gam_model$fitted.values - y_model)^2)
>
> MSE_MARS<-mean((mars_model$fitted.values - y_model)^2)
>
>  
> MSE_GAM
>
> MSE_MARS
>
>  
>
> mse <- function(data,i) {
>
> boot.gam <- gam(y_obs~s(x)+s(z)+s(a),data=data[i,])
>
> return(mean(boot.gam$residuals^2))
>
> }
>
> bootResults <-boot(data=data,statistic=mse,R=1000)
>
>  
>
> mse <- function(data,i) {
>
> boot.earth <- earth((y_obs~x+z+a),data=data[i,])
>
> return(mean(boot.earth$residuals^2))
>
> }
>
> bootResults <-boot(data=data,statistic=mse,R=1000)
> ##################################################
>
>  
>
>
>
>
>
>
>
>
>
> Le lundi 23 septembre 2019 à 21:42:56 UTC+2, varin sacha via R-help <[hidden email]> a écrit :
>
>
>
>
>
> Dear R-Experts,
>
> Here is my reproducible R code to get the Mean squared error of GAM and MARS after I = 50 iterations/replications.
> If I want to get the 95% bootstrap CIs around the MSE of GAM and around the MSE of MARS, how can I complete/modify my R code ?
>
> Many thanks for your precious help.
>
> ##################
>
> library(mgcv)
> library(earth)
> my.experiment <- function() {
> n<-500
> x <-runif(n, 0, 5)
> z <- rnorm(n, 2, 3)
> a <- runif(n, 0, 5)
> y_model <- 0.1*x^3 - 0.5*z^2 - a + x*z + x*a + 3*x*a*z + 10
> y_obs <- y_model +c( rnorm(n*0.97, 0, 0.1), rnorm(n*0.03, 0, 0.5) )
> gam_model<- gam(y_obs~s(x)+s(z)+s(a))
> mars_model<-earth(y_obs~x+z+a)
> MSE_GAM<-mean((gam_model$fitted.values - y_model)^2)
> MSE_MARS<-mean((mars_model$fitted.values - y_model)^2)
> return( c(MSE_GAM, MSE_MARS) )
> }
> my.data = t(replicate( 50, my.experiment() ))
> colnames(my.data) <- c("MSE_GAM", "MSE_MARS")
> summary(my.data)
>
> ##################
>
> ______________________________________________
> [hidden email] 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.
>
> ______________________________________________
> [hidden email] 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.
>

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: 95% bootstrap CIs

R help mailing list-2
Dear Rui,

Excellent ! Many thanks.








Le mercredi 25 septembre 2019 à 18:50:09 UTC+2, Rui Barradas <[hidden email]> a écrit :





Hello,

In your reproducible example you forget to define 'data'.
You should also

set.seed(<some_int_number>)


The following works.


data <- data.frame(a, x, z, y_obs)
boot.ci.type <- c("norm","basic", "perc")

mse_gam <- function(data,i) {
  boot.gam <- gam(y_obs~s(x)+s(z)+s(a),data=data[i,])
  mean(boot.gam$residuals^2)
}

bootResults_gam <-boot(data=data, statistic=mse_gam, R=1000)
boot.ci(bootResults_gam, type = boot.ci.type)


mse <- function(data,i) {
  boot.earth <- earth((y_obs~x+z+a),data=data[i,])
  mean(boot.earth$residuals^2)
}

bootResults <- boot(data=data, statistic=mse, R=1000)
boot.ci(bootResults, type = boot.ci.type)



Hope this helps,

Rui Barradas

Às 13:43 de 25/09/19, varin sacha via R-help escreveu:

> Dear R-experts,
>
> Below the reproducible example. I have tried to write a function that returns the statistic of interest (MSE in my case). I have run boot( ) where the function is included in the statistic argument. I have run boot.ci with the result from boot( ). I guess the error comes from the data : bootResults <- boot(data=?????,statistic=mse, R=1000)
> Many thanks for your help.
>
> ##################################################
> library(mgcv)
>
> library(earth)
>
> library(boot)
>

> n<-2000
>
> x <-runif(n, 0, 5)
>
> z <- rnorm(n, 2, 3)
>
> a <- runif(n, 0, 5)
>
>
> y_model<- 0.1*x^3 - 0.5 * z^2 - a + 10
>
> y_obs<-rnorm(n, y_model, 0.1)
>
> gam_model<- gam(y_obs~s(x)+s(z)+s(a))
>
> mars_model<-earth(y_obs~x+z+a)
>

> MSE_GAM<-mean((gam_model$fitted.values - y_model)^2)
>
> MSE_MARS<-mean((mars_model$fitted.values - y_model)^2)
>

> MSE_GAM
>
> MSE_MARS
>

>
> mse <- function(data,i) {
>
> boot.gam <- gam(y_obs~s(x)+s(z)+s(a),data=data[i,])
>
> return(mean(boot.gam$residuals^2))
>
> }
>
> bootResults <-boot(data=data,statistic=mse,R=1000)
>

>
> mse <- function(data,i) {
>
> boot.earth <- earth((y_obs~x+z+a),data=data[i,])
>
> return(mean(boot.earth$residuals^2))
>
> }
>
> bootResults <-boot(data=data,statistic=mse,R=1000)
> ##################################################
>

>
>
>
>
>
>
>
>
>
> Le lundi 23 septembre 2019 à 21:42:56 UTC+2, varin sacha via R-help <[hidden email]> a écrit :
>
>
>
>
>
> Dear R-Experts,
>
> Here is my reproducible R code to get the Mean squared error of GAM and MARS after I = 50 iterations/replications.
> If I want to get the 95% bootstrap CIs around the MSE of GAM and around the MSE of MARS, how can I complete/modify my R code ?
>
> Many thanks for your precious help.
>
> ##################
>
> library(mgcv)
> library(earth)
> my.experiment <- function() {
> n<-500
> x <-runif(n, 0, 5)
> z <- rnorm(n, 2, 3)
> a <- runif(n, 0, 5)
> y_model <- 0.1*x^3 - 0.5*z^2 - a + x*z + x*a + 3*x*a*z + 10
> y_obs <- y_model +c( rnorm(n*0.97, 0, 0.1), rnorm(n*0.03, 0, 0.5) )
> gam_model<- gam(y_obs~s(x)+s(z)+s(a))
> mars_model<-earth(y_obs~x+z+a)
> MSE_GAM<-mean((gam_model$fitted.values - y_model)^2)
> MSE_MARS<-mean((mars_model$fitted.values - y_model)^2)
> return( c(MSE_GAM, MSE_MARS) )
> }
> my.data = t(replicate( 50, my.experiment() ))
> colnames(my.data) <- c("MSE_GAM", "MSE_MARS")
> summary(my.data)
>
> ##################
>
> ______________________________________________
> [hidden email] 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.
>
> ______________________________________________
> [hidden email] 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.
>

______________________________________________
[hidden email] 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.