首页 文章

在tidyverse中按组滚动回归?

提问于
浏览
7

关于在R中滚动回归有很多问题,但在这里我特意寻找使用 dplyr ,_ _1184053_和(如果需要) purrr 的东西 .

这就是使这个问题与众不同的原因 . 我希望与_1844055一致 . 是否可以使用整洁的工具(例如 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 ):

  • 计算 yx 上的滚动回归,超过最后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

EDIT:

这里有一个有趣的后续问题rolling regression with confidence interval (tidyverse)

4 回答

  • 2

    这更像是一个想法而不是答案,但可能不是使用 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)
    
  • 1

    定义一个函数 Coef ,其参数由 cbind(y, x) 形成,并使用截距在x上对y进行回归,返回系数 . 然后使用每组上的当前行和先前行应用 rollapplyr . 如果最后你的意思是当前行的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
    
  • 8

    这会做你想要的吗?

    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
    

    Edit: 稍微修改过的代码,但是 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))))
    

    Edit 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)))
    
  • 2

    这是一个类似于G. Grothendieck's answer但使用 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
    

相关问题