Search code examples
rconfidence-intervalquantilestatistics-bootstrap

Confidence interval for quantile regression using bootstrap


I am trying to get the five types of bootstrap intervals for linear and quantile regression. I was able to bootstrap and find the 5 boostrap intervals (Quantile,Normal,Basic,Studentized and BCa) for the linear regression using Boot from car and boot.ci from boot. When i tried to do the same for quantile regression using rq from quantreg, it throws up an error. Here is the sample code

Creating the model

library(car)
library(quantreg)
library(boot)
newdata = Prestige[,c(1:4)]
education.c = scale(newdata$education, center=TRUE, scale=FALSE)
prestige.c = scale(newdata$prestige, center=TRUE, scale=FALSE)
women.c = scale(newdata$women, center=TRUE, scale=FALSE)
new.c.vars = cbind(education.c, prestige.c, women.c)
newdata = cbind(newdata, new.c.vars)
names(newdata)[5:7] = c("education.c", "prestige.c", "women.c" )
mod1 = lm(income ~ education.c + prestige.c + women.c, data=newdata)
mod2 = rq(income ~ education.c + prestige.c + women.c, data=newdata)

Booting linear and quantile regression

mod1.boot <- Boot(mod1, R=999)
boot.ci(mod1.boot, level = .95, type = "all")
dat2 <- newdata[5:7]
mod2.boot <- boot.rq(cbind(1,dat2),newdata$income,tau=0.5, R=10000)
boot.ci(mod2.boot, level = .95, type = "all")
Error in if (ncol(boot.out$t) < max(index)) { : 
argument is of length zero

1) Why does boot.ci not work for quantile regression

2)Using this solution I got from stackexchange, I was able to find the quantile CI.

Solution for quantile(percentile CI) for rq

t(apply(mod2.boot$B, 2, quantile, c(0.025,0.975)))

how do i obtain other CI for bootstrap (normal, basic, studentized, BCa).

3) Also, my boot.ci command for linear regression produces this warning

Warning message:
In sqrt(tv[, 2L]) : NaNs produced

What does this signify?


Solution

  • Thanks for everyone who helped. I was able to figure out the solution myself. I ran a loop calculating the coefficients of the quantile regression and then used boot and boot.ci respectively. Here is the code

    Booting commands only, model creation from question

    mod3 <- formula(income ~ education.c + prestige.c + women.c)
    coefsf <- function(data,ind){
    rq(mod3, data=newdata[ind,])$coef
    }
    boot.mod <- boot(newdata,coefsf,R=10000)
    myboot.ci <- list()
    for (i in 1:ncol(boot.mod$t)){
    myboot.ci[[i]] <- boot.ci(boot.mod, level = .95, type = 
    c("norm","basic","perc", "bca"),index = i)
      }
    

    I did this as I wanted CI on all variables not just the intercept.