R纠正不一致的数据记录
R correcting inconsistent data records
我每 5 分钟记录四个变量。当我在 R 中绘制四个变量的时间序列时，我意识到由于数据收集错误（记录设备 error/sensor 错误），变量 3 正在记录不一致的数据。如何更正数据记录？
变量3的数据记录有一些异常跳跃，这不是研究变量的物理效应。该图显示了一周的数据记录以及每天的振荡。连续两个读数之间不应该有如此高的跳跃。
前段时间我尝试了一些 R 离群值包，但没有得到结果...
当我绘制整个时间序列时，结果更糟。
任何帮助将不胜感激。
谢谢
我分享图片1的数据：
这是一个可能的解决方案，但首先我们需要生成一些代表您的问题的数据。您的方案的好处在于，虚假数据点是大尖峰，即使在视觉上也相当明显。
生成数据
set.seed(15161)
x <- seq(pi/10,10*pi,by=pi/100)
y <- sin(x) # using sin() generates some osciliating data
z <- sample(c(0,-5),length(y),
prob=c(0.99,0.01),replace=TRUE) # pepper the data with random spikes
y <- y + z
df <- data.frame(cbind(x,y,z))
length(which(df$z==-5)) # the number of spikes ~ 13
plot(df$x,df$y,type="l",ylim=c(-10,2),col="blue",xlab="x",ylab="y")
abline(h=0,lty=5)
删除虚假测量（清理数据）
在您提供的数据中，与良好测量的背景相比，虚假数据点非常大。那就是您的测量值以缓慢增加或减少的方式很好地移动，然后将 jump/drop 归咎于 > 20 个单位。所以我写了一个函数，它将发现并删除任何代表 increase/decrease 高于某个阈值的数据点（在你的例子中~20 个单位，在我上面的工作示例中~2 个单位就足够了）。
函数代码为：
f <- function(df,clean,threshold){
y <- df[,clean]
for(i in 1:length(y)){
if(is.na(y[i]) | is.na(y[i+1])){
next
}
if(abs(y[i+1]-y[i])>threshold){
y[i+1] <- NA
}
}
return(df[!is.na(y),])
}
cleaned.df <- f(df,clean="y",threshold=2) # Run the function to clean the data
length(which(cleaned.df$z==-5)) # number of spikes in cleaned data is now 0
绘制清理后的结果
plot(cleaned.df$x,cleaned.df$y,type="l",ylim=c(-10,2),col="blue",xlab="x",ylab="y")
abline(h=0,lty=5)
注意事项和注意事项
- 在运行启用函数（即按时间顺序排序的测量值）
之前，请确保您的数据按顺序
- 我建议您选择 20 个单位左右的阈值（仅通过目视检查您的图形，这似乎足够了。
- 清洁功能可能无法有效去除 2 个或更多连续的尖峰。 但是您可以运行多次通过清理功能来处理数据，这应该有效。
- 我们可以设计更严格的方法，但我认为此解决方案简单有效。如果您仍有问题，请告诉我们，我们可以制定更严格的解决方案。
编辑 1：
我刚刚看到您上传了一些实际数据。稍微调整函数以适应改变测量符号的尖峰。以下是适用于您的数据的结果，看起来对我有用。
df <- read.csv("figure1data.csv")
plot(df$X,df$three,type="l",col="blue",xlab="x",ylab="y",ylim=c(-150,50))
abline(h=0,lty=5)
cleaned.df1 <- f(df,clean="three",threshold=20)
plot(cleaned.df1$X,cleaned.df1$three,type="l",col="blue",xlab="x",ylab="y",
ylim=c(-150,50))
abline(h=0,lty=5)
编辑 2：对 OP 评论的回应
要消除连续出现峰值的情况，只需重新运行 清理数据上的函数即可。
cleaned.df2 <- f(cleaned.df1,clean="three",threshold=20)
要将所有行恢复到数据并将尖峰变量 "three" 点转换为 NA，只需按如下方式合并数据即可。
New.df <- merge(df[,colnames(df)!="three"],
cleaned.df2[,colnames(df) %in% c("X","three")],
by="X",all.x=TRUE)
检查是否按预期工作
df[which(!complete.cases(New.df)),]
New.df[which(!complete.cases(New.df)),]
您清楚地看到具有变量 "three" 尖峰的行现在位于 New.df
中的 NA
读完你的数据并绘图后，我看到了这个：
df <- read.csv("~/Whosebug/RaülOo.csv")
df$TIMESTAMP <- as.POSIXct(df$TIMESTAMP)
library(dplyr)
library(tidyr)
library(ggplot2)
gather(df, k, v, -X, -TIMESTAMP) %>%
ggplot(aes(TIMESTAMP, v, color=k)) +
geom_path()
是否像"anything above -50"一样简单？十分位数如下所示：
quantile(unlist(df[,3:6]), seq(0,1,len=11))
# 0% 10% 20% 30% 40% 50% 60% 70%
# -122.7000 -22.9600 -17.5500 -13.4200 -10.0700 -5.9615 3.4800 16.0500
# 80% 90% 100%
# 26.6040 35.6860 81.4000
IQR 约为 37。与箱线图中的 "whiskers" 类似，假设 "1.5 IQR" 可能是现实的，即：值低于“下四分位数以下 IQR 的 1.5 倍”（同样以上，虽然不存在于此数据中）可能被安全地视为异常值。
(q <- quantile(unlist(df[,3:6]), c(0.25, 0.75)))
# 25% 75%
# -15.4000 22.0025
unname( q[1] - 1.5*diff(q) ) # "unname" only to remove the now-misleading percentile label
# -71.50375
gather(df, k, v, -X, -TIMESTAMP) %>%
filter(v > q[1] - 1.5*diff(q)) %>%
ggplot(aes(TIMESTAMP, v, color=k)) +
geom_path()
所以 1.5 可能不足以真正识别异常值，但这取决于您的需要。如果您只需要一个清理过的图（并且一些异常值不会使人衰弱），那么我建议使用标准的“1.5 倍 IQR”就足够了。如果你想更好地控制它，也许使用更接近 1 的东西会起作用。
gather(df, k, v, -X, -TIMESTAMP) %>%
filter(v > q[1] - diff(q)) %>%
ggplot(aes(TIMESTAMP, v, color=k)) +
geom_path()
如果您需要以 "wide" 格式返回，您可以这样做：
gather(df, k, v, -X, -TIMESTAMP) %>%
filter(v > -50) %>%
spread(k, v) %>%
slice(37:43) # just for demonstration
# X TIMESTAMP four one three two
# 1 37 2018-07-15 03:05:00 -21.68 -32.04 -23.11 -12.87
# 2 38 2018-07-15 03:10:00 -21.79 -31.71 -23.11 -12.87
# 3 39 2018-07-15 03:15:00 -21.79 -31.71 -23.11 -12.87
# 4 40 2018-07-15 03:20:00 -21.79 -31.71 -23.11 -12.87
# 5 41 2018-07-15 03:25:00 -17.43 -25.37 NA -10.29
# 6 42 2018-07-15 03:30:00 -21.79 -31.71 -23.11 -12.87
# 7 43 2018-07-15 03:35:00 -21.79 -31.28 -23.11 -12.87
你的异常值现在在哪里 NA。更简洁的非 dplyr/tidyr 替代方案可能是：
df[,3:6] <- lapply(df[,3:6], function(a) ifelse(a < -50, NA, a))
然后无论您进行任何后续处理或绘图，都需要考虑（忽略）NA 值。
我会更进一步，因为（对您来说）了解坏数据进入的频率（或周期性）可能很有趣。
newdat <- df %>%
gather(k, v, -X, -TIMESTAMP) %>%
mutate(v = if_else(v < q[1] - diff(q), NA_real_, v))
baddat <- filter(newdat, is.na(v))
newdat <- filter(newdat, !is.na(v))
baddat$v <- min(newdat$v) - 5 # arbitrary
ggplot(newdat, aes(TIMESTAMP, v, color = k)) +
geom_path() +
geom_point(data = baddat)
在这里您可以看到问题数据点所在的位置，而无需扩展图表的其余部分。
备注
这似乎是一个快速入门的技巧。例如，如果您的四个不同测量值不是同质的，而是在截然不同的尺度上，则需要按列完成。
我使用 dplyr 进行数据处理，尽管它们不是严格要求的。这可以很容易地在 base-R 中用相对简单的函数完成。使用 ggplot2 强制长数据，因此 tidyr::gather（和 tidyr::spread）；如果您使用的是基础图形，那么您可能不需要重塑数据（这表明可能首选每列数据替换）。
我每 5 分钟记录四个变量。当我在 R 中绘制四个变量的时间序列时，我意识到由于数据收集错误（记录设备 error/sensor 错误），变量 3 正在记录不一致的数据。如何更正数据记录？
变量3的数据记录有一些异常跳跃，这不是研究变量的物理效应。该图显示了一周的数据记录以及每天的振荡。连续两个读数之间不应该有如此高的跳跃。 前段时间我尝试了一些 R 离群值包，但没有得到结果...
当我绘制整个时间序列时，结果更糟。
任何帮助将不胜感激。 谢谢
我分享图片1的数据：
这是一个可能的解决方案，但首先我们需要生成一些代表您的问题的数据。您的方案的好处在于，虚假数据点是大尖峰，即使在视觉上也相当明显。
生成数据
set.seed(15161)
x <- seq(pi/10,10*pi,by=pi/100)
y <- sin(x) # using sin() generates some osciliating data
z <- sample(c(0,-5),length(y),
prob=c(0.99,0.01),replace=TRUE) # pepper the data with random spikes
y <- y + z
df <- data.frame(cbind(x,y,z))
length(which(df$z==-5)) # the number of spikes ~ 13
plot(df$x,df$y,type="l",ylim=c(-10,2),col="blue",xlab="x",ylab="y")
abline(h=0,lty=5)
删除虚假测量（清理数据）
在您提供的数据中，与良好测量的背景相比，虚假数据点非常大。那就是您的测量值以缓慢增加或减少的方式很好地移动，然后将 jump/drop 归咎于 > 20 个单位。所以我写了一个函数，它将发现并删除任何代表 increase/decrease 高于某个阈值的数据点（在你的例子中~20 个单位，在我上面的工作示例中~2 个单位就足够了）。
函数代码为：
f <- function(df,clean,threshold){
y <- df[,clean]
for(i in 1:length(y)){
if(is.na(y[i]) | is.na(y[i+1])){
next
}
if(abs(y[i+1]-y[i])>threshold){
y[i+1] <- NA
}
}
return(df[!is.na(y),])
}
cleaned.df <- f(df,clean="y",threshold=2) # Run the function to clean the data
length(which(cleaned.df$z==-5)) # number of spikes in cleaned data is now 0
绘制清理后的结果
plot(cleaned.df$x,cleaned.df$y,type="l",ylim=c(-10,2),col="blue",xlab="x",ylab="y")
abline(h=0,lty=5)
注意事项和注意事项
- 在运行启用函数（即按时间顺序排序的测量值） 之前，请确保您的数据按顺序
- 我建议您选择 20 个单位左右的阈值（仅通过目视检查您的图形，这似乎足够了。
- 清洁功能可能无法有效去除 2 个或更多连续的尖峰。 但是您可以运行多次通过清理功能来处理数据，这应该有效。
- 我们可以设计更严格的方法，但我认为此解决方案简单有效。如果您仍有问题，请告诉我们，我们可以制定更严格的解决方案。
编辑 1：
我刚刚看到您上传了一些实际数据。稍微调整函数以适应改变测量符号的尖峰。以下是适用于您的数据的结果，看起来对我有用。
df <- read.csv("figure1data.csv")
plot(df$X,df$three,type="l",col="blue",xlab="x",ylab="y",ylim=c(-150,50))
abline(h=0,lty=5)
cleaned.df1 <- f(df,clean="three",threshold=20)
plot(cleaned.df1$X,cleaned.df1$three,type="l",col="blue",xlab="x",ylab="y",
ylim=c(-150,50))
abline(h=0,lty=5)
编辑 2：对 OP 评论的回应
要消除连续出现峰值的情况，只需重新运行 清理数据上的函数即可。
cleaned.df2 <- f(cleaned.df1,clean="three",threshold=20)
要将所有行恢复到数据并将尖峰变量 "three" 点转换为 NA，只需按如下方式合并数据即可。
New.df <- merge(df[,colnames(df)!="three"],
cleaned.df2[,colnames(df) %in% c("X","three")],
by="X",all.x=TRUE)
检查是否按预期工作
df[which(!complete.cases(New.df)),]
New.df[which(!complete.cases(New.df)),]
您清楚地看到具有变量 "three" 尖峰的行现在位于 New.df
NA
读完你的数据并绘图后，我看到了这个：
df <- read.csv("~/Whosebug/RaülOo.csv")
df$TIMESTAMP <- as.POSIXct(df$TIMESTAMP)
library(dplyr)
library(tidyr)
library(ggplot2)
gather(df, k, v, -X, -TIMESTAMP) %>%
ggplot(aes(TIMESTAMP, v, color=k)) +
geom_path()
是否像"anything above -50"一样简单？十分位数如下所示：
quantile(unlist(df[,3:6]), seq(0,1,len=11))
# 0% 10% 20% 30% 40% 50% 60% 70%
# -122.7000 -22.9600 -17.5500 -13.4200 -10.0700 -5.9615 3.4800 16.0500
# 80% 90% 100%
# 26.6040 35.6860 81.4000
IQR 约为 37。与箱线图中的 "whiskers" 类似，假设 "1.5 IQR" 可能是现实的，即：值低于“下四分位数以下 IQR 的 1.5 倍”（同样以上，虽然不存在于此数据中）可能被安全地视为异常值。
(q <- quantile(unlist(df[,3:6]), c(0.25, 0.75)))
# 25% 75%
# -15.4000 22.0025
unname( q[1] - 1.5*diff(q) ) # "unname" only to remove the now-misleading percentile label
# -71.50375
gather(df, k, v, -X, -TIMESTAMP) %>%
filter(v > q[1] - 1.5*diff(q)) %>%
ggplot(aes(TIMESTAMP, v, color=k)) +
geom_path()
所以 1.5 可能不足以真正识别异常值，但这取决于您的需要。如果您只需要一个清理过的图（并且一些异常值不会使人衰弱），那么我建议使用标准的“1.5 倍 IQR”就足够了。如果你想更好地控制它，也许使用更接近 1 的东西会起作用。
gather(df, k, v, -X, -TIMESTAMP) %>%
filter(v > q[1] - diff(q)) %>%
ggplot(aes(TIMESTAMP, v, color=k)) +
geom_path()
如果您需要以 "wide" 格式返回，您可以这样做：
gather(df, k, v, -X, -TIMESTAMP) %>%
filter(v > -50) %>%
spread(k, v) %>%
slice(37:43) # just for demonstration
# X TIMESTAMP four one three two
# 1 37 2018-07-15 03:05:00 -21.68 -32.04 -23.11 -12.87
# 2 38 2018-07-15 03:10:00 -21.79 -31.71 -23.11 -12.87
# 3 39 2018-07-15 03:15:00 -21.79 -31.71 -23.11 -12.87
# 4 40 2018-07-15 03:20:00 -21.79 -31.71 -23.11 -12.87
# 5 41 2018-07-15 03:25:00 -17.43 -25.37 NA -10.29
# 6 42 2018-07-15 03:30:00 -21.79 -31.71 -23.11 -12.87
# 7 43 2018-07-15 03:35:00 -21.79 -31.28 -23.11 -12.87
你的异常值现在在哪里 NA。更简洁的非 dplyr/tidyr 替代方案可能是：
df[,3:6] <- lapply(df[,3:6], function(a) ifelse(a < -50, NA, a))
然后无论您进行任何后续处理或绘图，都需要考虑（忽略）NA 值。
我会更进一步，因为（对您来说）了解坏数据进入的频率（或周期性）可能很有趣。
newdat <- df %>%
gather(k, v, -X, -TIMESTAMP) %>%
mutate(v = if_else(v < q[1] - diff(q), NA_real_, v))
baddat <- filter(newdat, is.na(v))
newdat <- filter(newdat, !is.na(v))
baddat$v <- min(newdat$v) - 5 # arbitrary
ggplot(newdat, aes(TIMESTAMP, v, color = k)) +
geom_path() +
geom_point(data = baddat)
在这里您可以看到问题数据点所在的位置，而无需扩展图表的其余部分。
备注
这似乎是一个快速入门的技巧。例如，如果您的四个不同测量值不是同质的，而是在截然不同的尺度上，则需要按列完成。
我使用
dplyr进行数据处理，尽管它们不是严格要求的。这可以很容易地在 base-R 中用相对简单的函数完成。使用ggplot2强制长数据，因此tidyr::gather（和tidyr::spread）；如果您使用的是基础图形，那么您可能不需要重塑数据（这表明可能首选每列数据替换）。