1

I would like to extract p-values and coefficients from a series of quantile regressions made with a grouping variable. I mostly use dplyr to manipulate data frames so a would like dplyr solution.

require(quantreg)
data("engel")
require(dplyr)
engel$grp <- trunc(runif(nrow(engel), min=0, max=3))

group_by(engel,grp) %>% do(summary(rq(foodexp~income,data=.,tau=c(.05, .25, .5, .75, .95)),se="boot"))

This result in a error

Error: Results are not data frames at positions: 1, 2, 3

I tried another version doing first the models and later the summary

rqm <- group_by(engel,grp) %>% do(mdl=rq(foodexp~income,data=.,tau=c(.05, .25, .5, .75, .95)))
summarise(rqm, coef(summary(mdl,se="boot")))

which also result in error

Error: not a vector

lmo
  • 37,904
  • 9
  • 56
  • 69
Leosar
  • 2,010
  • 4
  • 21
  • 32
  • `do` requires a data.frame as an output. You could try `group_by(engel,grp) %>% do(as.data.frame(lapply(summary(rq(foodexp~income,data=.,tau=c(.05, .25, .5, .75, .95)),se="boot"), coef)))` – jeremycg Aug 26 '15 at 18:13
  • That works! but insted of having all quantiles (taus) in one row how could I put the results for each tau in a separate row? – Leosar Aug 26 '15 at 18:17

1 Answers1

2

It's a mess, but it works:

library(dplyr)
group_by(engel,grp) %>% 
      do(as.data.frame(do.call(rbind,
          lapply(summary(rq(foodexp~income,data=.,tau=c(.05, .25, .5, .75, .95)), se="boot"), coef)
            ), row.names = NULL))

You would then want to label the rows with the tau values and if it is a coef or p value. I'll leave that to you.

jeremycg
  • 24,657
  • 5
  • 63
  • 74
  • I have been using R for a year and I have no idea what is going on here. but this helps me get slightly closer to solving my problem. – Brandon Kuczenski Jan 19 '21 at 10:37
  • 1
    @BrandonKuczenski you might have better luck using the `broom` package these days, it's made for running models across dataframes like this. – jeremycg Jan 19 '21 at 14:14
  • Thanks for the comment! I was actually using `broom::tidy()` and discovered that what I *really* wanted could be accomplished simply by adding the `se = "nid"` switch to my `tidy()` call, cf https://github.com/tidymodels/broom/issues/404 – Brandon Kuczenski Jan 19 '21 at 21:37
  • Really what this did was help me to better understand what the tidyverse tools are doing under the hood – Brandon Kuczenski Jan 19 '21 at 21:40