从相关系数计算中去除异常值

Posted

技术标签:

【中文标题】从相关系数计算中去除异常值【英文标题】:Remove outliers from correlation coefficient calculation 【发布时间】:2011-06-07 16:26:57 【问题描述】:

假设我们有两个数字向量 xyxy 之间的 Pearson 相关系数由下式给出

cor(x, y)

如何在计算中自动考虑 xy 的子集(比如 90%)以最大化相关系数?

【问题讨论】:

您认为这里的异常值是什么?偏离最小二乘拟合线(即最大残差),或xy 二元分布的极端值? @Gavin 这里我认为最大的残差是异常值。 【参考方案1】:

您可以尝试引导数据以找到最高的相关系数,例如:

x <- cars$dist
y <- cars$speed
percent <- 0.9         # given in the question above
n <- 1000              # number of resampling
boot.cor <- replicate(n, tmp <- sample(round(length(x)*percent), replace=FALSE); cor(x[tmp], y[tmp]))

在运行max(boot.cor) 之后。如果所有相关系数都相同,请不要失望:)

【讨论】:

【参考方案2】:

如果您真的想要这样做(去除最大(绝对)残差),那么我们可以使用线性模型来估计最小二乘解和相关残差,然后选择中间的 n%的数据。这是一个例子:

首先,生成一些虚拟数据:

require(MASS) ## for mvrnorm()
set.seed(1)
dat <- mvrnorm(1000, mu = c(4,5), Sigma = matrix(c(1,0.8,1,0.8), ncol = 2))
dat <- data.frame(dat)
names(dat) <- c("X","Y")
plot(dat)

接下来,我们拟合线性模型并提取残差:

res <- resid(mod <- lm(Y ~ X, data = dat))

quantile() 函数可以为我们提供所需的残差分位数。您建议保留 90% 的数据,因此我们需要上下 0.05 分位数:

res.qt <- quantile(res, probs = c(0.05,0.95))

选择那些残差在中间 90% 数据中的观察值:

want <- which(res >= res.qt[1] & res <= res.qt[2])

然后我们可以将其可视化,红点是我们将保留的点:

plot(dat, type = "n")
points(dat[-want,], col = "black", pch = 21, bg = "black", cex = 0.8)
points(dat[want,], col = "red", pch = 21, bg = "red", cex = 0.8)
abline(mod, col = "blue", lwd = 2)

完整数据和所选子集的相关性是:

> cor(dat)
          X         Y
X 1.0000000 0.8935235
Y 0.8935235 1.0000000
> cor(dat[want,])
          X         Y
X 1.0000000 0.9272109
Y 0.9272109 1.0000000
> cor(dat[-want,])
         X        Y
X 1.000000 0.739972
Y 0.739972 1.000000

请注意,在这里我们可能会丢弃非常好的数据,因为我们只选择具有最大正残差的 5% 和具有最大负残差的 5%。另一种方法是选择 absolute 残差最小的 90%:

ares <- abs(res)
absres.qt <- quantile(ares, prob = c(.9))
abswant <- which(ares <= absres.qt)
## plot - virtually the same, but not quite
plot(dat, type = "n")
points(dat[-abswant,], col = "black", pch = 21, bg = "black", cex = 0.8)
points(dat[abswant,], col = "red", pch = 21, bg = "red", cex = 0.8)
abline(mod, col = "blue", lwd = 2)

有了这个略有不同的子集,相关性略低:

> cor(dat[abswant,])
          X         Y
X 1.0000000 0.9272032
Y 0.9272032 1.0000000

另外一点是,即便如此,我们也会丢掉好的数据。您可能希望将 Cook 距离视为异常值强度的度量,并仅丢弃高于特定阈值 Cook 距离的那些值。 Wikipedia 有关于库克距离和建议阈值的信息。 cooks.distance() 函数可用于从mod 中检索值:

> head(cooks.distance(mod))
           1            2            3            4            5            6 
7.738789e-04 6.056810e-04 6.375505e-04 4.338566e-04 1.163721e-05 1.740565e-03

如果您计算 Wikipedia 上建议的阈值并仅删除超过阈值的阈值。对于这些数据:

> any(cooks.distance(mod) > 1)
[1] FALSE
> any(cooks.distance(mod) > (4 * nrow(dat)))
[1] FALSE

库克的距离没有一个超过建议的阈值(考虑到我生成数据的方式,这并不奇怪。)

说了这么多,你为什么要这么做?如果您只是想摆脱数据以改善相关性或生成重要关系,那听起来有点可疑,对我来说有点像数据挖掘。

【讨论】:

