Скользящая регрессия по группам в тидиверсе?
Есть много вопросов о скользящей регрессии в R, но здесь я специально ищу что-то, что использует dplyr, broom и (при необходимости) purrr,
Это то, что делает этот вопрос другим. я хочу быть tidyverse последовательны. Можно ли сделать правильную регрессию бега с помощью таких инструментов, как purrr:map а также dplyr?
Пожалуйста, рассмотрите этот простой пример:
library(dplyr)
library(purrr)
library(broom)
library(zoo)
library(lubridate)
mydata = data_frame('group' = c('a','a', 'a','a','b', 'b', 'b', 'b'),
'y' = c(1,2,3,4,2,3,4,5),
'x' = c(2,4,6,8,6,9,12,15),
'date' = c(ymd('2016-06-01', '2016-06-02', '2016-06-03', '2016-06-04',
'2016-06-03', '2016-06-04', '2016-06-05','2016-06-06')))
group y x date
<chr> <dbl> <dbl> <date>
1 a 1.00 2.00 2016-06-01
2 a 2.00 4.00 2016-06-02
3 a 3.00 6.00 2016-06-03
4 a 4.00 8.00 2016-06-04
5 b 2.00 6.00 2016-06-03
6 b 3.00 9.00 2016-06-04
7 b 4.00 12.0 2016-06-05
8 b 5.00 15.0 2016-06-06
Для каждой группы (в этом примере a или же b):
- вычислить скользящую регрессию
yнаxза последние 2 наблюдения. - сохранить коэффициент этой скользящей регрессии в столбце данных.
Конечно, как вы можете видеть, скользящая регрессия может быть вычислена только для последних 2 строк в каждой группе.
Я пытался использовать следующее, но безуспешно.
data %>% group_by(group) %>%
mutate(rolling_coef = do(tidy(rollapply(. ,
width=2,
FUN = function(df) {t = lm(formula=y ~ x,
data = as.data.frame(df),
na.rm=TRUE);
return(t$coef) },
by.column=FALSE, align="right"))))
Error in mutate_impl(.data, dots) :
Evaluation error: subscript out of bounds.
In addition: There were 21 warnings (use warnings() to see them)
Есть идеи?
Ожидаемый результат для двух последних строк первого a группа составляет 0,5 и 0,5 (действительно существует идеальная линейная корреляция между y а также x в этом примере)
Более конкретно:
mydata_1 <- mydata %>% filter(group == 'a',
row_number() %in% c(1,2))
# A tibble: 2 x 3
group y x
<chr> <dbl> <dbl>
1 a 1.00 2.00
2 a 2.00 4.00
> tidy(lm(y ~ x, mydata_1))['estimate'][2,]
[1] 0.5
а также
mydata_2 <- mydata %>% filter(group == 'a',
row_number() %in% c(2,3))
# A tibble: 2 x 3
group y x
<chr> <dbl> <dbl>
1 a 2.00 4.00
2 a 3.00 6.00
> tidy(lm(y ~ x, mydata_2))['estimate'][2,]
[1] 0.5
РЕДАКТИРОВАТЬ:
интересное продолжение этого вопроса здесь скользящая регрессия с доверительным интервалом (tidyverse)
4 ответа
Определить функцию Coef чей аргумент сформирован из cbind(y, x) и который регрессирует y на x с перехватом, возвращая коэффициенты. Тогда подать заявку rollapplyr используя текущие и предыдущие строки по каждой группе. Если под последним вы подразумевали 2 предыдущие строки для текущей строки, т.е. исключили текущую строку, то замените 2 на list(-seq(2)) в качестве аргумента rollapplyr,
Coef <- . %>% as.data.frame %>% lm %>% coef
mydata %>%
group_by(group) %>%
do(cbind(reg_col = select(., y, x) %>% rollapplyr(2, Coef, by.column = FALSE, fill = NA),
date_col = select(., date))) %>%
ungroup
давая:
# A tibble: 8 x 4
group `reg_col.(Intercept)` reg_col.x date
<chr> <dbl> <dbl> <date>
1 a NA NA 2016-06-01
2 a 0 0.500 2016-06-02
3 a 0 0.500 2016-06-03
4 a 0 0.500 2016-06-04
5 b NA NA 2016-06-03
6 b 0.00000000000000126 0.333 2016-06-04
7 b - 0.00000000000000251 0.333 2016-06-05
8 b 0 0.333 2016-06-06
варьирование
Вариант вышеупомянутого будет:
mydata %>%
group_by(group) %>%
do(select(., date, y, x) %>%
read.zoo %>%
rollapplyr(2, Coef, by.column = FALSE, fill = NA) %>%
fortify.zoo(names = "date")
) %>%
ungroup
Только уклон
Если необходим только уклон, возможны дальнейшие упрощения. Мы используем тот факт, что наклон равен cov(x, y) / var(x),
slope <- . %>% { cov(.[, 2], .[, 1]) / var(.[, 2])}
mydata %>%
group_by(group) %>%
mutate(slope = rollapplyr(cbind(y, x), 2, slope, by.column = FALSE, fill = NA)) %>%
ungroup
Это делает то, что вы после?
data %>%
group_by(group) %>%
do(data.frame(., rolling_coef = c(NA, rollapply(data = ., width = 2, FUN = function(df_) {
d = data.frame(df_)
d[, 2:3] <- apply(d[,2:3], MARGIN = 2, FUN = as.numeric)
mod = lm(y ~ x, data = d)
return(coef(mod)[2])
}, by.column = FALSE, align = "right"))))
Предоставление:
# A tibble: 8 x 4
# Groups: group [2]
group y x rolling_coef
<chr> <dbl> <dbl> <dbl>
1 a 1. 2. NA
2 a 2. 4. 0.500
3 a 3. 6. 0.500
4 a 4. 8. 0.500
5 b 2. 6. NA
6 b 3. 9. 0.333
7 b 4. 12. 0.333
8 b 5. 15. 0.333
Изменить: Слегка измененный код, но data_frame не примет . групповой заполнитель в качестве аргумента - не уверен, как это исправить.
data %>%
group_by(group) %>%
do(data.frame(., rolling_coef = c(NA, rollapplyr(data = ., width = 2, FUN = function(df_) {
mod = lm(y ~ x, data = .)
return(coef(mod)[2])
}, by.column = FALSE))))
Изменить 2: Использование fill = NA вместо того, чтобы использовать c(NA, ...) достигает того же результата.
data %>%
group_by(group) %>%
do(data.frame(., rolling_coef = rollapplyr(data = ., width = 2, FUN = function(df_) {
mod = lm(y ~ x, data = .)
return(coef(mod)[2])
}, by.column = FALSE, fill = NA)))
Это скорее идея, чем ответ, но, возможно, вместо использования group_by попробуйте использовать map и ваш список групп:
FUN <- function(g, df = NULL) {
tmp <- tidy(rollapply(
zoo(filter(df, group == g)),
width = 2,
FUN = function(z) {
t <- lm(y ~ x, data = as.data.frame(z)) ; return(t$coef)
},
by.column = FALSE,
align = "right"
))
tmp$series <- c(rep('intercept', nrow(tmp) / 2), rep('slope', nrow(tmp) / 2))
spread(tmp, series, value) %>% mutate(group = g)
}
map_dfr(list('a', 'b'), FUN, df = data)
Вот решение, аналогичное ответу Г. Гротендика, но с использованием rollRegres пакет. Я должен увеличить width аргумент 3, чтобы избежать ошибки (кстати, зачем вам регрессия с таким небольшим количеством наблюдений?)
library(rollRegres)
Coef <- . %>% { roll_regres.fit(x = cbind(1, .$x), y = .$y, width = 2L)$coefs }
mydata %>%
group_by(group) %>%
do(cbind(reg_col = select(., y, x) %>% Coef,
date_col = select(., date))) %>%
ungroup
#R Error in mydata %>% group_by(group) %>% do(cbind(reg_col = select(., y, :
#R Assertion on 'width' failed: All elements must be >= 3.
# change width to avoid error
Coef <- . %>% { roll_regres.fit(x = cbind(1, .$x), y = .$y, width = 3L)$coefs }
mydata %>%
group_by(group) %>%
do(cbind(reg_col = select(., y, x) %>% Coef,
date_col = select(., date))) %>%
ungroup
#R # A tibble: 8 x 4
#R group reg_col.1 reg_col.2 date
#R <chr> <dbl> <dbl> <date>
#R 1 a NA NA 2016-06-01
#R 2 a NA NA 2016-06-02
#R 3 a 1.54e-15 0.500 2016-06-03
#R 4 a -5.13e-15 0.5 2016-06-04
#R 5 b NA NA 2016-06-03
#R 6 b NA NA 2016-06-04
#R 7 b -3.08e-15 0.333 2016-06-05
#R 8 b -4.62e-15 0.333 2016-06-06
#R Warning messages:
#R 1: In evalq((function (..., call. = TRUE, immediate. = FALSE, noBreaks. = FALSE, :
#R low sample size relative to number of parameters
#R 2: In evalq((function (..., call. = TRUE, immediate. = FALSE, noBreaks. = FALSE, :
#R low sample size relative to number of parameters