Wednesday, 4 September 2013

R: Calculate moving maximum slope by week accounting for factors

R: Calculate moving maximum slope by week accounting for factors

I have a data.frame that includes heating degree day (HDD) below.
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")
I'm trying to calculate the moving maximum slope over 7 days beginning on
4 November of each year without using a loop. This moving maximum slope
needs to account for WinterID, SiteID, and SubstrateConcat.
For clarification, the calculation I'm trying to obtain is this:
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
So, using a real example starting with 2002-11-19 as i:
(0+0.24+7.13+10.61+2.97+0+0.97) - 0)/7 = 3.13
I tried using zoo package rollmean, however, I could not figure out how to
account for WinterID, SiteID, and SubstrateConcat. This gave me an
"order.by" error where my Date values were not unique, since I have dates
with different SubstrateConcat and WinterID criteria. As I enter more data
into the database, there will eventually be dates with multiple SiteID
criteria as well.
I thought maybe xts, TTR and ROC would be what I could use as in this
question: Maximum slope for a given interval each day. But again, I don't
understand how to specify the multiple group factors, as well as going
three days forward and three days back as in align=center with rollmean.
Will someone please point me in the right direction here? Will one of the
above functions combined with ddply work?
Thank you!

No comments:

Post a Comment