Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Calculating 95% confidence intervals in quantile regression in R using rq function

I would like to get 95% confidence intervals for the regression coefficients of a quantile regression. You can calculate quantile regressions using the rq function of the quantreg package in R (compared to an OLS model):

library(quantreg)
LM<-lm(mpg~disp, data = mtcars)
QR<-rq(mpg~disp, data = mtcars, tau=0.5)

I am able to get 95% confidence intervals for the linear model using the confint function:

confint(LM)

When I use quantile regression I understand that the following code produces bootstrapped standard errors:

summary.rq(QR,se="boot")

But actually I would like something like 95% confidence intervals. That is, something to interprete like: "with a probability of 95%, the interval [...] includes the true coefficient". When I calculate standard errors using summary.lm() I would just multiply SE*1.96 and get similar results as from confint(). But this is not possible using bootstrapped standard errors. So my question is how get 95% confidence intervals for quantile regression coefficients?

like image 881
ehi Avatar asked Jun 29 '16 17:06

ehi


People also ask

What is the quantile for 95 confidence interval?

We round j and k up to the next integer. Then the 95% confidence interval is between the jth and the kth observations in the ordered data. The 95% confidence interval is thus from the 22nd to the 36th observation, 3.75 to 4.30 litres from the Table.

What is tau in quantile regression?

tau. the quantile(s) to be estimated, this is generally a number strictly between 0 and 1, but if specified strictly outside this range, it is presumed that the solutions for all values of tau in (0,1) are desired. In the former case an object of class "rq" is returned, in the latter, an object of class "rq.


2 Answers

You can use the boot.rq function directly to bootstrap the coefficients:

x<-1:50
y<-c(x[1:48]+rnorm(48,0,5),rnorm(2,150,5))

QR <- rq(y~x, tau=0.5)
summary(QR, se='boot')

LM<-lm(y~x)

QR.b <- boot.rq(cbind(1,x),y,tau=0.5, R=10000)

t(apply(QR.b$B, 2, quantile, c(0.025,0.975)))
confint(LM)


plot(x,y)
abline(coefficients(LM),col="green")
abline(coefficients(QR),col="blue")

for(i in seq_len(nrow(QR.b$B))) {
  abline(QR.b$B[i,1], QR.b$B[i,2], col='#0000ff01')
}

You may want to use the boot package to compute intervals other than the percentile interval.

like image 52
Greg Snow Avatar answered Nov 03 '22 02:11

Greg Snow


You could also simply retrieve the vcov from the object, setting covariance=TRUE. This amounts to using boostrapped standard errors in your CI:

vcov.rq <- function(x, se = "iid") {
 vc <- summary.rq(x, se=se, cov=TRUE)$cov
 dimnames(vc) <- list(names(coef(x)), names(coef(x)))
 vc
}

confint(QR)

But yes, the better way to do this is to use a studentized bootstrap.

like image 30
Matifou Avatar answered Nov 03 '22 01:11

Matifou