, Не разочаровывайтесь, если все коэффициенты корреляции будут одинаковыми :)

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

кор (х, у)

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

 Leo12 янв. 2011 г., 10:18
@ Гавин Здесь я считаю, что самые большие остатки являются выбросами.
 Gavin Simpson12 янв. 2011 г., 10:13
Что вы считаете выбросом здесь? Отклонение от линии соответствия наименьших квадратов (т.е. наибольшие остатки) или значения в крайних точках двумерного распределенияx а такжеy?

Ответы на вопрос(5)

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
 Prasad Chalasani12 янв. 2011 г., 14:34
+1 за указание наspearman
 cashoes07 авг. 2014 г., 23:57
spearman будет устойчивым к некоторым типам загрязнения, а именно, к точным точным соотношениям отдельных точек высокого значения, приводящим кpearson корреляция. Это не будет полностью устойчиво к загрязнению выбросами в нижней части шкалы.

это уже было очевидно для ОП, но только для того, чтобы убедиться ... Вы должны быть осторожны, потому что попытка максимизировать корреляцию может на самом деле иметь тенденциювключают останцы. (@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 пакет.

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

Решение Вопроса

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

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

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

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

 Leo12 янв. 2011 г., 11:15
Большое спасибо за такой отличный ответ! Причина, по которой я хочу это сделать, заключается в следующем. Я тестирую различные методы прогнозирования экспериментальных наблюдений (изменения энергии связи при мутации белкового комплекса) на основе экспериментальных структур комплексов. Целевые значения поступают из различных источников с различным качеством. И ошибки в структурах могут серьезно повлиять на прогнозы. Таким образом, у меня есть несколько выбросов, но просмотр «сокращенной» корреляции для различных методов позволит мне легче выбрать метод, который лучше всего подходит для благоприятных случаев.

аналогичную 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:

 bill_08016 янв. 2011 г., 02:47
Я написал по электронной почте сопровождающему mvoutlier о проблеме, с которой я столкнулся с альфой в приведенных выше выражениях aq.plot (). С тех пор он исправил проблему и обновил mvoutlier до версии 1.6 (обновлено 14 января 2011 г.)cran.r-project.org/web/packages/mvoutlier/index.html

Ваш ответ на вопрос