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?
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.
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.
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.
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.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With