如何根据多年数据创建昼夜周期的 运行 中位数？
How can I create a running median of diel cycle from multiyear data?
我认为这个问题可能对处理长期环境变量的数据平滑的其他人感兴趣。
我的数据集结构如下：
列：
Date Hour_Min Y(response variable)
这些数据是每小时的，我需要创建昼夜周期的移动平均线，但按 Hour_Min 分类。换句话说，如果我要使用 31 天 window，对于给定的一天，Hour_Min 00:00 的 运行 平均数据点将取当天的平均值问题来自 Hour_Min 00:00 的前 15 天和后 15 天的数据点。然后，这将通过数据框在当天的小时 1:00 等重复。
不幸的是，数据也有很多 NA，这对于移动 window 平均值是有问题的，尽管我认为这可以使用 zoo 包中的 rollapply 来解决。
我尝试的一种方法是使用 tidyr 的传播功能从长格式切换到宽格式，以创建这样的数据帧：
Date Y_Hour_Min_0000 Y_Hour_Min_0100 Y_Hour_Min_0200 etc...
如果我能以这种方式更改格式，那么我就可以创建新的列，其中包含每个 Y_Hour_Min_.... 列的 运行 平均值。然后我需要将所有内容收集起来恢复为长格式（这是我不确定如何处理的另一项任务）。
但是，我无法使传播函数起作用，因此它将 Date 作为与每个 Y_Hour_Min_.... 列关联的分组变量。
另一个可能更优雅的解决方案是，如果有一种方法可以使用 rollapply 和自定义函数的某种组合，一步创建一个新列。
任何有关如何为此任务实现代码的想法都将不胜感激。下面我有一个简单的代码来模拟我的数据集：
模拟数据：
### Create vector of hours/dates:
date <- seq(as.POSIXct("2016-01-01 00:00"), as.POSIXct("2016-12-30
23:00"), by="hour")
### Create vector of noisy sine function:
d <- 365
n <- 24*d # number of data points
t <- seq(from = 0, to = 2*d*pi, length.out=24*d)
a <- 6
b <- 1
c.norm <- rnorm(n)
amp <- 3
y <- a*sin(b*t)+c.norm*amp+15
### Randomly insert NAs into data:
ind <- which(y %in% sample(y, 1000))
y[ind]<-NA
### Create test dataframe:
df <- data.frame(dt = date, y = y) %>%
separate(dt, c("date", "hour_min"), sep=" ") %>%
mutate(date = as.Date(date))
我会尝试一下，但它并不完美。希望有人能进来顶我。
TL:DR;
df2 <- df %>% slice(-7441) %>% spread(hour_min, y)
mov_avg <- function(x) {c(rep(NA, 15), rollapply(x, width = list(-15:15), FUN = mean, align="center", na.rm=T), rep(NA, 15))}
avgs <- as.data.frame(matrix(unlist(lapply(df2[,2:ncol(df2)], mov_avg)), nrow = nrow(df2), byrow = FALSE))
colnames(avgs) <- paste0("avg_", colnames(df2[,2:ncol(df2)]))
final_df <- cbind(df2, avgs) %>%
gather(2:ncol(.), key = "hour_min", value = "value") %>%
arrange(date, hour_min)
深度：
从你的起点开始。我添加了 set.seed(1) 这样我们就可以一起跟进了。
您的初始起点：
### Create vector of hours/dates:
set.seed(1)
date <- seq(as.POSIXct("2016-01-01 00:00"), as.POSIXct("2016-12-30
23:00"), by="hour")
### Create vector of noisy sine function:
d <- 365
n <- 24*d # number of data points
t <- seq(from = 0, to = 2*d*pi, length.out=24*d)
a <- 6
b <- 1
c.norm <- rnorm(n)
amp <- 3
y <- a*sin(b*t)+c.norm*amp+15
### Randomly insert NAs into data:
ind <- which(y %in% sample(y, 1000))
y[ind]<-NA
### Create test dataframe:
df <- data.frame(dt = date, y = y) %>%
separate(dt, c("date", "hour_min"), sep=" ") %>%
mutate(date = as.Date(date))
第一件事就是按你说的做，试试长格式。通常我认为这个问题最好通过在 hour_min 列上使用 dplyr 的 group_by 并在那里进行滚动平均，但我不确定该怎么做。
我注意到的第一件事是给定日期的一行有重复值。凌晨 1 点有两个观察结果，这打破了我们的 spread，所以我使用 slice(-7441)
删除了那个观察结果
所以让我们传播你的df。
df2 <- df %>% slice(-7441) %>% spread(hour_min, y)
正如我们所见，数据框现在有 365 个观察值长（日期），25 列宽（日期 + 24 小时）
dim(df2)
[1] 365 25
我做的下一件事是使用 rollapply，这是不完美的地方。使用 rollapply 时，我们可以给它一个 width = list(-15:15)。这将回顾过去的 15 天和未来的 15 天，并对所有 31 天进行平均。问题是前 15 天没有过去的 15，最后 15 天没有未来的 15。所以我用 NAs 填充了这些。我希望有人可以解决我的这部分答案。
我创建了一个自定义函数来执行此操作：
mov_avg <- function(x) {c(rep(NA, 15), rollapply(x, width = list(-15:15), FUN = mean, align="center", na.rm=T), rep(NA, 15))}
如果我们只做 rollapply，我们将得到一个长度为 335 的向量。我在前面和后面填充了 15 以得到我们需要的 365。
接下来我们希望 lapply 在整个数据帧中发挥作用。这将为我们提供一个包含 24 个长度为 365 的向量的列表。然后我们想将其转换为数据帧并将其绑定到我们当前的数据帧。
最后我们 gather 将所有列恢复为长格式并且 arrange
avgs <- as.data.frame(matrix(unlist(lapply(df2[,2:ncol(df2)], mov_avg)), nrow = nrow(df2), byrow = FALSE))
colnames(avgs) <- paste0("avg_", colnames(df2[,2:ncol(df2)]))
final_df <- cbind(df2, avgs) %>%
gather(2:ncol(.), key = "hour_min", value = "value") %>%
arrange(date, hour_min)
希望对您有所帮助。
我认为这可行：
编辑： 按照注释中的建议，通过将 fill = NA 参数添加到 rollapply() 函数来简化代码。
# add a complete date + time stamp
df$date_time <- paste(df$date, df$hour_min)
# make new column to store median data
df$median_y <- NA
# set rolling median width
width_roll <- 31
# do a rolling median for each hour, one at a time
# add NAs where no median can be calculated
for (i in levels(factor(df$hour_min))) {
df[df$hour_min == i, "median_y"] <- rollapply(df[df$hour_min == i, "y"],
width = width_roll,
median,
na.rm = TRUE,
fill = NA))
}
该方法只是按照您的建议使用 rollapply() 函数，但一次只能在一个特定的小时内使用。然后将它们中的每一个依次放回新的列中。
这里有一个全年只有一个小时的例子，这使得中值平滑更容易形象化。
# Examples:
# plot one hour plus rolling median over time
# here i = "23:00:00"
plot(x = as.POSIXct(df[df$hour_min == i, "date_time"]),
y = df[df$hour_min == i, "y"],
type = "l",
col = "blue",
ylab = "y values",
xlab = i)
lines(x = as.POSIXct(df[df$hour_min == i, "date_time"]),
y = df[df$hour_min == i, "median_y"],
lwd = 3)
legend("topleft",
legend = c("raw", "median"),
col = c("blue", "black"),
lwd = 3)
Plot for a single hour
这适用于所有内容（大量数据所以不太容易看到，但看起来很有效）。
# plot all the data
plot(x = as.POSIXct(df$date_time),
y = df$y,
type = "l",
col = "blue",
ylab = "y values",
xlab = "Date")
lines(x = as.POSIXct(df$date_time),
y = df$median_y,
lwd = 3)
legend("topleft",
legend = c("raw", "median"),
col = c("blue", "black"),
lwd = 3)
Plot for all data
我认为这个问题可能对处理长期环境变量的数据平滑的其他人感兴趣。
我的数据集结构如下：
列：
Date Hour_Min Y(response variable)
这些数据是每小时的，我需要创建昼夜周期的移动平均线，但按 Hour_Min 分类。换句话说，如果我要使用 31 天 window，对于给定的一天，Hour_Min 00:00 的 运行 平均数据点将取当天的平均值问题来自 Hour_Min 00:00 的前 15 天和后 15 天的数据点。然后，这将通过数据框在当天的小时 1:00 等重复。
不幸的是，数据也有很多 NA，这对于移动 window 平均值是有问题的，尽管我认为这可以使用 zoo 包中的 rollapply 来解决。
我尝试的一种方法是使用 tidyr 的传播功能从长格式切换到宽格式，以创建这样的数据帧：
Date Y_Hour_Min_0000 Y_Hour_Min_0100 Y_Hour_Min_0200 etc...
如果我能以这种方式更改格式，那么我就可以创建新的列，其中包含每个 Y_Hour_Min_.... 列的 运行 平均值。然后我需要将所有内容收集起来恢复为长格式（这是我不确定如何处理的另一项任务）。
但是，我无法使传播函数起作用，因此它将 Date 作为与每个 Y_Hour_Min_.... 列关联的分组变量。
另一个可能更优雅的解决方案是，如果有一种方法可以使用 rollapply 和自定义函数的某种组合，一步创建一个新列。
任何有关如何为此任务实现代码的想法都将不胜感激。下面我有一个简单的代码来模拟我的数据集：
模拟数据：
### Create vector of hours/dates:
date <- seq(as.POSIXct("2016-01-01 00:00"), as.POSIXct("2016-12-30
23:00"), by="hour")
### Create vector of noisy sine function:
d <- 365
n <- 24*d # number of data points
t <- seq(from = 0, to = 2*d*pi, length.out=24*d)
a <- 6
b <- 1
c.norm <- rnorm(n)
amp <- 3
y <- a*sin(b*t)+c.norm*amp+15
### Randomly insert NAs into data:
ind <- which(y %in% sample(y, 1000))
y[ind]<-NA
### Create test dataframe:
df <- data.frame(dt = date, y = y) %>%
separate(dt, c("date", "hour_min"), sep=" ") %>%
mutate(date = as.Date(date))
我会尝试一下，但它并不完美。希望有人能进来顶我。
TL:DR;
df2 <- df %>% slice(-7441) %>% spread(hour_min, y)
mov_avg <- function(x) {c(rep(NA, 15), rollapply(x, width = list(-15:15), FUN = mean, align="center", na.rm=T), rep(NA, 15))}
avgs <- as.data.frame(matrix(unlist(lapply(df2[,2:ncol(df2)], mov_avg)), nrow = nrow(df2), byrow = FALSE))
colnames(avgs) <- paste0("avg_", colnames(df2[,2:ncol(df2)]))
final_df <- cbind(df2, avgs) %>%
gather(2:ncol(.), key = "hour_min", value = "value") %>%
arrange(date, hour_min)
深度：
从你的起点开始。我添加了 set.seed(1) 这样我们就可以一起跟进了。
您的初始起点：
### Create vector of hours/dates:
set.seed(1)
date <- seq(as.POSIXct("2016-01-01 00:00"), as.POSIXct("2016-12-30
23:00"), by="hour")
### Create vector of noisy sine function:
d <- 365
n <- 24*d # number of data points
t <- seq(from = 0, to = 2*d*pi, length.out=24*d)
a <- 6
b <- 1
c.norm <- rnorm(n)
amp <- 3
y <- a*sin(b*t)+c.norm*amp+15
### Randomly insert NAs into data:
ind <- which(y %in% sample(y, 1000))
y[ind]<-NA
### Create test dataframe:
df <- data.frame(dt = date, y = y) %>%
separate(dt, c("date", "hour_min"), sep=" ") %>%
mutate(date = as.Date(date))
第一件事就是按你说的做，试试长格式。通常我认为这个问题最好通过在 hour_min 列上使用 dplyr 的 group_by 并在那里进行滚动平均，但我不确定该怎么做。
我注意到的第一件事是给定日期的一行有重复值。凌晨 1 点有两个观察结果，这打破了我们的 spread，所以我使用 slice(-7441)
所以让我们传播你的df。
df2 <- df %>% slice(-7441) %>% spread(hour_min, y)
正如我们所见，数据框现在有 365 个观察值长（日期），25 列宽（日期 + 24 小时）
dim(df2)
[1] 365 25
我做的下一件事是使用 rollapply，这是不完美的地方。使用 rollapply 时，我们可以给它一个 width = list(-15:15)。这将回顾过去的 15 天和未来的 15 天，并对所有 31 天进行平均。问题是前 15 天没有过去的 15，最后 15 天没有未来的 15。所以我用 NAs 填充了这些。我希望有人可以解决我的这部分答案。
我创建了一个自定义函数来执行此操作：
mov_avg <- function(x) {c(rep(NA, 15), rollapply(x, width = list(-15:15), FUN = mean, align="center", na.rm=T), rep(NA, 15))}
如果我们只做 rollapply，我们将得到一个长度为 335 的向量。我在前面和后面填充了 15 以得到我们需要的 365。
接下来我们希望 lapply 在整个数据帧中发挥作用。这将为我们提供一个包含 24 个长度为 365 的向量的列表。然后我们想将其转换为数据帧并将其绑定到我们当前的数据帧。
最后我们 gather 将所有列恢复为长格式并且 arrange
avgs <- as.data.frame(matrix(unlist(lapply(df2[,2:ncol(df2)], mov_avg)), nrow = nrow(df2), byrow = FALSE))
colnames(avgs) <- paste0("avg_", colnames(df2[,2:ncol(df2)]))
final_df <- cbind(df2, avgs) %>%
gather(2:ncol(.), key = "hour_min", value = "value") %>%
arrange(date, hour_min)
希望对您有所帮助。
我认为这可行：
编辑： 按照注释中的建议，通过将 fill = NA 参数添加到 rollapply() 函数来简化代码。
# add a complete date + time stamp
df$date_time <- paste(df$date, df$hour_min)
# make new column to store median data
df$median_y <- NA
# set rolling median width
width_roll <- 31
# do a rolling median for each hour, one at a time
# add NAs where no median can be calculated
for (i in levels(factor(df$hour_min))) {
df[df$hour_min == i, "median_y"] <- rollapply(df[df$hour_min == i, "y"],
width = width_roll,
median,
na.rm = TRUE,
fill = NA))
}
该方法只是按照您的建议使用 rollapply() 函数，但一次只能在一个特定的小时内使用。然后将它们中的每一个依次放回新的列中。
这里有一个全年只有一个小时的例子，这使得中值平滑更容易形象化。
# Examples:
# plot one hour plus rolling median over time
# here i = "23:00:00"
plot(x = as.POSIXct(df[df$hour_min == i, "date_time"]),
y = df[df$hour_min == i, "y"],
type = "l",
col = "blue",
ylab = "y values",
xlab = i)
lines(x = as.POSIXct(df[df$hour_min == i, "date_time"]),
y = df[df$hour_min == i, "median_y"],
lwd = 3)
legend("topleft",
legend = c("raw", "median"),
col = c("blue", "black"),
lwd = 3)
Plot for a single hour
这适用于所有内容（大量数据所以不太容易看到，但看起来很有效）。
# plot all the data
plot(x = as.POSIXct(df$date_time),
y = df$y,
type = "l",
col = "blue",
ylab = "y values",
xlab = "Date")
lines(x = as.POSIXct(df$date_time),
y = df$median_y,
lwd = 3)
legend("topleft",
legend = c("raw", "median"),
col = c("blue", "black"),
lwd = 3)
Plot for all data