7

我希望在应用模型之前删除异常值。我正在使用黄土曲线来划定趋势线并设置异常值限制。我想删除超出定义限制的行。除了使用自定义函数执行此操作外，该函数一次获取每个点并检查当地的黄土坡度等……还有更简单的方法吗？

有界限的黄土趋势线 (1.2)

# Code generating image above
scatter.smooth( idam$T_d, idam$T_x10d)
loessline <- loess.smooth( idam$T_d, idam$T_x10d)
lines(loessline$x, loessline$y, lwd=3)
lines(loessline$x, loessline$y*1.2, lwd=3, col='red')
lines(loessline$x, loessline$y/1.2, lwd=3, col='red')
4

3 回答 3

8

您可以使用approxfun

这是一个带有“异常值”的示例

plot(wt ~ mpg, data = mtcars)
lo <- loess.smooth(mtcars$mpg, mtcars$wt)
lines(lo$x, lo$y, lwd = 3)
lines(lo$x, lo$y * 1.2, lwd = 3, col = 2)
lines(lo$x, lo$y / 1.2, lwd = 3, col = 2)

在此处输入图像描述

approxfun返回一个使用观察到的 x 值的函数，我们可以使用它来插入一组新的 y 值。

然后，您可以设置将点称为异常值的阈值；在这里，我使用1.2 * y原始问题中的方法来识别极端观察结果。

f1 <- approxfun(lo$x, lo$y * 1.2)
(wh1 <- which(mtcars$wt > f1(mtcars$mpg)))
# [1]  8 17 18

f2 <- approxfun(lo$x, lo$y / 1.2)
(wh2 <- which(mtcars$wt < f2(mtcars$mpg)))
# [1] 28

## identify points to exclude
mt <- mtcars[c(wh1, wh2), ]
points(mt$mpg, mt$wt, pch = 4, col = 2, cex = 2)

在此处输入图像描述

## plot without points
plot(wt ~ mpg, data = mt2 <- mtcars[-c(wh1, wh2), ])
lo <- loess.smooth(mt2$mpg, mt2$wt)
lines(lo$x, lo$y, lwd = 3)
lines(lo$x, lo$y * 1.2, lwd = 3, col = 2)
lines(lo$x, lo$y / 1.2, lwd = 3, col = 2)

在此处输入图像描述

由于这里有几个步骤，您可以将其打包成一个函数以使事情变得更容易：

par(mfrow = c(2,2))
with(mtcars, {
  plot_lo(mpg, wt)
  plot_lo(mpg, wt, limits = c(1 / 1.5, 1.5))
  dd <<- plot_lo(mpg, wt, limits = c(1 / 1.2, 1.2))
  plot_lo(mpg, wt, pch = 16, las = 1, tcl = .5, bty = 'l')
})

str(dd)
# List of 2
# $ x: num [1:28] 21 21 22.8 21.4 18.7 18.1 14.3 22.8 19.2 17.8 ...
# $ y: num [1:28] 2.62 2.88 2.32 3.21 3.44 ...

在此处输入图像描述

plot_lo <- function(x, y, limits = c(-Inf, Inf), ...) {
  lo <- loess.smooth(x, y)
  fx <- approxfun(lo$x, lo$y * limits[1L])
  fy <- approxfun(lo$x, lo$y * limits[2L])

  idx <- which(y < fx(x) | y > fy(x))
  if (length(idx)) {
    x  <- x[-idx]
    y  <- y[-idx]
    lo <- loess.smooth(x, y)
  }

  op <- par(..., no.readonly = TRUE)
  on.exit(par(op))

  plot(x, y)
  lines(lo$x, lo$y, lwd = 3)
  lines(lo$x, lo$y * limits[1L], lwd = 3, col = 2L)
  lines(lo$x, lo$y * limits[2L], lwd = 3, col = 2L)

  invisible(list(x = x, y = y))
}
于 2016-01-19T16:00:49.707 回答
8

检测异常值可以在 DBSCAN R 包的帮助下完成，DBSCAN R 包是用于集群识别的著名算法（有关更多详细信息，请参阅WIKIPEDIA）。

这个函数有三个重要的输入：

  • x：您的数据（仅数值）
  • eps：目标最大距离
  • minPts：将它们视为集群的最小点数

