[R] Bootstrap and average median squared error

[Rui Barradas]

Hello,

If you want to bootstrap a statistic, I suggest you use base package boot.
You would need the data in a data.frame, see how you could do it.


library(boot)

bootMedianSE <- function(data, indices){
     d <- data[indices, ]
     fit <- rq(crp ~ bmi + glucose, tau = 0.5, data = d)
     ypred <- predict(fit)
     y <- d$crp
     median(y - ypred)^2
}

dat <- data.frame(crp, bmi, glucose)
nboot <- 100

medse <- boot(dat, bootMedianSE, R = nboot)

medse$t0
mean(medse$t)    # This is the value you want


Hope this helps,

Rui Barradas



On 5/22/2018 12:19 AM, varin sacha via R-help wrote:
> Dear R-experts,
> 
> I am trying to bootstrap (and average) the median squared error evaluation metric for a robust regression. I can't get it. What is going wrong ?
> 
> Here is the reproducible example.
> 
> #############################
> install.packages( "quantreg" )
> library(quantreg)
> 
> crp <-c(12,14,13,24,25,34,45,56,25,34,47,44,35,24,53,44,55,46,36,67)
> bmi <-c(34,32,12,76,54,34,21,18,92,32,11,13,45,46,56,57,67,87,12,13)
> glucose <-c(23,54,11,12,13,21,32,12,45,54,65,87,21,23,12,12,23,23,43,54)
> 
> # Create a list to store the results
> lst<-list()
> 
> # Numbers of bootstrap samples
> nboot=100
> bootstrap.MedAESQ =rep(NA,nboot)
> 
> for(i in 1 :nboot)
> {
> 
> fit <- rq( crp ~ bmi+glucose, tau = 0.5)
> ypred=predict(fit)
> y=new$crp
> bootstrap.MedAESQ [i]=median(y-ypred)^2
> lst[i]<-bootstrap.MedAESQ
> 
> }
> mean(unlist(lst))
> ###################################
> 
>   
> 
> ______________________________________________
> 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.
>