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

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)) 

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% с самым большим отрицательный. Альтернатива заключается в выборе 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 

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

и если вы вычислить порог (s) предложил в Википедии и удалить только те, которые превышают порог. Для этих данных:

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

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

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

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

Answer 3

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

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

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

Теперь давайте добавим точку Outlier (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

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