首页 文章

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

提问于
浏览
8

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

cor(x,y)

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

5 回答

  • 9

    如果你真的想这样做(删除最大(绝对)残差),那么我们可以使用线性模型来估计最小二乘解和相关残差,然后选择中间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)
    

    The plot produced from the dummy data showing the selected points with the smallest residuals

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

    > 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% . 另一种方法是选择具有最小绝对残差的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
    

    另一点是,即使这样,我们也会抛出好的数据 . 你可能想看看库克_900151的距离 . 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
    

    如果您计算维基百科上建议的阈值,只删除那些超过阈值的阈值 . 对于这些数据:

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

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

    说完这一切之后,你为什么要这样做呢?如果你只是试图摆脱数据以改善相关性或产生重要关系,那听起来有点可疑,有点像数据疏浚给我 .

  • 23

    cor 中使用 method = "spearman" 将对污染很稳健,并且易于实施,因为它只涉及用 cor(x, y, method = "spearman") 替换 cor(x, y) .

    重复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
    
  • 15

    对于OP来说这可能已经很明显了,但只是为了确保......你必须要小心,因为尝试最大化相关可能实际上往往包括异常值 . (@Gavin在他的回答/评论中触及了这一点 . )我将首先删除异常值,然后计算相关性 . 更一般地,我们想要计算对异常值稳健的相关性(并且在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
    

    如果要估计对异常值不敏感的"true"相关性,可以尝试 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 的参数来决定如何修剪数据 . UPDATE: MASS 包中还有 rlm (稳健线性回归) .

  • 4

    这是捕获异常值的另一种可能性 . 使用与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数据中的图像:

    alt text

    alt text

    alt text

  • 2

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

    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) 之后 . 如果所有相关系数都相同,请不要失望:)

相关问题