Удалить выбросы из расчета коэффициента корреляции

Question

Предположим, у нас есть два числовых вектора x и y. Коэффициент корреляции Пирсона между x и y определяется как

кор (х, у)

Как я могу автоматически учитывать только подмножество x и y в расчете (скажем, 90%), чтобы максимально увеличить коэффициент корреляции?

Answer 1

Если вы действительно хотите сделать это (удалить самые большие (абсолютные) невязки), то мы можем использовать линейную модель для оценки решения наименьших квадратов и связанных невязок, а затем выбрать среднее значение 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))

The *Функция 1011 * может дать нам требуемые квантили остатков.Вы предложили сохранить 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% с наибольшим отрицательным.Альтернативой является выбор 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

ДругойДело в том, что даже тогда мы выбрасываем хорошие данные.Возможно, вы захотите посмотреть на расстояние Кука как меру силы выбросов и отбросить только те значения, которые превышают определенное пороговое расстояние Кука. Википедия содержит информацию о расстоянии Кука и предлагаемых порогах.Функцию 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

ни одно из расстояний Кука не превышает предложенных порогов (неудивительно, учитывая то, как я генерировал данные.)

Сказав все это, зачем вамсделать это?Если вы просто пытаетесь избавиться от данных, чтобы улучшить корреляцию или создать значительную связь, это звучит немного странно и немного похоже на то, как данные мне нужны.

Answer 2

Использование method = "spearman" в cor будет устойчиво к загрязнению и легко реализуемо, поскольку оно предполагает только замену cor(x, y) на cor(x, y, method = "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

Answer 3

Возможно, это уже было очевидно для ОП, но просто чтобы убедиться ... Вы должны быть осторожны, потому что попытка максимизировать корреляцию может на самом деле иметь тенденцию к включать выбросы. (@Gavin затронул этот вопрос в своем ответе / комментариях.) Я бы сначала удалял выбросы, затем вычислял корреляцию. В более общем смысле, мы хотим вычислить корреляцию, которая является устойчивой к выбросам (и в R существует много таких методов).

Просто чтобы проиллюстрировать это драматически, давайте создадим два некоррелированных вектора x и y:

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, чтобы решить, как обрезать данные. ОБНОВЛЕНИЕ: Существует также rlm (устойчивая линейная регрессия) в пакете MASS.

Answer 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:

Answer 5

Вы можете попробовать загрузить свои данные, чтобы найти самый высокий коэффициент корреляции, например ::100100

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). Не разочаровывайтесь, если все коэффициенты корреляции будут одинаковыми:)