5

I'd like to extract coefficients and upper and lower bounds from a quantile regression using the quantreg package. Here's an example from the help file.

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

I know I can use coefficients() to get the 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

But I can't figure out how to get the upper/lower bounds that appear in summary(). I looked at str(sf), but I did not see how to extract.

Ultimately, I'd like to put taus, coefficients, and upper/lower bounds in a dataframe for further processing.

lmo
  • 37,904
  • 9
  • 56
  • 69
Eric Green
  • 7,385
  • 11
  • 56
  • 102

4 Answers4

3

I'm assuming you just want the coefficients on the non-intercept term. How about this

sapply(sf, function(x) c(tau=x$tau, x$coefficients[-1, ]))

That will iterate over the different levels of tau and extract the intervals for the coefficients

                  [,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
  • 195,160
  • 17
  • 277
  • 295
  • an updated version of quantreg doesn't give lower and upper bounds as output in the summary function – NBK Feb 06 '22 at 12:53
1

You can use the function coef with the object returned by summary to extract the values.

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

Here, coef returns a matrix. The lower and upper bounds are in the second and third column, respectively.

Sven Hohenstein
  • 80,497
  • 17
  • 145
  • 168
1

Here's a tidy solution that requires a broom:

    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

See: https://cran.r-project.org/web/packages/broom/vignettes/broom.html

PatrickT
  • 10,037
  • 9
  • 76
  • 111
0

This might help.

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

or

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
Mark Miller
  • 12,483
  • 23
  • 78
  • 132