变窗移动平均法

Question

我想计算一个变化窗口大小的R中的变量的移动平均值。更确切地说:移动平均应该在三年内计算，但是数据(时间序列)在更高的频率下可用，而且每个三年窗口的窗口大小都可能不同。

假设以下数据集：

library(data.table)
set.seed(1)   # reproduceable data
dataset <- data.table(ID=c(rep("A",2208),rep("B",2208)),
x = c(rnorm(2208*2)), time=c(seq(as.Date("1988/03/15"),
as.Date("2000/04/16"), "day"),seq(as.Date("1988/03/15"),
as.Date("2000/04/16"), "day")))

变量x的三年移动平均值应该为两个ID(个人A和B)计算。最好用zoo和datatable来完成这个任务吗？但任何解决办法都可以。

请注意，我知道如何使用一个固定的窗口大小，这里的问题是变化的窗口大小。

Answer 1

如果我理解正确的话，OP想要跨越整整3年。由于这可能包括闰年，窗口大小可以是1095天，也可以是1096天。

这可以通过在非赤道连接中聚合和lubridate的回滚日期算法来解决。

library(data.table)
library(lubridate)
# create 3 years windows for each ID for later non-equi join
win <- dataset[, CJ(ID = ID, start = time, unique = TRUE)][
  # make sure to pick
  , end := start %m+% years(3) - days(1)][
    # remove windows which end out of date range
    end <= max(start)]
win

ID start end 1: A 1988-03-15 1991-03-14 2: A 1988-03-16 1991-03-15 3: A 1988-03-17 1991-03-16 4: A 1988-03-18 1991-03-17 5: A 1988-03-19 1991-03-18 --- 6638: B 1997-04-13 2000-04-12 6639: B 1997-04-14 2000-04-13 6640: B 1997-04-15 2000-04-14 6641: B 1997-04-16 2000-04-15 6642: B 1997-04-17 2000-04-16

# check window lengths
win[, .N, by = .(days = end - start + 1L)]

days N 1: 1095 days 2166 2: 1096 days 4476

# see what happens in leap years
win[leap_year(start) & month(start) == 2 & day(start) %in% 28:29, 
  .(start, end, days = end - start + 1L)]

start end days 1: 1992-02-28 1995-02-27 1096 days 2: 1992-02-29 1995-02-27 1095 days 3: 1996-02-28 1999-02-27 1096 days 4: 1996-02-29 1999-02-27 1095 days 5: 1992-02-28 1995-02-27 1096 days 6: 1992-02-29 1995-02-27 1095 days 7: 1996-02-28 1999-02-27 1096 days 8: 1996-02-29 1999-02-27 1095 days

win[leap_year(end) & month(end) == 2 & day(end) %in% 28:29,
    .(start, end, days = end - start + 1L)]

start end days 1: 1989-03-01 1992-02-29 1096 days 2: 1993-03-01 1996-02-29 1096 days 3: 1997-03-01 2000-02-29 1096 days 4: 1989-03-01 1992-02-29 1096 days 5: 1993-03-01 1996-02-29 1096 days 6: 1997-03-01 2000-02-29 1096 days

# aggregate in a non-equi-join
dataset[win, on = .(ID, time >= start, time <= end), by = .EACHI, .(avg = mean(x))]

ID time time avg 1: A 1988-03-15 1991-03-14 -0.01184078 2: A 1988-03-16 1991-03-15 -0.01317813 3: A 1988-03-17 1991-03-16 -0.01179571 4: A 1988-03-18 1991-03-17 -0.01006100 5: A 1988-03-19 1991-03-18 -0.01221798 --- 6638: B 1997-04-13 2000-04-12 -0.03412214 6639: B 1997-04-14 2000-04-13 -0.03604176 6640: B 1997-04-15 2000-04-14 -0.03556291 6641: B 1997-04-16 2000-04-15 -0.03392185 6642: B 1997-04-17 2000-04-16 -0.03393674

Answer 2

正如@A.Suliman在注释中提到的，您的示例数据有固定的窗口宽度，但是让我们假设您的真实数据还没有，因为您的文本就是这么说的。

来自rollapply的参数rollapply不必是常数，所以我们可以先计算宽度，确保窗口与左边对齐，然后运行rollaply。

library(zoo)
library(tidyverse)
dataset %>% 
  arrange(ID,time) %>%
  group_by(ID) %>%
  mutate(avg = rollapply(x, FUN = mean, align = "left",
                         width = map_dbl(time, ~which.max(time[time < .x + 3*365.25])) - row_number()+1))
#                          
# # A tibble: 8,832 x 4
# # Groups:   ID [2]
#         ID       x       time    avg
#     <fctr>   <dbl>     <date>  <dbl>
#   1      A -0.0258 1988-03-15 0.0109
#   2      A -0.0258 1988-03-15 0.0109
#   3      A -0.1562 1988-03-16 0.0115
#   4      A -0.1562 1988-03-16 0.0115
#   5      A  0.8193 1988-03-17 0.0115
#   6      A  0.8193 1988-03-17 0.0112
#   7      A -1.1136 1988-03-18 0.0102
#   8      A -1.1136 1988-03-18 0.0107
#   9      A -0.9336 1988-03-19 0.0105
#  10      A -0.9336 1988-03-19 0.0109