非常感谢您提供如此出色的答案!我想这样做的原因如下。我正在对基于复合物的实验结构预测实验观察结果(蛋白质复合物突变时结合能的变化)的各种方法进行基准测试。目标值来自不同质量的各种来源。结构中的错误会严重影响预测。所以我有几个异常值,但是查看各种方法的“修剪”相关性将使我能够更轻松地选择最适合有利情况的方法。【参考方案3】:

这对 OP 来说可能已经很明显了,但只是为了确保...您必须小心,因为尝试最大化相关性实际上可能会包含异常值。 (@Gavin 在他的回答/cmets 中谈到了这一点。)我会首先删除异常值,然后计算相关性。更一般地说,我们希望计算对异常值具有鲁棒性的相关性(R 中有许多这样的方法)。

为了生动地说明这一点,让我们创建两个不相关的向量 xy

set.seed(1)
x <- rnorm(1000)
y <- rnorm(1000)
> cor(x,y)
[1] 0.006401211

现在让我们添加一个异常点(500,500)

x <- c(x, 500)
y <- c(y, 500)

现在任何子集包含异常点的相关性将接近 100%,并且任何足够大的子集排除异常点的相关性将接近于零。特别是,

> cor(x,y)
[1] 0.995741

如果你想估计一个对异常值不敏感的“真实”相关性,你可以试试robust 包:

require(robust)
> covRob(cbind(x,y), corr = TRUE)
Call:
covRob(data = cbind(x, y), corr = TRUE)

Robust Estimate of Correlation: 
            x           y
x  1.00000000 -0.02594260
y -0.02594260  1.00000000

您可以使用covRob 的参数来决定如何修剪数据。 更新:MASS 包中还有 rlm(稳健线性回归)。

【讨论】:

【参考方案4】:

cor 中使用method = "spearman" 对污染具有很强的抵抗力,并且易于实施,因为它只涉及将cor(x, y) 替换为cor(x, y, method = "spearman")

重复 Prasad 的分析,但使用 Spearman 相关性,我们发现 Spearman 相关性确实对此处的污染具有鲁棒性,恢复了潜在的零相关性:

set.seed(1)

# x and y are uncorrelated
x <- rnorm(1000)
y <- rnorm(1000)
cor(x,y)
## [1] 0.006401211

# add contamination -- now cor says they are highly correlated
x <- c(x, 500)
y <- c(y, 500)
cor(x, y)
## [1] 0.995741

# but with method = "spearman" contamination is removed & they are shown to be uncorrelated
cor(x, y, method = "spearman")
## [1] -0.007270813

【讨论】:

spearman 对某些类型的污染具有鲁棒性,即单个高值点完全相关,导致pearson 相关性膨胀。但是,它不会完全抵抗规模较低端的异常值的污染。【参考方案5】:

这是捕获异常值的另一种可能性。使用与 Prasad 类似的方案:

library(mvoutlier)    
set.seed(1)    
x <- rnorm(1000)    
y <- rnorm(1000)    
xy <- cbind(x, y)    
outliers <- aq.plot(xy, alpha=0.975) #The documentation/default says alpha=0.025.  I think the functions wants 0.975   
cor.plot(x, y)    
color.plot(xy)   
dd.plot(xy)   
uni.plot(xy)    

在其他答案中,500 被卡在 x 和 y 的末尾作为异常值。可能,或可能不会导致机器的内存问题,所以我将其丢弃到4以避免。

x1 <- c(x, 4)     
y1 <- c(y, 4)    
xy1 <- cbind(x1, y1)    
outliers1 <- aq.plot(xy1, alpha=0.975) #The documentation/default says alpha=0.025.  I think the functions wants 0.975
cor.plot(x1, y1)    
color.plot(xy1)    
dd.plot(xy1)    
uni.plot(xy1)    

这是来自 x1, y1, xy1 数据的图像:

【讨论】:

我向 mvoutlier 的维护者发送了电子邮件,说明我在上述 aq.plot() 语句中遇到的 alpha 问题。自修复了问题并更新了MVOUTLIER到1.6版(更新了2011年1月14日)cran.r-project.org/web/packages/mvoutlier/index.html span>

以上是关于从相关系数计算中去除异常值的主要内容,如果未能解决你的问题,请参考以下文章

SPSS中pearson(皮尔逊相关系数)看r值还是P值,确定相关性

计算自相关系数acf和偏相关系数pacf

相关系数r的计算公式是啥?

MATLAB如何求相关系数

初学python,怎样用python做pearson相关系数的检验呢,求指导啊

如何理解皮尔逊相关系数