从 R 中的多年数据运行每个日历日的百分位值答案

【问题标题】：Running percentile value for each calendar day from multi-year data in R从 R 中的多年数据运行每个日历日的百分位值
【发布时间】：2023-01-19 07:23:02
【问题描述】：

我需要根据多年数据计算每个日历日的 30 天运行（窗口）第 90 个百分位最高温度值。例如，要计算 1 月 1 日的第 90 个百分位数值，我必须选择一个以 1 月 1 日为中心的 30 天窗口，即 12 月 16 日到 1 月 15 日的所有 42 年的数据。所以，我每天会有 1260 (30*42) 个数据点。我需要 366 天的值。我有从 1980 年到 2022 年的 42 年每日数据集，如下所示：

date    tmax    tmin
1981-01-01  19.2    5.4
1981-01-02  18.2    5
1981-01-03  16.1    3.8
1981-01-04  17.2    4.4
1981-01-05  15.7    2.4
1981-01-06  15.6    5.4
1981-01-07  11.2    4.1
1981-01-08  14.8    -1
1981-01-09  15  0.8
1981-01-10  16.2    -0.4

.........................
.........................
.........................
2022-12-25  17.4    4.4
2022-12-26  16.5    4.1
2022-12-27  17  5.4
2022-12-28  15.2    3.6
2022-12-29  8.1 7.7
2022-12-30  13.5    6
2022-12-31  14.8    4.5

我怎样才能在 R 中做到这一点？最初，我认为它会像这样简单。

temp_data <- read.csv("temperature.csv")

#as the date and tmax data are being read as characters by R
temp_data$tmax <- as.numeric(temp_data$tmax)
temp_data$date <- as.Date(temp_data$date, "%Y-%m-%d")
#Create a day of year variable for the day of the year
temp_data$doy <- as.numeric(format(temp_data$date,"%j"))

#load libraries
library(dplyr)
library(zoo)

temp_data_90th <- temp_data %>% 
  group_by(doy) %>% 
  summarize(rolling_90th = rollapply(tmax, width = 30, FUN = quantile, prob = 0.9, align = "center", na.rm=T))

但我不认为它给出了正确的结果，因为 temp_data_90th 有 4,470 行，一年中的每一天都有 13 个数据。

请你能建议我哪里做错了吗？预先感谢您对我们的支持。

【问题讨论】：

当您 group_by(doy) 时，您是在告诉 dplyr 将您的数据分解为每个 doy 值的单独组，并且仅执行后续操作之内那些小团体。您想要对 30 个唯一的 doy 值应用滚动函数，所以您肯定不要想group_by(doy)。我猜你可能也想要 mutate 而不是 summarize。
嗨@GregorThomas。即使我们在不分组的情况下使用 rollapply，它也会计算所有时间序列数据 (nrows = 15065) 的百分位值，而不是一年中的某一天。我需要一年中每一天的百分位值，即最终结果应该是 nrows = 366。
关键在于您的计算需要从不同的doy 值访问数据，而group_by(doy) 将使这成为不可能。您不能使用 width = 30，因为那将是 30 次观察，而您希望每年进行 30 天的观察。我认为 slider 包 might make this easier 使用 slide_index 函数，但我以前从未使用过它，所以我无法提供比快速指针更多的帮助。

标签： r percentile

【解决方案1】：

为了说明这一点，我们将需要可重现的数据，因此请使用末尾注释中可重现显示的 DF。

现在创建年份和值列（后者如果是 rollapply 输出）然后使用 read.zoo 将其转换为动物园系列，每年一列和月/日索引 0101、0102、...、1231。这将对齐每年的同一天创建专栏。然后取每一行的平均值，给出一年中每一天的所需平均值。 z 将是 366 x 3 -- 1 行代表一年中的每一天，3 列包括 1 列代表两年中的每一年加上平均值列。索引也存在，但存储为属性，而不是动物园对象中的列。 index(z) 可以用来查看。

如果您需要数据框，fortify.zoo(z) 会将 z 转换为数据框。

library(zoo)

z <- DF |>
  transform(year = as.integer(as.yearmon(date)),
            value = rollapply(value, 30, quantile, prob = 0.5, fill = NA)) |>
  read.zoo(split = "year", FUN = function(x) format(x, "%m%d"))
  transform(mean = rowMeans(na.rm = TRUE))

笔记

d <- seq(as.Date("2023-01-01"), as.Date("2024-12-31"), "day")
DF <- data.frame(date = d, value = seq_along(d))

【讨论】：