如何从分位数回归rq()中提取系数的上/下限界?

5

我想从使用 quantreg 软件包的分位数回归中提取系数以及上限和下限。以下是来自帮助文件的示例。

data(engel)
attach(engel)
taus <- c(.05,.1,.25,.75,.9,.95)
f <- rq((foodexp)~(income),tau=taus)
sf <- summary(f)
sf[1]
#[[1]]

#Call: rq(formula = (foodexp) ~ (income), tau = taus)

#tau: [1] 0.05

#Coefficients:
#            coefficients lower bd  upper bd 
#(Intercept) 124.88004     98.30212 130.51695
#income        0.34336      0.34333   0.38975

我知道我可以使用 coefficients() 函数来获取系数。
cf <- t(data.frame(coefficients(f)))    # transpose for better arrangement
cf
#              (Intercept)    income
#tau..0.05   124.88004 0.3433611
#tau..0.10   110.14157 0.4017658
#tau..0.25    95.48354 0.4741032
#tau..0.75    62.39659 0.6440141
#tau..0.90    67.35087 0.6862995
#tau..0.95    64.10396 0.7090685

但我不知道如何获取出现在summary()中的上/下限。我查看了str(sf),但我没有看到如何提取。

最终,我希望将tau、系数和上/下限放入数据框以进行进一步处理。

4个回答

3
我假设您只想要非截距项的系数。这样怎么样?
sapply(sf, function(x) c(tau=x$tau, x$coefficients[-1, ]))

那会迭代不同的 tau 级别,并提取系数的间隔。
                  [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
tau          0.0500000 0.1000000 0.2500000 0.7500000 0.9000000 0.9500000
coefficients 0.3433611 0.4017658 0.4741032 0.6440141 0.6862995 0.7090685
lower bd     0.3433270 0.3420992 0.4203298 0.5801552 0.6493680 0.6739000
upper bd     0.3897500 0.4507941 0.4943288 0.6904127 0.7422294 0.7344405

很棒的解决方案,@MrFlick。 - Eric Green
quantreg 的更新版本在 summary 函数中不会输出下限和上限。 - NBK

1
这里有一个需要使用扫帚的整洁解决方案:
    library("broom")
    tidy(rq(data = engel, foodexp ~ income, tau = c(.05,.1,.25,.75,.9,.95)))
              term    estimate   conf.low   conf.high  tau
    1  (Intercept) 124.8800408 96.9260199 131.1526572 0.05
    2       income   0.3433611  0.3256521   0.3957653 0.05
    3  (Intercept) 110.1415742 74.7176192 151.7365540 0.10
    4       income   0.4017658  0.3399444   0.4542858 0.10
    5  (Intercept)  95.4835396 67.5662860 134.9199638 0.25
    6       income   0.4741032  0.4168149   0.5112239 0.25
    7  (Intercept)  62.3965855 27.8317452 120.0708811 0.75
    8       income   0.6440141  0.5738877   0.7051561 0.75
    9  (Intercept)  67.3508721 34.4846868 114.7983532 0.90
    10      income   0.6862995  0.6408633   0.7430180 0.90
    11 (Intercept)  64.1039632 44.4751575 107.6600046 0.95
    12      income   0.7090685  0.6444755   0.7381036 0.95

请参见:https://cran.r-project.org/web/packages/broom/vignettes/broom.html


1
您可以使用summary返回的对象和coef函数来提取值。
library(quantreg)
f <- rq(stack.loss ~ stack.x,.5)

sf <- summary(f)
sf
# Call: rq(formula = stack.loss ~ stack.x, tau = 0.5)

# tau: [1] 0.5

# Coefficients:
#                   coefficients lower bd  upper bd 
# (Intercept)       -39.68986    -41.61973 -29.67754
# stack.xAir.Flow     0.83188      0.51278   1.14117
# stack.xWater.Temp   0.57391      0.32182   1.41090
# stack.xAcid.Conc.  -0.06087     -0.21348  -0.02891

coef(sf)
#                   coefficients    lower bd     upper bd
# (Intercept)       -39.68985507 -41.6197317 -29.67753515
# stack.xAir.Flow     0.83188406   0.5127787   1.14117115
# stack.xWater.Temp   0.57391304   0.3218235   1.41089812
# stack.xAcid.Conc.  -0.06086957  -0.2134829  -0.02891341

在这里,coef 返回一个矩阵。下限和上限分别位于第二和第三列。

谢谢。但是对于我的例子,coef(sf)返回NULL。你也是这样吗? - Eric Green
这对我也没有在OP的例子中起作用。我认为这与多个tau水平有关。 - MrFlick
@MrFlick 没错,我应该使用原帖作者的例子。 - Sven Hohenstein

0

这可能会有所帮助。

matrix(unlist(lapply(sf, function(x) data.frame(unlist(x)[3:8]))), nrow=6, byrow=TRUE)
          [,1]      [,2]     [,3]      [,4]      [,5]      [,6]
[1,] 124.88004 0.3433611 98.30212 0.3433270 130.51695 0.3897500
[2,] 110.14157 0.4017658 79.88753 0.3420992 146.18875 0.4507941
[3,]  95.48354 0.4741032 73.78608 0.4203298 120.09847 0.4943288
[4,]  62.39659 0.6440141 32.74488 0.5801552 107.31362 0.6904127
[5,]  67.35087 0.6862995 37.11802 0.6493680 103.17399 0.7422294
[6,]  64.10396 0.7090685 46.26495 0.6739000  83.57896 0.7344405

或者

matrix(unlist(lapply(sf, function(x) data.frame(unlist(x)[3:8]))), nrow=6, byrow=FALSE)
            [,1]        [,2]        [,3]        [,4]        [,5]       [,6]
[1,] 124.8800408 110.1415742  95.4835396  62.3965855  67.3508721 64.1039632
[2,]   0.3433611   0.4017658   0.4741032   0.6440141   0.6862995  0.7090685
[3,]  98.3021170  79.8875277  73.7860765  32.7448768  37.1180206 46.2649456
[4,]   0.3433270   0.3420992   0.4203298   0.5801552   0.6493680  0.6739000
[5,] 130.5169478 146.1887545 120.0984745 107.3136214 103.1739898 83.5789592
[6,]   0.3897500   0.4507941   0.4943288   0.6904127   0.7422294  0.7344405

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接