三个月内每五个月滚动一次
Rolling mean every five months over three months
我想计算滚动平均值，规格如下：
- 从给定月份的月底开始，例如五月
- 使用最近三个月的（每日）数据计算这段时间的平均值
- 注意：特定月份的某些日期可能存在缺失值，并且每个月的天数可能不同，这使得每次计算的观测值数量通常可变
- 向前 5 个月重复此计算，例如如果在 5 月是最后一次计算，在 10 月底等，因此 window 每 5 个月滑动一次，并使用最近 3 个可用月份的数据 [假设数据从 2018 年 3 月开始，首先 window 将是：3 月-4 月-5 月 18 日，然后是 8 月-9 月-10 月 18 日，依此类推]
- dataset/memory 的大小对我来说很重要，因为我的真实数据集非常大
找了半天也没找到明确的解决办法，当width参数是可变的，window是滑动的。我特别在 zoo 中寻找解决方案。 datatable 和 plyr（或 xts）对于完成也很有趣。
示例数据（注意：这里没有缺失值，因为我不能轻易删除数据表中的行）
set.seed(44)
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")))
数据集包含 2 个个体 A 和 B 的数据点 'x'，可用于计算平均值。
下面我们使用最后注释中显示的数据，而不是问题中的样本数据。
1) 2 rollapply创建一个year/month变量ym然后对每个ID的值求和year/month也统计个数每个 ID 中的值和 year/month。然后对总和进行滚动总和，然后除以相应的计数滚动总和，按 ID 进行计算。
library(data.table)
library(zoo)
ym <- as.yearmon(dataset$time)
roll <- function(x) rollapplyr(x, 3, by = 5, sum, fill = NA)
ds <- na.omit(dataset[, list(x = sum(x), n = .N), by = list(ID, time = ym)][
, list(time, mean = roll(x) / roll(n)), by = ID])
给予：
> ds
ID time mean
1: A May 1988 -0.118017121
2: A Oct 1988 -0.045631016
3: A Mar 1989 -0.035498703
4: A Aug 1989 -0.055121507
5: A Jan 1990 0.018735210
6: A Jun 1990 0.091084791
7: A Nov 1990 -0.183955430
8: A Apr 1991 0.011909178
9: A Sep 1991 -0.040233435
10: A Feb 1992 0.051567634
11: A Jul 1992 0.006015941
12: A Dec 1992 0.253320798
13: A May 1993 -0.037722177
14: A Oct 1993 -0.145811906
15: A Mar 1994 0.134181429
16: A Aug 1994 -0.119081185
17: A Jan 1995 0.001921224
18: A Jun 1995 0.232193754
19: A Nov 1995 -0.077158954
20: A Apr 1996 -0.070271862
21: A Sep 1996 0.033858600
22: A Feb 1997 -0.053623676
23: A Jul 1997 -0.201388554
24: A Dec 1997 0.051488747
25: A May 1998 -0.073193772
26: A Oct 1998 -0.094019699
27: A Mar 1999 -0.078863959
28: A Aug 1999 0.110231533
29: A Jan 2000 0.141657202
30: B May 1988 0.130180515
31: B Oct 1988 0.025095818
32: B Mar 1989 -0.032415997
33: B Aug 1989 0.041286368
34: B Jan 1990 0.219208544
35: B Jun 1990 -0.023717715
36: B Nov 1990 -0.049073449
37: B Apr 1991 -0.051479646
38: B Sep 1991 0.124340203
39: B Feb 1992 0.040786822
40: B Jul 1992 0.019159682
41: B Dec 1992 0.083195470
42: B May 1993 0.006695704
43: B Oct 1993 0.119093846
44: B Mar 1994 0.077608445
45: B Aug 1994 0.132860266
46: B Jan 1995 -0.225050074
47: B Jun 1995 -0.091877628
48: B Nov 1995 -0.157798169
49: B Apr 1996 -0.219238136
50: B Sep 1996 0.289506566
51: B Feb 1997 0.118216626
52: B Jul 1997 0.186950994
53: B Dec 1997 -0.035447587
54: B May 1998 -0.159754318
55: B Oct 1998 -0.066470703
56: B Mar 1999 0.230782925
57: B Aug 1999 -0.052620748
58: B Jan 2000 -0.190938190
ID time mean
2) 1 rollapply 下面是上面的变体。它使用 by.column = FALSE 以便 mean2 可以同时处理 x 和 n。
library(data.table)
library(zoo)
ym <- as.yearmon(dataset$time)
mean2 <- function(xn) sum(xn[, 1]) / sum(xn[, 2])
roll2 <- function(x) rollapplyr(x, 3, by = 5, mean2, by.column = FALSE, fill = NA)
ds2 <- na.omit(dataset[, list(x = sum(x), n = .N), by = list(ID, time = ym)][
, list(time, mean = roll2(.SD)), .SDcols = c("x", "n"), by = ID])
3) 矢量宽度
我们可以定义一个矢量宽度并像这样滚动应用它。我们将宽度设置为大于那些不在月底的日期的元素数量的数字，以便它不计算这些日期的平均值。然后我们计算每个月末的平均值，并在最后一行代码中将其细分为每 5 个月一次。
library(data.table)
library(zoo)
ds3 <- dataset[, list(ID, time = as.yearmon(time), x)][,
list(time, x, width = seq_len(.N) - match(time - 2/12, time) + 1,
is_last = !duplicated(time, fromLast = TRUE)), by = ID][,
list(time, x, width = na.fill(ifelse(is_last, width, .N + 1), .N+1)), by = ID][,
list(time, mean = rollapplyr(x, width, mean, fill = NA_real_)),
by = ID][, na.omit(.SD)[seq(1, .N, 5), ], by = ID]
4) data.table join 这使用了 data.table join 而不是 rollapply。 eom 是 data.table 仅包含月末行。它还有一列 time2 代表 2 个月前的年月。我们将其与 datasetym 连接起来并提取适当的行和列。
library(data.table)
library(zoo)
datasetym <- dataset[, list(ID, time = as.yearmon(time), x)]
eom <- datasetym[, .SD[!duplicated(time, fromLast = TRUE), ], by = ID][
, cbind(.SD, time2 = time - 2/12)]
ds4 <- datasetym[eom, list(mean = mean(x)),
on = .(ID, time >= time2, time <= time), by = .EACHI][
, .SD[seq(3, .N, 5), -2], by = ID]
5) sqldf 您可能更喜欢使用更熟悉的 SQL 语法来表达连接。创建 datasetym 并按照 (4) 中的方式完成每第 5 行。
library(data.table)
library(sqldf)
library(zoo)
datasetym <- dataset[, list(ID, time = as.yearmon(time), x)]
s <- sqldf("select a.ID, a.time, avg(b.x) mean
from (select ID, time from datasetym group by ID, time) a
left join datasetym b
on a.ID = b.ID and b.time between a.time - 2.0/12.0 and a.time
group by a.ID, a.time")
ds5 <- data.table(s)[, .SD[seq(3, .N, 5), ], by = ID]
6) zoo 如果我们使用宽格式，我们可以只使用 zoo 来解决这个问题。如果需要，我们总是可以在之后转换回长格式（如注释行所示）。
library(zoo)
z <- read.zoo(dataset, index = "time", split = "ID")
zsum <- aggregate(z, as.yearmon, sum)
zlength <- aggregate(z, as.yearmon, length)
zroll <- rollapplyr(zsum, 3, by = 5, sum) / rollapplyr(zlength, 3, by = 5, sum)
# fortify(zroll, melt = TRUE) # if long form wanted
给予：
> zroll
A B
May 1988 -0.118017121 0.130180515
Oct 1988 -0.045631016 0.025095818
Mar 1989 -0.035498703 -0.032415997
Aug 1989 -0.055121507 0.041286368
Jan 1990 0.018735210 0.219208544
Jun 1990 0.091084791 -0.023717715
Nov 1990 -0.183955430 -0.049073449
Apr 1991 0.011909178 -0.051479646
Sep 1991 -0.040233435 0.124340203
Feb 1992 0.051567634 0.040786822
Jul 1992 0.006015941 0.019159682
Dec 1992 0.253320798 0.083195470
May 1993 -0.037722177 0.006695704
Oct 1993 -0.145811906 0.119093846
Mar 1994 0.134181429 0.077608445
Aug 1994 -0.119081185 0.132860266
Jan 1995 0.001921224 -0.225050074
Jun 1995 0.232193754 -0.091877628
Nov 1995 -0.077158954 -0.157798169
Apr 1996 -0.070271862 -0.219238136
Sep 1996 0.033858600 0.289506566
Feb 1997 -0.053623676 0.118216626
Jul 1997 -0.201388554 0.186950994
Dec 1997 0.051488747 -0.035447587
May 1998 -0.073193772 -0.159754318
Oct 1998 -0.094019699 -0.066470703
Mar 1999 -0.078863959 0.230782925
Aug 1999 0.110231533 -0.052620748
Jan 2000 0.141657202 -0.190938190
备注
请注意，问题中定义的 dataset 有 8832 行，但用于定义 ID 列的向量只有 4416 个元素，因此它被回收，结果前 2216 个日期在 A 中结束两次并且在 B 中根本没有，接下来的 2216 个日期在 B 中结束两次，而在 A 中根本没有。大概这不是预期的，我们通过在数据集定义中将每次出现的 2208 替换为 4416 来解决这个问题，以便每个日期在 A 中出现一次，在 B 中出现一次：
set.seed(44)
dataset <- data.table(ID = c(rep("A", 4416), rep("B", 4416)),
x = rnorm(4416 * 2),
time = c(seq(as.Date("1988/03/15"), as.Date("2000/04/16"), "day")))
我想计算滚动平均值，规格如下：
- 从给定月份的月底开始，例如五月
- 使用最近三个月的（每日）数据计算这段时间的平均值
- 注意：特定月份的某些日期可能存在缺失值，并且每个月的天数可能不同，这使得每次计算的观测值数量通常可变
- 向前 5 个月重复此计算，例如如果在 5 月是最后一次计算，在 10 月底等，因此 window 每 5 个月滑动一次，并使用最近 3 个可用月份的数据 [假设数据从 2018 年 3 月开始，首先 window 将是：3 月-4 月-5 月 18 日，然后是 8 月-9 月-10 月 18 日，依此类推]
- dataset/memory 的大小对我来说很重要，因为我的真实数据集非常大
找了半天也没找到明确的解决办法，当width参数是可变的，window是滑动的。我特别在 zoo 中寻找解决方案。 datatable 和 plyr（或 xts）对于完成也很有趣。
示例数据（注意：这里没有缺失值，因为我不能轻易删除数据表中的行）
set.seed(44)
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")))
数据集包含 2 个个体 A 和 B 的数据点 'x'，可用于计算平均值。
下面我们使用最后注释中显示的数据，而不是问题中的样本数据。
1) 2 rollapply创建一个year/month变量ym然后对每个ID的值求和year/month也统计个数每个 ID 中的值和 year/month。然后对总和进行滚动总和，然后除以相应的计数滚动总和，按 ID 进行计算。
library(data.table)
library(zoo)
ym <- as.yearmon(dataset$time)
roll <- function(x) rollapplyr(x, 3, by = 5, sum, fill = NA)
ds <- na.omit(dataset[, list(x = sum(x), n = .N), by = list(ID, time = ym)][
, list(time, mean = roll(x) / roll(n)), by = ID])
给予：
> ds
ID time mean
1: A May 1988 -0.118017121
2: A Oct 1988 -0.045631016
3: A Mar 1989 -0.035498703
4: A Aug 1989 -0.055121507
5: A Jan 1990 0.018735210
6: A Jun 1990 0.091084791
7: A Nov 1990 -0.183955430
8: A Apr 1991 0.011909178
9: A Sep 1991 -0.040233435
10: A Feb 1992 0.051567634
11: A Jul 1992 0.006015941
12: A Dec 1992 0.253320798
13: A May 1993 -0.037722177
14: A Oct 1993 -0.145811906
15: A Mar 1994 0.134181429
16: A Aug 1994 -0.119081185
17: A Jan 1995 0.001921224
18: A Jun 1995 0.232193754
19: A Nov 1995 -0.077158954
20: A Apr 1996 -0.070271862
21: A Sep 1996 0.033858600
22: A Feb 1997 -0.053623676
23: A Jul 1997 -0.201388554
24: A Dec 1997 0.051488747
25: A May 1998 -0.073193772
26: A Oct 1998 -0.094019699
27: A Mar 1999 -0.078863959
28: A Aug 1999 0.110231533
29: A Jan 2000 0.141657202
30: B May 1988 0.130180515
31: B Oct 1988 0.025095818
32: B Mar 1989 -0.032415997
33: B Aug 1989 0.041286368
34: B Jan 1990 0.219208544
35: B Jun 1990 -0.023717715
36: B Nov 1990 -0.049073449
37: B Apr 1991 -0.051479646
38: B Sep 1991 0.124340203
39: B Feb 1992 0.040786822
40: B Jul 1992 0.019159682
41: B Dec 1992 0.083195470
42: B May 1993 0.006695704
43: B Oct 1993 0.119093846
44: B Mar 1994 0.077608445
45: B Aug 1994 0.132860266
46: B Jan 1995 -0.225050074
47: B Jun 1995 -0.091877628
48: B Nov 1995 -0.157798169
49: B Apr 1996 -0.219238136
50: B Sep 1996 0.289506566
51: B Feb 1997 0.118216626
52: B Jul 1997 0.186950994
53: B Dec 1997 -0.035447587
54: B May 1998 -0.159754318
55: B Oct 1998 -0.066470703
56: B Mar 1999 0.230782925
57: B Aug 1999 -0.052620748
58: B Jan 2000 -0.190938190
ID time mean
2) 1 rollapply 下面是上面的变体。它使用 by.column = FALSE 以便 mean2 可以同时处理 x 和 n。
library(data.table)
library(zoo)
ym <- as.yearmon(dataset$time)
mean2 <- function(xn) sum(xn[, 1]) / sum(xn[, 2])
roll2 <- function(x) rollapplyr(x, 3, by = 5, mean2, by.column = FALSE, fill = NA)
ds2 <- na.omit(dataset[, list(x = sum(x), n = .N), by = list(ID, time = ym)][
, list(time, mean = roll2(.SD)), .SDcols = c("x", "n"), by = ID])
3) 矢量宽度
我们可以定义一个矢量宽度并像这样滚动应用它。我们将宽度设置为大于那些不在月底的日期的元素数量的数字，以便它不计算这些日期的平均值。然后我们计算每个月末的平均值，并在最后一行代码中将其细分为每 5 个月一次。
library(data.table)
library(zoo)
ds3 <- dataset[, list(ID, time = as.yearmon(time), x)][,
list(time, x, width = seq_len(.N) - match(time - 2/12, time) + 1,
is_last = !duplicated(time, fromLast = TRUE)), by = ID][,
list(time, x, width = na.fill(ifelse(is_last, width, .N + 1), .N+1)), by = ID][,
list(time, mean = rollapplyr(x, width, mean, fill = NA_real_)),
by = ID][, na.omit(.SD)[seq(1, .N, 5), ], by = ID]
4) data.table join 这使用了 data.table join 而不是 rollapply。 eom 是 data.table 仅包含月末行。它还有一列 time2 代表 2 个月前的年月。我们将其与 datasetym 连接起来并提取适当的行和列。
library(data.table)
library(zoo)
datasetym <- dataset[, list(ID, time = as.yearmon(time), x)]
eom <- datasetym[, .SD[!duplicated(time, fromLast = TRUE), ], by = ID][
, cbind(.SD, time2 = time - 2/12)]
ds4 <- datasetym[eom, list(mean = mean(x)),
on = .(ID, time >= time2, time <= time), by = .EACHI][
, .SD[seq(3, .N, 5), -2], by = ID]
5) sqldf 您可能更喜欢使用更熟悉的 SQL 语法来表达连接。创建 datasetym 并按照 (4) 中的方式完成每第 5 行。
library(data.table)
library(sqldf)
library(zoo)
datasetym <- dataset[, list(ID, time = as.yearmon(time), x)]
s <- sqldf("select a.ID, a.time, avg(b.x) mean
from (select ID, time from datasetym group by ID, time) a
left join datasetym b
on a.ID = b.ID and b.time between a.time - 2.0/12.0 and a.time
group by a.ID, a.time")
ds5 <- data.table(s)[, .SD[seq(3, .N, 5), ], by = ID]
6) zoo 如果我们使用宽格式，我们可以只使用 zoo 来解决这个问题。如果需要，我们总是可以在之后转换回长格式（如注释行所示）。
library(zoo)
z <- read.zoo(dataset, index = "time", split = "ID")
zsum <- aggregate(z, as.yearmon, sum)
zlength <- aggregate(z, as.yearmon, length)
zroll <- rollapplyr(zsum, 3, by = 5, sum) / rollapplyr(zlength, 3, by = 5, sum)
# fortify(zroll, melt = TRUE) # if long form wanted
给予：
> zroll
A B
May 1988 -0.118017121 0.130180515
Oct 1988 -0.045631016 0.025095818
Mar 1989 -0.035498703 -0.032415997
Aug 1989 -0.055121507 0.041286368
Jan 1990 0.018735210 0.219208544
Jun 1990 0.091084791 -0.023717715
Nov 1990 -0.183955430 -0.049073449
Apr 1991 0.011909178 -0.051479646
Sep 1991 -0.040233435 0.124340203
Feb 1992 0.051567634 0.040786822
Jul 1992 0.006015941 0.019159682
Dec 1992 0.253320798 0.083195470
May 1993 -0.037722177 0.006695704
Oct 1993 -0.145811906 0.119093846
Mar 1994 0.134181429 0.077608445
Aug 1994 -0.119081185 0.132860266
Jan 1995 0.001921224 -0.225050074
Jun 1995 0.232193754 -0.091877628
Nov 1995 -0.077158954 -0.157798169
Apr 1996 -0.070271862 -0.219238136
Sep 1996 0.033858600 0.289506566
Feb 1997 -0.053623676 0.118216626
Jul 1997 -0.201388554 0.186950994
Dec 1997 0.051488747 -0.035447587
May 1998 -0.073193772 -0.159754318
Oct 1998 -0.094019699 -0.066470703
Mar 1999 -0.078863959 0.230782925
Aug 1999 0.110231533 -0.052620748
Jan 2000 0.141657202 -0.190938190
备注
请注意，问题中定义的 dataset 有 8832 行，但用于定义 ID 列的向量只有 4416 个元素，因此它被回收，结果前 2216 个日期在 A 中结束两次并且在 B 中根本没有，接下来的 2216 个日期在 B 中结束两次，而在 A 中根本没有。大概这不是预期的，我们通过在数据集定义中将每次出现的 2208 替换为 4416 来解决这个问题，以便每个日期在 A 中出现一次，在 B 中出现一次：
set.seed(44)
dataset <- data.table(ID = c(rep("A", 4416), rep("B", 4416)),
x = rnorm(4416 * 2),
time = c(seq(as.Date("1988/03/15"), as.Date("2000/04/16"), "day")))