评估 eps 可以在 knndist(...) 和 knndistplot(...) 函数的帮助下完成：

  • knndistplot 将为给定的 k（即 minPts）绘制数据集上的 eps 值 ==> 您可以直观地选择一个有效的 eps 值（通常在膝盖曲线部分）
  • knndist 将评估 eps 值并将它们从矩阵中返回。k 输入将生成 1:1:k 估值，您可以使用结果以编程方式确定准确的 eps 和 k 值

接下来，您只需使用 dbscan(yourdata, eps, k) 来获取具有以下组件的 dbscan 对象：

  • eps：用于计算的 eps
  • minPts：识别一个簇的最小点数
  • cluster：一个整数向量，用于标识属于（=1）或不属于（=0）的点。最后一个对应于您要消除的异常值。

请注意 dbscan 的以下限制：

  • dbscan 使用欧几里得距离，因此它被提交给“维度诅咒”。这可以通过使用 PCA 来避免
  • dbscan 消除了可能产生未知点的叠加点。这可以通过使用左外连接将结果与数据合并或使用 jitter(...) 函数来解决，该函数会为数据添加噪声。根据您显示的数据，我认为您的数据可能是这种情况

知道这个限制，dbscan 包提供了两种替代方法：LOF 和 OPTICS（DBSCAN 的扩展）

编辑于 2016 年 1 月 25 日

在@rawr 回答之后，我将给出一个基于mtcars数据集的示例，以展示如何dbscan用于识别异常值。请注意，我的示例将使用优秀的data.table包而不是经典的data.frame.

首先，我开始复制 rawr 的方法来说明 data.table 的使用

require(data.table)
require(ggplot2)
require(dbscan)
data(mtcars)
dt_mtcars <- as.data.table(mtcars)

# based on rawr's approach
plot(wt~mpg, data=dt_mtcars)
lo <- loess.smooth(dt_mtcars[,mpg], dt_mtcars[,wt])
lines(lo$x,lo$y, lwd=3)
lines(lo$x,lo$y * 1.2, lwd=3 , col=2 )
lines(lo$x,lo$y / 1.2, lwd=3 , col=2 )

在此处输入图像描述

因此，我们可以评估我们是否获得了独立于基础支持的相同结果。

其次，以下代码说明了 DBSCAN 方法，该方法从确定eps和开始k，即识别集群的必要点数：

res_knn = kNNdist( dt_mtcars[, .(wt, mpg)] , k = 10)
dim_knn = dim(res_knn)
x_knn =  seq(1, dim_knn[1])
ggplot() + 
   geom_line( aes( x = x_knn , y = sort(res_knn[, 1])  , col = 1 ) ) + 
   geom_line( aes( x = x_knn , y = sort(res_knn[, 2])  , col = 2 ) ) + 
   geom_line( aes( x = x_knn , y = sort(res_knn[, 3])  , col = 3 ) ) + 
   geom_line( aes( x = x_knn , y = sort(res_knn[, 4])  , col = 4 ) ) + 
   geom_line( aes( x = x_knn , y = sort(res_knn[, 5])  , col = 5 ) ) + 
   geom_line( aes( x = x_knn , y = sort(res_knn[, 6])  , col = 6 ) ) + 
   geom_line( aes( x = x_knn , y = sort(res_knn[, 7])  , col = 7 ) ) + 
   geom_line( aes( x = x_knn , y = sort(res_knn[, 8])  , col = 8 ) ) + 
   geom_line( aes( x = x_knn , y = sort(res_knn[, 9])  , col = 9 ) ) + 
   geom_line( aes( x = x_knn , y = sort(res_knn[, 10]) , col = 10 ) )  +
   xlab('sorted results') + 
   ylab('kNN distance')

结果绘制在下图中：

在此处输入图像描述

它表明计算出的 kNN 距离对该因子敏感，k但是分离异常值的准确eps值位于曲线的膝部 ==> 合适eps的位于 2 到 4 之间。这是一种视觉评估，可以通过适当的自动化搜索算法（例如，请参阅此链接）。关于k，必须定义权衡，知道 k 越低，结果越不严格。

eps = 3在下一部分中，我们将使用（基于视觉估计）对dbscan 进行参数化，并k = 4获得稍微严格的结果。我们将借助 rawr 的代码绘制这些结果：

