Извлечь коэффициент с p из квантильной регрессии по группе с dplyr

Я хотел бы извлечь p-значения и коэффициенты из серии квантильных регрессий, сделанных с помощью группирующей переменной. В основном я использую dplyr для управления фреймами данных, поэтому мне нужно решение dplyr.

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

Это приводит к ошибке

Ошибка: результаты не являются фреймами данных в позициях: 1, 2, 3

Я попробовал другую версию, сделав сначала модели, а затем сводку.

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

что также приводит к ошибке

Ошибка: не вектор


r quantile dplyr quantreg
person Leosar    schedule 26.08.2015    source источник
comment
do требует data.frame в качестве вывода. Вы можете попробовать group_by(engel,grp) %>% do(as.data.frame(lapply(summary(rq(foodexp~income,data=.,tau=c(.05, .25, .5, .75, .95)),se="boot"), coef)))   -  person jeremycg    schedule 26.08.2015
comment
Это работает! но вместо того, чтобы иметь все квантили (тау) в одной строке, как я могу поместить результаты для каждого тау в отдельную строку?   -  person Leosar    schedule 26.08.2015

Ответы (1)


arrow_upward
2
arrow_downward

Это беспорядок, но он работает:

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

Затем вы захотите пометить строки значениями тау и если это значение coef или p. Я оставлю это вам.

person jeremycg    schedule 26.08.2015
comment
Я использую R в течение года, и я понятия не имею, что здесь происходит. но это помогает мне немного приблизиться к решению моей проблемы. - person Brandon Kuczenski; 19.01.2021
comment
@BrandonKuczenski, возможно, вам повезет больше, используя пакет broom в наши дни, он создан для запуска моделей в таких фреймах данных, как этот. - person jeremycg; 19.01.2021
comment
Спасибо за комментарий! На самом деле я использовал broom::tidy() и обнаружил, что то, что я на самом деле хотел, можно выполнить, просто добавив переключатель se = "nid" к моему вызову tidy(), см. github.com/tidymodels/broom/issues/404 - person Brandon Kuczenski; 20.01.2021
comment
На самом деле это помогло мне лучше понять, что инструменты tidyverse делают под капотом. - person Brandon Kuczenski; 20.01.2021