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))
###################################
?
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 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.
>
I forgot, you should also set.seed() before calling boot() to make the results reproducible. Rui Barradas On 5/22/2018 10:00 AM, Rui Barradas wrote:> 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 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.