eps = 3
k = 4
res_dbscan = dbscan( dt_mtcars[, .(wt, mpg)] , eps , k )
plot(wt~mpg, data=dt_mtcars, col = res_dbscan$cluster)
lo <- loess.smooth(dt_mtcars[res_dbscan$cluster>0,mpg], dt_mtcars[res_dbscan$cluster>0,wt])
lines(lo$x,lo$y, lwd=3)
lines(lo$x,lo$y * 1.2, lwd=3 , col=2 )
lines(lo$x,lo$y / 1.2, lwd=3 , col=2 )

在此处输入图像描述

我们得到这个数字，我们可以评估我们从 rawr 的方法得到不同的结果，其中位于的点mpg = [10,13]被视为异常值。

与 rawr 的解决方案相比，这些结果可能被认为是奇怪的，后者在具有双变量数据 (Y ~ X) 的假设下工作。然而mtcars，变量之间的关系可能是（或不是）线性的多维数据集......为了评估这一点，我们可以散点图这个数据集，例如过滤数值

pairs(dt_mtcars[, .(mpg, disp, hp, drat, wt, qsec)])

在此处输入图像描述

如果我们只关注结果wt ~ mpg，我们乍一看可能会认为它是一种反线性关系。但是对于其他绘制的关系，情况可能并非如此，并且在 N-Dim 环境中查找异常值有点棘手。实际上，当在特定的 2D 比较中投影时，一个点可能会被视为异常值……但如果我们添加一个新的比较维度，则相反。实际上，我们可能具有可以识别的共线性，从而加强或不加强集群关系。

我的朋友们，我同意它很多，if为了说明这种情况，我们将继续对 的dbscan数值进行分析mtcars。

所以我将复制前面介绍的过程，让我们从 kNN 距离分析开始：

res_knn = kNNdist( dt_mtcars[, .(mpg, disp, hp, drat, wt, qsec)] , k = 10)
dim_knn = dim(res_knn)
x_knn =  seq(1, dim_knn[1])
ggplot() + 
    geom_line( aes( x = x_knn , y = sort(res_knn[, 1])  , col = 1 ) ) + 
    geom_line( aes( x = x_knn , y = sort(res_knn[, 2])  , col = 2 ) ) + 
    geom_line( aes( x = x_knn , y = sort(res_knn[, 3])  , col = 3 ) ) + 
    geom_line( aes( x = x_knn , y = sort(res_knn[, 4])  , col = 4 ) ) + 
    geom_line( aes( x = x_knn , y = sort(res_knn[, 5])  , col = 5 ) ) + 
    geom_line( aes( x = x_knn , y = sort(res_knn[, 6])  , col = 6 ) ) + 
    geom_line( aes( x = x_knn , y = sort(res_knn[, 7])  , col = 7 ) ) + 
    geom_line( aes( x = x_knn , y = sort(res_knn[, 8])  , col = 8 ) ) + 
    geom_line( aes( x = x_knn , y = sort(res_knn[, 9])  , col = 9 ) ) + 
    geom_line( aes( x = x_knn , y = sort(res_knn[, 10]) , col = 10 ) )  +
    xlab('sorted results') + 
    ylab('kNN distance')

排序的 kNN 距离

与 上产生的分析相比wt ~ mpg，我们可以看到kNNdist(...)产生了更重要的 kNN 距离（k = 10例如，直到 200）。然而，我们仍然有膝盖部分，它可以帮助我们估计一个合适的eps值。

在下一部分中，我们将使用eps = 75andk = 5和

# optimal eps value is between 40 (k=1) and 130 (k=10)
eps = 75
k = 5
res_dbscan = dbscan( dt_mtcars[, .(mpg, disp, hp, drat, wt, qsec)] , eps , k )
pairs(dt_mtcars[, .(mpg, disp, hp, drat, wt, qsec)] , col = res_dbscan$cluster+2L)

在此处输入图像描述

因此，该分析的散点图突出表明，由于变量之间的复杂关系，在 N-Dim 环境中识别异常值可能很棘手。但请注意，在大多数情况下，异常值位于 2D 投影的角落部分，这加强了用wt ~ mpg

于 2016-01-19T16:25:08.650 回答
0

我的建议是去看看outliers package。该软件包允许在分析发生之前进行识别。这是一个非常简单的例子：

library(outliers)
series<-c(runif(100,1,2),1000)
round(scores(series,prob=1,type="chisq"),3)

使用此功能，可以执行大量测试，并且您可以设置您感到满意的异常值的概率水平。

series<-series[which(series<0.95),]
于 2016-01-22T18:38:16.613 回答