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)? ##################
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) > > ################## > > ______________________________________________ > R-help at 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. >
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 <r-help at r-project.org> 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)? ################## ______________________________________________ R-help at 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.
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 <r-help at r-project.org> 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) > > ################## > > ______________________________________________ > R-help at 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. > > ______________________________________________ > R-help at 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. >