我如何计算特定位置多年来每个化学品的加权移动平均浓度？

Question

我正在尝试计算每个analyte_station的加权移动平均值。我已经在R中尝试了很多方法，但是一直得到NA返回值，我不知道为什么。

移动平均值将根据每个analyte_station的前5年数据(不包括本年度)计算。权重将是样本的数量或“计数”列。

因此，每个"ANALYTE_STATION“每年都有自己的加权平均值。

我一直得到的错误是:列Mean.5必须是长度1(组大小)，而不是35320

或者，该函数将返回整个列的NA值

我是R的新手，并且尝试过使用zoo和deplyr库以及滞后、变异和rollapply函数。

我经常遇到的一个错误是:列Mean.5的长度必须为1(组大小)，而不是35320

或者，该函数将返回整个列的NA值。

dput(head(mean2))
structure(list(
YEAR_ANALYTE_STATION = c("2006 4,4'-DDD CCE-01", 
"2007 4,4'-DDD CCE-01", "2008 4,4'-DDD CCE-01", "2009 4,4'-DDD CCE-01", 
"2010 4,4'-DDD CCE-01", "2013 4,4'-DDD CCE-01"), 
Year = c(2006L, 
2007L, 2008L, 2009L, 2010L, 2013L), 
CCEMean = c(1.96969696966667, 
0.635237880333333, 10.108880584, 8.91217270166667, 18.54267311, 
1.409054917), Count = c("3", "3", "3", "3", "3", "3"), 
ANALYTE_STATION = 
c("4,4'-DDD CCE-01", "4,4'-DDD CCE-01", "4,4'-DDD CCE-01", "4,4'-DDD CCE- 
01", "4,4'-DDD CCE-01", "4,4'-DDD CCE-01")), 
row.names = c(NA, -6L), 
class = 
c("grouped_df", "tbl_df", "tbl", "data.frame"), 
vars = "ANALYTE_STATION", 
drop = TRUE, indices = list(0:5), group_sizes = 6L, biggest_group_size = 6L, 
labels = structure(list(ANALYTE_STATION = "4,4'-DDD CCE-01"), row.names = 
c(NA, -1L), class = "data.frame", vars = "ANALYTE_STATION", drop = TRUE))

head(mean2)
# A tibble: 6 x 5
# Groups:   ANALYTE_STATION [1]
YEAR_ANALYTE_STATION  Year CCEMean Count ANALYTE_STATION
<chr>                <int>   <dbl> <chr> <chr>          
1 2006 4,4'-DDD CCE-01  2006   1.97  3     4,4'-DDD CCE-01
2 2007 4,4'-DDD CCE-01  2007   0.635 3     4,4'-DDD CCE-01
3 2008 4,4'-DDD CCE-01  2008  10.1   3     4,4'-DDD CCE-01
4 2009 4,4'-DDD CCE-01  2009   8.91  3     4,4'-DDD CCE-01
5 2010 4,4'-DDD CCE-01  2010  18.5   3     4,4'-DDD CCE-01
6 2013 4,4'-DDD CCE-01  2013   1.41  3     4,4'-DDD CCE-01

我能够让下面的代码工作，但现在我想给滚动均值增加一个权重。

前几年的移动平均值(不包括本年度)(例如5表示前5年的平均值)

mean5 = mean %>%
  mutate(Mean.lag1 = lag(Mean, n = 1)) %>%
  mutate(Mean.5.previous = rollapply(data = Mean.lag1, 
                                 width = 5, 
                                 FUN = mean, 
                                 align = "right", 
                                 fill = NA, 
                                 na.rm = T))

最终，每个"ANALYTE_STATION“都会有它自己的加权平均值(来自前5年)。

Answer 1

你可以使用一个自定义函数，它计算过去5年的平均值而不是行数的平均值。

myRollmean <- function(x, tm, rge, excl.last=0) sapply(1:length(x), function(i) {
  period <- (tm[i] - rge - excl.last):(tm[i] - excl.last)
  return(mean(x[which(tm %in% period)]))
})

x表示我们感兴趣的列。tm是时间列，rge是时间范围，即5年，excl.last表示应该排除多少年。

df1$Mean.5 <- myRollmean(df1$CCEMean, df1$Year, 5)
df1$Mean.5.previous <- myRollmean(df1$CCEMean, df1$Year, 5, 1)

结果

> df1
  YEAR_ANALYTE_STATION Year    CCEMean Count ANALYTE_STATION   Mean.5 Mean.5.previous
1 2006 4,4'-DDD CCE-01 2006  1.9696970     3 4,4'-DDD CCE-01 1.969697             NaN
2 2007 4,4'-DDD CCE-01 2007  0.6352379     3 4,4'-DDD CCE-01 1.302467        1.969697
3 2008 4,4'-DDD CCE-01 2008 10.1088806     3 4,4'-DDD CCE-01 4.237938        1.302467
4 2009 4,4'-DDD CCE-01 2009  8.9121727     3 4,4'-DDD CCE-01 5.406497        4.237938
5 2010 4,4'-DDD CCE-01 2010 18.5426731     3 4,4'-DDD CCE-01 8.033732        5.406497
6 2013 4,4'-DDD CCE-01 2013  1.4090549     3 4,4'-DDD CCE-01 9.743195        9.549741

您能检查一下函数是否产生了预期的结果吗？

Data

mean2 <- structure(list(YEAR_ANALYTE_STATION = c("2006 4,4'-DDD CCE-01", 
"2007 4,4'-DDD CCE-01", "2008 4,4'-DDD CCE-01", "2009 4,4'-DDD CCE-01", 
"2010 4,4'-DDD CCE-01", "2013 4,4'-DDD CCE-01"), Year = c(2006L, 
2007L, 2008L, 2009L, 2010L, 2013L), CCEMean = c(1.96969696966667, 
0.635237880333333, 10.108880584, 8.91217270166667, 18.54267311, 
1.409054917), Count = c("3", "3", "3", "3", "3", "3"), ANALYTE_STATION = c("4,4'-DDD CCE-01", 
"4,4'-DDD CCE-01", "4,4'-DDD CCE-01", "4,4'-DDD CCE-01", "4,4'-DDD CCE-01", 
"4,4'-DDD CCE-01")), row.names = c(NA, -6L), class = c("grouped_df", 
"tbl_df", "tbl", "data.frame"), vars = "ANALYTE_STATION", drop = TRUE, indices = list(
    0:5), group_sizes = 6L, biggest_group_size = 6L, labels = structure(list(
    ANALYTE_STATION = "4,4'-DDD CCE-01"), row.names = c(NA, -1L
), class = "data.frame", vars = "ANALYTE_STATION", drop = TRUE))

df1 <- as.data.frame(mean2)

Answer 2

假设问题说明问题末尾的代码给出了所需的未加权答案，则似乎需要的是数据帧的前5行的平均值。实际上，该代码不可能正确，因为6意味着它产生除了最后一行以外的所有行都应该是NA，因为只有最后一行有5个前面的行；然而，实际上，代码为第5行给出了一个不正确的非NA值-它应该是NA，因为在第5行之前只有4行，而不是5行。

要解决此问题，请尝试使用以下代码。注意，在?rollapply中，width参数可以是一个单元素的偏移量列表，因此list(-seq(5))表示前面的5个元素，还请注意，可以使用结尾带有r的rollapplyr来代替align="right"。在问题所示的数据中，CCEMean列中没有NAs，如果是这样的话，我们就不需要na.rm参数了，尽管这不会有什么坏处。

mean2 %>%
  mutate(Mean5prev = rollapplyr(CCEMean, list(-seq(5)), mean, fill = NA, na.rm = TRUE))

这个问题并没有精确地定义要使用的权重，而只是将上面的mean替换为一个函数，该函数可以以任何您想要的方式计算加权平均值。为此，尝试适当地使用weighted.mean函数。