2016-06-29 45 views
5

Ich möchte 95% Konfidenzintervalle für die Regressionskoeffizienten einer Quantil-Regression erhalten. Sie können mit der rq Funktion des quantreg Pakets in R (im Vergleich zu einem OLS-Modell) Quantilsregressionen berechnen:Berechnen 95% Konfidenzintervalle in Quantil Regression in R mit RQ-Funktion

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

Ich bin in der Lage 95% Konfidenzintervall für das lineare Modell erhalten mit der confint Funktion:

confint(LM) 

Wenn ich Quantilsregression verwenden ich verstehe, dass der folgende Code Standardfehler Bootstrap produziert:

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

Aber eigentlich hätte ich gerne 95% Konfidenzintervalle. Das heißt, etwas zu interpretieren wie: "Mit einer Wahrscheinlichkeit von 95% enthält das Intervall [...] den wahren Koeffizienten". Wenn ich Standardfehler mit summary.lm() berechne, würde ich einfach SE * 1.96 multiplizieren und ähnliche Ergebnisse wie von confint() erhalten. Dies ist jedoch nicht mit Bootstrapped-Standardfehlern möglich. Also meine Frage ist, wie erhalten 95% Konfidenzintervalle für Quantil-Regressionskoeffizienten?

Antwort

4

können Sie die boot.rq Funktion verwenden, um direkt die Koeffizienten Bootstrap:

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

Sie können das Boot-Paket verwenden Intervalle andere als die Perzentil-Intervall zu berechnen.

+0

Schöne Antwort. Im t (apply ...) 'Fehler in QR.b $ B: $ Operator ist ungültig für atomare Vektoren'. Ich denke keine Spaltenauswahl für das QR.b: 't (apply (QR.b, 2, Quantil, c (0,025,0,975))'? – Minnow

1

Sie können den vcov auch einfach vom Objekt abrufen, indem Sie covariance=TRUE einstellen. Dies entspricht mit boostrapped Standardfehler in Ihrem 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) 

Aber ja, der bessere Weg, dies zu tun, ist ein studentisierten Bootstrap zu verwenden.