在tidyverse中按组进行滚动回归?

10

关于在R中进行滚动回归的问题有很多,但我这里特别要求使用dplyrbroom和(如有必要)purrr来完成。这就是让这个问题与众不同的地方。我想要保持tidyverse的一致性。是否可能使用诸如purrr:mapdplyr等整洁工具进行恰当的滚动回归呢?

请考虑这个简单的例子:

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

对于每个组(在本例中,ab):

  1. 计算y在过去2个观测期内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(在这个例子中,yx之间确实存在完美的线性相关性)

更具体地说:

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)滚动回归置信区间


可能是使用rollapply和lm处理多列数据的重复问题。 - Maurits Evers
我知道那个SO问题。我正在寻找的是一个purrr、tidyverse解决方案,这在你提供的链接中并不可用。 - ℕʘʘḆḽḘ
4个回答

13

定义一个函数Coef,其参数由cbind(y, x)形成,该函数拟合带有截距的y和x的回归,并返回系数。然后使用当前行和前一行对每个组进行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

真的很酷,但我对这里的超小系数感到困惑?此外,是否可能有适当的列命名? - ℕʘʘḆḽḘ
1
小系数基本上为0,是由于浮点近似所导致的,与rollapply或dplyr无关。我已经更新了代码并输出了几次,现在名称应该没问题了。 - G. Grothendieck
谢谢Grothendieck。我现在意识到这里有一些重要的缺失(我的错)。目前的代码不能让我将回归结果合并回原始数据集,因为行索引信息丢失了(在你的代码中只选择了(y,x))。我已经更新了我的问题中的示例数据框。你能否看到一种简单的方法来更新你的好答案,以便输出数据框中同时出现组和日期?这样我们就可以在需要时进行合并。谢谢! - ℕʘʘḆḽḘ
我尝试过这个,但我想知道数据是否总是正确地对齐在这里 mydata %>% group_by(group) %>% do(cbind(reg_col = select(., y, x) %>% rollapply(2, Coef, by.column = FALSE, align = 'right', fill = NA), date_col = select(., date))) %>% ungroup - ℕʘʘḆḽḘ
1
是的,如果你使用fill=NA,数据将会对齐。这将用NA值填充原本缺失的行。(还建议使用rollapplyr而不是align = "right"来使代码更简洁。) - G. Grothendieck
疯狂的事情,我想知道在处理大数据集时哪个更快、更节省内存。 - ℕʘʘḆḽḘ

2
这更像是一个想法而不是一个答案,但也许可以尝试使用 map 和你的组列表,而不是使用 group_by
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)

那很有趣。你能否尝试在这个上下文中使它工作? - ℕʘʘḆḽḘ
@ℕʘʘḆḽḘ 我添加了一些额外的代码行以返回更干净的结果。我对zoo不太熟悉,结果看起来很奇怪。斜率不应该是0.5和0.33吗? - johnson-shuffle
请查看更新后的问题。谢谢! - ℕʘʘḆḽḘ

2
这是否符合您的需求?
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)))

我只有几个问题。这里的NA是什么目的?为什么要使用margin?你不能在代码中使用dplyr::data_framepurrr:map代替data.frameapply吗? - ℕʘʘḆḽḘ
@ℕʘʘḆḽḘ - NA 的存在是因为第一个观测没有“前2个”观测。我原以为您想要在第1行和第2行、第2行和第3行、第3行和第4行上进行回归。但是,由于未对第0行和第1行执行回归,因此第一行没有这样的回归。我无法弄清楚如何让 data_frame 正常工作,但实际上 apply 是不必要的(我忘记从第一次尝试中删除它了)。请参见编辑。 - Luke C
啊,太糟糕了...我认为手动添加这个NA不安全。我希望这可以自动化。看一下 Grothendieck 的回答,也许?无论如何,谢谢! - ℕʘʘḆḽḘ
@ℕʘʘḆḽḘ 我非常喜欢Grothendieck的回答,它很简洁,我认为提供了更容易灵活的输出!就我所理解的而言,在任何从第一行开始滚动的情况下,使用c(NA, ...)在功能上与使用fill = NA是相同的-如上面的编辑2所示。 - Luke C

2
这是一个类似于G. Grothendieck的答案的解决方案,但使用了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

你想到如何并行化这个程序了吗?我已经让它在一个有68k条记录的数据框上运行成功,但是要花费一个小时的时间! - philiporlando
你可以使用parallel包,特别是clusterCallclusterEvalQ,就像这里一样,将数据集分割成每个进程上唯一的groups集合。 - Benjamin Christoffersen
@philiporlando 你可以尝试在 data.table 中使用新的 frollapply 函数,它支持并行计算。请注意,这个功能目前还没有合并到主分支,你需要安装 PR 分支来使用。 - undefined

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接