2016-06-29 14 views
5

Я бы хотел получить 95% доверительные интервалы для коэффициентов регрессии регрессии квантилей. Вы можете рассчитать квантилей регрессии с помощью rq функции quantreg пакета в R (по сравнению с моделью МНК):Расчет 95% доверительных интервалов в квантильной регрессии в R с использованием функции rq

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

Я могу получить 95% доверительные интервалы для линейной модели с использованием функции confint:

confint(LM) 

Когда я использую квантиль регрессию Я понимаю, что следующий код производит бутстрапированный стандартные ошибки:

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

Но на самом деле я хотел бы что-то вроде 95% доверительных интервалов. То есть что-то интерпретировать как: «с вероятностью 95% интервал [...] включает истинный коэффициент». Когда я вычисляю стандартные ошибки с помощью summary.lm(), я бы просто умножал SE * 1.96 и получал аналогичные результаты, как из confint(). Но это невозможно, используя загрузочные стандартные ошибки. Итак, мой вопрос в том, как получить доверительные интервалы 95% для коэффициентов регрессии квантилей?

ответ

4

Вы можете использовать функцию boot.rq непосредственно самонастройки коэффициентов:

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') 
} 

Вы можете использовать пакет загрузки для вычисления отличных от процентиля интервала интервалов.

+0

Хороший ответ. В t (применить ...) 'Ошибка в QR.b $ B: $ оператор недействителен для атомных векторов'. Я думаю, что никакой выбор столбцов для QR.b: 't (apply (QR.b, 2, quantile, c (0.025,0.975)))? – Minnow

1

Вы также можете просто извлечь vcov из объекта, установив covariance=TRUE. Это равносильно использованию boostrapped стандартных ошибок в 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) 

Но да, лучший способ сделать это состоит в использовании стьюдентизированной начальной загрузки.