面板数据的滚动平均值(有一些细节)
Rolling average for panel data (with a few details)
我想出了一些代码来计算面板数据滚动平均值(数据中的一行包含一天中一个主题的值)。由于我有一些更具体的要求,代码变得相当复杂。在我看来,对于一个不太罕见的应用来说太复杂了。
这是我需要的:
滚动平均值(前 3 天 (a) 值的平均值,不包括 "current" 天,( b)只有在这个window)
中至少有2个非缺失值才计算
尊重面板结构
不会太复杂吧?
对于 1. 我决定使用 rollapplyr()
和 mean( , na.rm = T)
来排除当天.对于 2。为了尊重面板结构,我将所有内容都包装在 tapply()
(使用 unlist()
)中。
代码示例如下:
library(zoo)
# example data (with missings)
set.seed(1)
df = data.frame(subject = rep(c("a", "b"), each = 10), day = rep(1:10, 2), value = rnorm(20))
df$value[15:17] = NA
# lag function (sensitive to "single day" subjects)
lag <- function(x, l = 1) {
if (length(x) > 1) (c(rep(NA, l), x[1:(length(x)-l)])) else (NA)
}
# calculate rolling mean
df$roll_mean3 = unlist(tapply(df$value, df$subject,
FUN = function(x) lag(rollapplyr(x, width = 3, fill = NA, partial = T,
FUN = function(x) ifelse(sum(!is.na(x)) > 1, mean(x, na.rm = T), NA)))))
df
正如我所说,对于我认为不远的情况,此解决方案似乎过于复杂。
关于如何以更简单(不易出错)的方式执行此操作,您有什么建议吗?
我是否错过了一些可以更轻松地处理面板数据的基本功能?
为了说明,我的代码的输出是:
subject day value roll_mean3
1 a 1 -0.6264538 NA
2 a 2 0.1836433 NA
3 a 3 -0.8356286 -0.221405243
4 a 4 1.5952808 -0.426146366
5 a 5 0.3295078 0.314431838
6 a 6 -0.8204684 0.363053321
7 a 7 0.4874291 0.368106730
8 a 8 0.7383247 -0.001177187
9 a 9 0.5757814 0.135095124
10 a 10 -0.3053884 0.600511703
11 b 1 1.5117812 NA
12 b 2 0.3898432 NA
13 b 3 -0.6212406 0.950812202
14 b 4 -2.2146999 0.426794608
15 b 5 NA -0.815365744
16 b 6 NA -1.417970234
17 b 7 NA NA
18 b 8 0.9438362 NA
19 b 9 0.8212212 NA
20 b 10 0.5939013 0.882528703
除了我上面的评论,我不完全确定你期望的输出应该是什么,但也许以下是一个很好的起点:
df %>%
group_by(subject) %>%
mutate(roll_mean3 = rollapplyr(
lag(value),
width = 3,
fill = NA,
FUN = function(x) ifelse(sum(!is.na(x)) > 1, mean(x, na.rm = T), NA)))
## A tibble: 20 x 4
## Groups: subject [2]
# subject day value roll_mean3
# <fct> <int> <dbl> <dbl>
# 1 a 1 -0.626 NA
# 2 a 2 0.184 NA
# 3 a 3 -0.836 -0.221
# 4 a 4 1.60 -0.426
# 5 a 5 0.330 0.314
# 6 a 6 -0.820 0.363
# 7 a 7 0.487 0.368
# 8 a 8 0.738 -0.00118
# 9 a 9 0.576 0.135
#10 a 10 -0.305 0.601
#11 b 1 1.51 NA
#12 b 2 0.390 NA
#13 b 3 -0.621 0.951
#14 b 4 -2.21 0.427
#15 b 5 NA -0.815
#16 b 6 NA -1.42
#17 b 7 NA NA
#18 b 8 0.944 NA
#19 b 9 0.821 NA
#20 b 10 0.594 0.883
或使用data.table
custom_mean <- function(x) ifelse(sum(!is.na(x)) > 1, mean(x, na.rm = T), NA)
setDT(df)[, roll_mean3 := rollapplyr(shift(value), width = 3, fill = NA, FUN = custom_mean), by = subject]
df
# subject day value roll_mean3
#1: a 1 -0.6264538 NA
#2: a 2 0.1836433 NA
#3: a 3 -0.8356286 -0.221405243
#4: a 4 1.5952808 -0.426146366
#5: a 5 0.3295078 0.314431838
#6: a 6 -0.8204684 0.363053321
#7: a 7 0.4874291 0.368106730
#8: a 8 0.7383247 -0.001177187
#9: a 9 0.5757814 0.135095124
#10: a 10 -0.3053884 0.600511703
#11: b 1 1.5117812 NA
#12: b 2 0.3898432 NA
#13: b 3 -0.6212406 0.950812202
#14: b 4 -2.2146999 0.426794608
#15: b 5 NA -0.815365744
#16: b 6 NA -1.417970234
#17: b 7 NA NA
#18: b 8 0.9438362 NA
#19: b 9 0.8212212 NA
#20: b 10 0.5939013 0.882528703
这可能不是最优雅或可扩展的解决方案,但它确实提供了所需的结果:
df %>%
group_by(subject) %>%
mutate(n_values = 3 - is.na(lag(value, 1)) - is.na(lag(value, 2)) - is.na(lag(value, 3)),
roll_mean = ifelse(
n_values >= 2,
(coalesce(lag(value), 0) + coalesce(lag(value, 2), 0) + coalesce(lag(value, 3), 0)) / n_values,
NA)
)
说明:这是一个 dplyr
管道,它首先按主题分组,以便尊重分组。接下来,mutate
中有两个计算值:
n_values
统计前三行非NA值的个数,每一个NA值等于3减1。使用 lag
.
访问前面的行
roll_mean
是有条件的,使用ifelse
:如果n_values
至少等于2,则可以计算均值。它将前 3 个值相加,使用 coalesce
将 NA 替换为 0。总和除以 n_values
得到平均值。如果n_values < 2
,返回NA。
对每个主题分别使用ave
到运行rollapply
。然后,当使用 rollapply
时,请注意 width
可以是一个包含偏移向量(或向量)的列表,因此 list(-seq(3))
表示前 3 个元素。有关参数的更多信息,请参阅 ?rollapply
。
Mean <- function(x) if (sum(!is.na(x)) >= 2) mean(x, na.rm = TRUE) else NA
roll <- function(x) rollapply(x, list(-seq(3)), Mean, fill = NA, partial = TRUE)
transform(df, roll = ave(value, subject, FUN = roll))
我想出了一些代码来计算面板数据滚动平均值(数据中的一行包含一天中一个主题的值)。由于我有一些更具体的要求,代码变得相当复杂。在我看来,对于一个不太罕见的应用来说太复杂了。
这是我需要的:
滚动平均值(前 3 天 (a) 值的平均值,不包括 "current" 天,( b)只有在这个window)
中至少有2个非缺失值才计算
尊重面板结构
不会太复杂吧?
对于 1. 我决定使用 rollapplyr()
和 mean( , na.rm = T)
来排除当天.对于 2。为了尊重面板结构,我将所有内容都包装在 tapply()
(使用 unlist()
)中。
代码示例如下:
library(zoo)
# example data (with missings)
set.seed(1)
df = data.frame(subject = rep(c("a", "b"), each = 10), day = rep(1:10, 2), value = rnorm(20))
df$value[15:17] = NA
# lag function (sensitive to "single day" subjects)
lag <- function(x, l = 1) {
if (length(x) > 1) (c(rep(NA, l), x[1:(length(x)-l)])) else (NA)
}
# calculate rolling mean
df$roll_mean3 = unlist(tapply(df$value, df$subject,
FUN = function(x) lag(rollapplyr(x, width = 3, fill = NA, partial = T,
FUN = function(x) ifelse(sum(!is.na(x)) > 1, mean(x, na.rm = T), NA)))))
df
正如我所说,对于我认为不远的情况,此解决方案似乎过于复杂。
关于如何以更简单(不易出错)的方式执行此操作,您有什么建议吗? 我是否错过了一些可以更轻松地处理面板数据的基本功能?
为了说明,我的代码的输出是:
subject day value roll_mean3
1 a 1 -0.6264538 NA
2 a 2 0.1836433 NA
3 a 3 -0.8356286 -0.221405243
4 a 4 1.5952808 -0.426146366
5 a 5 0.3295078 0.314431838
6 a 6 -0.8204684 0.363053321
7 a 7 0.4874291 0.368106730
8 a 8 0.7383247 -0.001177187
9 a 9 0.5757814 0.135095124
10 a 10 -0.3053884 0.600511703
11 b 1 1.5117812 NA
12 b 2 0.3898432 NA
13 b 3 -0.6212406 0.950812202
14 b 4 -2.2146999 0.426794608
15 b 5 NA -0.815365744
16 b 6 NA -1.417970234
17 b 7 NA NA
18 b 8 0.9438362 NA
19 b 9 0.8212212 NA
20 b 10 0.5939013 0.882528703
除了我上面的评论,我不完全确定你期望的输出应该是什么,但也许以下是一个很好的起点:
df %>%
group_by(subject) %>%
mutate(roll_mean3 = rollapplyr(
lag(value),
width = 3,
fill = NA,
FUN = function(x) ifelse(sum(!is.na(x)) > 1, mean(x, na.rm = T), NA)))
## A tibble: 20 x 4
## Groups: subject [2]
# subject day value roll_mean3
# <fct> <int> <dbl> <dbl>
# 1 a 1 -0.626 NA
# 2 a 2 0.184 NA
# 3 a 3 -0.836 -0.221
# 4 a 4 1.60 -0.426
# 5 a 5 0.330 0.314
# 6 a 6 -0.820 0.363
# 7 a 7 0.487 0.368
# 8 a 8 0.738 -0.00118
# 9 a 9 0.576 0.135
#10 a 10 -0.305 0.601
#11 b 1 1.51 NA
#12 b 2 0.390 NA
#13 b 3 -0.621 0.951
#14 b 4 -2.21 0.427
#15 b 5 NA -0.815
#16 b 6 NA -1.42
#17 b 7 NA NA
#18 b 8 0.944 NA
#19 b 9 0.821 NA
#20 b 10 0.594 0.883
或使用data.table
custom_mean <- function(x) ifelse(sum(!is.na(x)) > 1, mean(x, na.rm = T), NA)
setDT(df)[, roll_mean3 := rollapplyr(shift(value), width = 3, fill = NA, FUN = custom_mean), by = subject]
df
# subject day value roll_mean3
#1: a 1 -0.6264538 NA
#2: a 2 0.1836433 NA
#3: a 3 -0.8356286 -0.221405243
#4: a 4 1.5952808 -0.426146366
#5: a 5 0.3295078 0.314431838
#6: a 6 -0.8204684 0.363053321
#7: a 7 0.4874291 0.368106730
#8: a 8 0.7383247 -0.001177187
#9: a 9 0.5757814 0.135095124
#10: a 10 -0.3053884 0.600511703
#11: b 1 1.5117812 NA
#12: b 2 0.3898432 NA
#13: b 3 -0.6212406 0.950812202
#14: b 4 -2.2146999 0.426794608
#15: b 5 NA -0.815365744
#16: b 6 NA -1.417970234
#17: b 7 NA NA
#18: b 8 0.9438362 NA
#19: b 9 0.8212212 NA
#20: b 10 0.5939013 0.882528703
这可能不是最优雅或可扩展的解决方案,但它确实提供了所需的结果:
df %>%
group_by(subject) %>%
mutate(n_values = 3 - is.na(lag(value, 1)) - is.na(lag(value, 2)) - is.na(lag(value, 3)),
roll_mean = ifelse(
n_values >= 2,
(coalesce(lag(value), 0) + coalesce(lag(value, 2), 0) + coalesce(lag(value, 3), 0)) / n_values,
NA)
)
说明:这是一个 dplyr
管道,它首先按主题分组,以便尊重分组。接下来,mutate
中有两个计算值:
n_values
统计前三行非NA值的个数,每一个NA值等于3减1。使用lag
. 访问前面的行
roll_mean
是有条件的,使用ifelse
:如果n_values
至少等于2,则可以计算均值。它将前 3 个值相加,使用coalesce
将 NA 替换为 0。总和除以n_values
得到平均值。如果n_values < 2
,返回NA。
对每个主题分别使用ave
到运行rollapply
。然后,当使用 rollapply
时,请注意 width
可以是一个包含偏移向量(或向量)的列表,因此 list(-seq(3))
表示前 3 个元素。有关参数的更多信息,请参阅 ?rollapply
。
Mean <- function(x) if (sum(!is.na(x)) >= 2) mean(x, na.rm = TRUE) else NA
roll <- function(x) rollapply(x, list(-seq(3)), Mean, fill = NA, partial = TRUE)
transform(df, roll = ave(value, subject, FUN = roll))