【问题标题】：Confidence interval for quantile regression using bootstrap使用 bootstrap 进行分位数回归的置信区间
【发布时间】：2018-03-05 20:58:54
【问题描述】：

我正在尝试获取线性和分位数回归的五种引导间隔。我能够使用 car 的 Boot 和 boot 的 boot.ci 引导并找到线性回归的 5 个引导间隔（分位数、正常、基本、Studentized 和 BCa） .当我尝试使用 quantreg 中的 rq 对分位数回归执行相同操作时，它会引发错误。这是示例代码

创建模型

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)

引导线性和分位数回归

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) 为什么 boot.ci 不能用于分位数回归

2）使用我从 stackexchange 获得的这个解决方案，我能够找到分位数 CI。

rq 的分位数（百分位数 CI）解

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

我如何为 bootstrap 获取其他 CI（正常、基本、学生化、BCa）。

3) 另外，我用于线性回归的 boot.ci 命令会产生此警告

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

这意味着什么？

【问题讨论】：

标签： r confidence-interval quantile statistics-bootstrap

【解决方案1】：

使用summary.rq 可以计算模型系数的自举标准误差。
有五种 boostrap 方法 (bsmethods) 可用（请参阅 ?boot.rq）。

summary(mod2, se = "boot", bsmethod= "xy")

# Call: rq(formula = income ~ education.c + prestige.c + women.c, data = newdata)
# 
# tau: [1] 0.5
#  
# Coefficients:
#             Value      Std. Error t value    Pr(>|t|)  
# (Intercept) 6542.83599  139.54002   46.88860    0.00000
# education.c  291.57468  117.03314    2.49139    0.01440
# prestige.c    89.68050   22.03406    4.07009    0.00010
# women.c      -48.94856    5.79470   -8.44712    0.00000

要计算引导置信区间，您可以使用以下技巧：

mod1.boot <- Boot(mod1, R=999)
set.seed(1234)
boot.ci(mod1.boot, level = .95, type = "all")

dat2 <- newdata[5:7]
set.seed(1234)
mod2.boot <- boot.rq(cbind(1,dat2),newdata$income,tau=0.5, R=10000)

# Create an object with the same structure of mod1.boot
# but with boostrap replicates given by boot.rq
mod3.boot <- mod1.boot
mod3.boot$R <- 10000
mod3.boot$t0 <- coef(mod2)
mod3.boot$t <- mod2.boot$B
boot.ci(mod3.boot, level = .95, type = "all")

# BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
# Based on 10000 bootstrap replicates
# 
# CALL : 
# boot.ci(boot.out = mod3.boot, type = "all", level = 0.95)
# 
# Intervals : 
# Level      Normal              Basic             Studentized     
# 95%   (6293, 6838 )   (6313, 6827 )   (6289, 6941 )  
# 
# Level     Percentile            BCa          
# 95%   (6258, 6772 )   (6275, 6801 )

【讨论】：

感谢您的解决方案，但我已经想出了一个更通用的解决方案来计算所有变量的 CI。再次感谢！

【解决方案2】：

感谢所有提供帮助的人。我能够自己找出解决方案。我运行了一个循环来计算分位数回归的系数，然后分别使用 boot 和 boot.ci。这是代码

仅引导命令，从问题中创建模型

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

我这样做是因为我想要对所有变量进行 CI，而不仅仅是截距。

【讨论】：