代码之家  ›  专栏  ›  技术社区  ›  ℕʘʘḆḽḘ

tidyverse中按组滚动回归?

rolling-computation broom purrr dplyr r
  •  10
  • ℕʘʘḆḽḘ  · 技术社区  · 7 年前

    关于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 ):

    1. 计算 滚动 的回归 y 在…上 x 在过去2次观察中 。
    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 和 十、 在本例中)

    更具体地说: 

    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

    编辑:

    有趣的后续问题 rolling regression with confidence interval (tidyverse)

    4 回复  |  直到 6 年前
        1
  •  13
  •   G. Grothendieck    7 年前

    定义函数 Coef 其论点由 cbind(y, x) 它用截距对x上的y进行回归,返回系数。然后应用 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
        2
  •  2
  •   johnson-shuffle    7 年前

    这与其说是一个答案,不如说是一个想法,但也许不是使用 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)
        3
  •  2
  •   Luke C    7 年前

    这是你想要的吗? 

    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)))
        4
  •  2
  •   Benjamin Christoffersen    6 年前

    这里有一个类似的解决方案 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
    推荐文章
    MCP_infiltrator  ·  在“tible”中添加一列,给出其列表位置
    2 年前
    sbac  ·  如何对R中数据帧的几列应用统计测试
    3 年前
    John-Henry  ·  使用“purrr::map”将“lm”对象循环到“broom::tidy”`
    3 年前
    deschen  ·  将列乘以向量的tidyverse解决方案
    3 年前
    Nuller  ·  无法对嵌套数据中的列进行子集设置。框架
    3 年前
    chrischi  ·  R Purrr-系数最高
    6 年前
    Richard Herron  ·  从purr的pmap()调用var()返回NA
    6 年前
    Davide Lorino  ·  嵌套数据帧上的函数向量算法
    6 年前
    ℕʘʘḆḽḘ  ·  如何使用purrr从两个元素的列表中提取元素?
    6 年前
    andrew_reece  ·  map\u dfr:填充。具有字符串标签的id列,而不是没有字符串标签的索引。x已命名
    6 年前
    关于   移动版
    代码之家 - 一站式码农服务社区
    沪ICP备11025650号