R：按周计算移动最大斜率，考虑因素答案

【问题标题】：R: Calculate moving maximum slope by week accounting for factorsR：按周计算移动最大斜率，考虑因素
【发布时间】：2013-09-04 22:34:02
【问题描述】：

我有一个 data.frame，其中包括下面的采暖度日 (HDD)。

structure(list(WinterID = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L), .Label = c("2002", "2002_2003", "2003", "2003_2004", 
"2004", "2004_2005", "2005", "2005_2006", "2006", "2006_2007", 
"2007", "2007_2008", "2008"), class = "factor"), Date = structure(c(11968, 
11969, 11970, 11971, 11972, 11973, 11974, 11975, 11976, 11977, 
11978, 11979, 11980, 11981, 11982, 11983, 11984, 11985, 11986, 
11987, 11988, 11989, 11990, 11991, 11992, 11993, 11994, 11995, 
11996, 11997, 11998, 11999, 12000, 12001, 12002, 12003, 12004, 
12005, 12006, 12007, 12008, 12009, 12010, 12011, 12012, 12013, 
12014, 12015, 12016, 12017, 12018, 12019, 12020, 12021, 12022, 
12023, 12024, 12025, 12026, 12027, 12028, 12029, 12030, 12031, 
12032, 12033, 12034, 12035, 12036, 12037, 12038, 12039, 12040, 
12041, 12042, 12043, 12044, 12045, 12046, 12047, 12048, 12049, 
12050, 12051, 12052, 12053, 12054, 12055, 12056, 12057, 12058, 
12059, 12060, 12061, 12062, 12063, 12064, 12065, 12066, 12067, 
12068, 12069, 12070, 12071, 12072, 12073, 12074, 12075, 12076, 
12077, 12078, 12079, 12080, 12081, 12082, 12083, 12084, 12085, 
12086, 12087, 12088, 12089, 12090, 12091, 12092, 12093, 12094, 
12095, 12096, 12097, 12098, 12099, 12100, 12101, 12102, 12103, 
12104, 12105, 12106, 12107, 12108, 12109, 12110, 12111, 12112, 
12113, 12114, 12115, 12116, 12117, 12118, 12119, 12120, 12121, 
12122, 12123, 12124, 12125, 12126, 12127, 12128, 12129, 12130, 
12131, 12132, 12133, 12134, 12135, 12136, 12137, 12138, 12139, 
12140, 12141, 12142, 12010, 12011, 12014, 12015, 12017, 12023, 
12024, 12025, 12026, 12027, 12028, 12029, 12030, 12042, 12070, 
12071, 12075, 12076, 12077, 12078, 12079, 12080, 12082, 12083, 
12084), class = "Date"), SiteID = structure(c(1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L), .Label = "NW_SB", class = "factor"), SubstrateConcat = structure(c(2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("B_A", "B_B", "B_E"), class = "factor"), 
    HDD = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0.246666666666667, 7.12666666666667, 10.6133333333333, 
    2.96666666666667, 0, 0.0933333333333337, 7.31333333333334, 
    10.7133333333333, 6.20000000000001, 2.70666666666667, 6.20000000000001, 
    3.88666666666667, 16.5866666666667, 28.3933333333333, 12.98, 
    21.6133333333333, 19.14, 12.6666666666667, 7.52, 3.33333333333334, 
    18.2933333333333, 4.14666666666667, 2.17333333333334, 26.08, 
    1.38, 7.48000000000001, 36.5733333333333, 53.4666666666667, 
    98.4533333333333, 109.093333333333, 104.14, 80.2466666666667, 
    47.0333333333333, 14.7133333333333, 15.7266666666667, 21.1066666666667, 
    5.07333333333334, 0.613333333333334, 6.18000000000001, 29.5666666666667, 
    45.5333333333333, 59.5666666666667, 91.44, 85.38, 51.1, 25.9666666666667, 
    14.8266666666667, 34.48, 79.16, 90.08, 66.3533333333333, 
    75.14, 97.1733333333333, 83.3066666666667, 50.0133333333333, 
    37.2733333333333, 88.9133333333334, 101.926666666667, 100.56, 
    99.2933333333334, 97.66, 89.6466666666667, 110.613333333333, 
    79.1466666666667, 92.6066666666667, 71.7133333333333, 31.32, 
    27.02, 39.02, 98.14, 62.5866666666667, 46.7933333333333, 
    47.5133333333333, 48.3666666666667, 25.5333333333333, 13.6, 
    17.9133333333333, 14.16, 7.98666666666667, 3.44, 1.86666666666667, 
    12.66, 0, 7.09333333333334, 21.3266666666667, 40.52, 18.8466666666667, 
    37.8466666666667, 33.42, 33.7133333333333, 15.6133333333333, 
    0.720000000000001, 2.31333333333334, 12.3066666666667, 8.48666666666667, 
    2.86, 0, 0, 0, 6.98666666666667, 6.67333333333334, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6.58000000000001, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 13.42, 30.5266666666667, 
    1.12, 28.5066666666667, 6.82666666666667, 10.3933333333333, 
    3.18, 11.0466666666667, 0, 0, 0)), .Names = c("WinterID", 
"Date", "SiteID", "SubstrateConcat", "HDD"), row.names = c(NA, 
200L), class = "data.frame")

