在 tidyverse 中按组滚动回归？
rolling regression by group in the tidyverse?
关于 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 的 rolling 回归 .
- 将滚动回归的系数存储在数据框的一列中。
当然，如您所见，只能对每组中的最后 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
编辑：
这里对这个问题进行了有趣的跟进
这与其说是一个答案，不如说是一个想法，但也许可以不使用 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)
这是否符合您的要求？
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)))
定义一个函数 Coef，其参数由 cbind(y, x) 构成，并用截距对 x 进行 y 回归，返回系数。然后使用每个组的当前行和先前行应用 rollapplyr。如果 last 是指当前行的前 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
这是一个类似于 的解决方案，但使用 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
关于 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的 rolling 回归 . - 将滚动回归的系数存储在数据框的一列中。
当然，如您所见，只能对每组中的最后 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
编辑：
这里对这个问题进行了有趣的跟进
这与其说是一个答案，不如说是一个想法，但也许可以不使用 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)
这是否符合您的要求？
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)))
定义一个函数 Coef，其参数由 cbind(y, x) 构成，并用截距对 x 进行 y 回归，返回系数。然后使用每个组的当前行和先前行应用 rollapplyr。如果 last 是指当前行的前 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
这是一个类似于 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