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@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@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@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@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.
