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@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland 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-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland 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.