我正在尝试计算从每年 11 月 4 日开始的 7 天内的最大移动斜率，而不使用循环。这个移动的最大坡度需要考虑 WinterID、SiteID 和 SubstrateConcat。

为了澄清，我试图获得的计算是这样的：

Slope=(max-min)/7, where:
Max= (i-3)+(i-2)+(i-1)+i+(i+1)+(i+2)+(i+3) 
Min= (i-3)

(((i-3)+(i-2)+(i-1)+i+(i+1)+(i+2)+(i+3)) - (i-3))/7

所以，使用从 2002-11-19 开始的真实示例作为 i：

(0+0.24+7.13+10.61+2.97+0+0.97) - 0)/7 = 3.13

我尝试使用zoo 包rollmean，但是，我不知道如何解释 WinterID、SiteID 和 SubstrateConcat。这给了我一个“order.by”错误，其中我的 Date 值不是唯一的，因为我的日期具有不同的 SubstrateConcat 和 WinterID 标准。随着我在数据库中输入更多数据，最终也会出现具有多个 SiteID 标准的日期。

我想也许xts、TTR 和ROC 将是我可以在这个问题中使用的：Maximum slope for a given interval each day。但同样，我不明白如何指定多个组因素，以及在 align=center 和 rollmean 中前进三天和三天后退。

有人可以在这里指出正确的方向吗？上述功能之一与ddply 结合是否有效？

谢谢！

编辑在@eddi 提供的答案之后包含答案。

dt <- data.table(df)
dt[, MaxSlope := if(length(HDD)<7) {rep(NA_real_, length(HDD))} else {filter(HDD, c(1,1,1,1,1,1,0)/7)}, by=list(Winter, Site, Substrate)]

此代码非常适用于连续日期。谁能推荐如何为缺少日期的数据调整此代码？例如，我有：

   Date  Temp 
 Nov 21  14 
 Nov 23  10 
 Nov 24  12 
 Nov 27  11 
 Nov 28  7 
 Nov 29  9 
 Nov 30  10 
 Dec 01  12 
 Dec 02  8  
 Dec 03  7

我不希望计算 11 月 21 日、11 月 23 日和 11 月 24 日的最大坡度，因为没有用于计算的连续数据。相反，我想插入“NA”。是否可以修改上述现有代码以适应这种情况？

【问题讨论】：

@naomik 几乎不是一个不赞成的理由 imo...
如果要粘贴很长的代码，至少尝试格式化它。但从表面上看，您并没有尝试减少代码粘贴的大小或隔离您的问题。
@naomik：OP 提供了可重现的数据。我认为您误解了带有长代码的dput :)

标签： r data.table

【解决方案1】：

听起来您需要filter（或者您也可以使用滚动均值/求和函数之一）。而分组部分最容易用data.table：

library(data.table)
dt = data.table(your_df)

dt[, filter(HDD, c(1,1,1,1,1,1,0))/7,
     by = list(WinterID, SiteID, SubstrateConcat)]

【讨论】：

你能解释一下过滤器在做什么吗？我认为 filter() 是子集的另一种方式。我不完全理解 "c(1,1,1,1,1,1,0))/7 正在做什么。它是在计算总和并除以 7 吗？我需要斜率，即 (Max- Min)/7.
另外，通过上面的代码，我看到我可以指定“sides=2”来获取中间值的前进和后退。但是，我如何指定中间值应该是什么？例如，如果我希望窗口从 11 月 4 日开始，这样我的窗口是 11 月 1 日（后 3 日）到 11 月 7 日（前 3 日），我该怎么做？
RBD.DT[, filter(HDD, c(1,1,1,1,1,1,0))/7, by=list(WinterID, SiteID, SubstrateConcat)] 过滤器错误（HDD, c(1, 1, 1, 1, 1, 1, 1)) ：“过滤器”比时间序列长
@BgnR 请参阅 ?filter 的详细信息部分以了解它的作用 - 在这种情况下，它只是根据 OP 添加 (i-2)+(i-1)+i+(i+1)+(i+2)+(i+3)；该错误仅意味着您有一个少于 7 个元素的子组 - 例如filter(c(1,2,3), c(1,1,1,1,1,1)) - 您可以添加 if 来说明这一点并采取您喜欢的任何操作，例如：dt[, if (length(HDD) < 7) { rep(NA_real_, length(HDD))} else {filter...}, ...]。
NA_real_ 是NA 的真实版本，我曾经明确表示它不是一个合乎逻辑的NA（默认值，当你输入NA 时）； rep 只是重复 NA length(HDD) 次
@BgnR 使用dt[, newcol := ...] 将该列添加到您的data.table

【解决方案2】：

我无法使用ddply 获得有效的解决方案，尽管我没有花太多时间进行调试。这是一个使用基函数的解决方案（假设您的对象名为hdd）。

# split your object into groups
shdd <- split(hdd, hdd[,c("WinterID","SiteID","SubstrateConcat")], drop=TRUE)
# create a function to apply to each group
f <- function(d) transform(d, MaxSlopeHDD=rollmax(c(NA,diff(d$HDD)),7,fill=NA))
# apply the function to each group and rbind the results together
shdd <- do.call(rbind, lapply(shdd, f))

【讨论】：