cor и cov2cor разные результаты с использованием = pairwise.complete.obs

Question

Бег:

cor(x, use = "pairwise.complete.obs")`

против бега

c <- cov(x, use = "pairwise.complete.obs")
cov2cor(c)

дают разные результаты. Кто-нибудь знает, почему и какой дает правильные результаты? Обе функции вызывают код C ++, который я не понял, как анализировать.

Воспроизводимые данные:

x <- data.frame(a1 = rnorm(10), a2 = rnorm(10), a3 = rnorm(10))
x$a1[c(1,3)] <- NA

c <- cov(x, use = "pairwise.complete.obs")
cov2cor(c)
cor(x, use = "pairwise.complete.obs")

Answer 1

Я не совсем уверен, но мое первое предположение, что это зависит от количества точек данных, которые вы используете для расчета. Ваши данные:

x <- data.frame(a1 = rnorm(10), a2 = rnorm(10), a3 = rnorm(10))
x
           a1          a2           a3
1          NA -0.13838924 -0.321692757
2  -0.4508542  0.94765320  1.951455501
3          NA  2.31819262 -1.411309267
4   0.7828306 -0.53437879 -1.073991019
5   0.2298984 -0.46396636 -0.008471517
6   0.6543559 -1.67425582  0.163433255
7   0.1454043  0.37971651  1.197296537
8   0.2055262  0.67526617 -0.913544378
9   0.3360491  1.84172088  0.366159132
10  0.2507107  0.08284055  1.819004908

Здесь вы видите, что a1 имеет 2 пропущенных значения. Таким образом, корреляция между a1 и a3 и a1 и a2 основана только на 8 точках данных, а корреляция между a2 и a3 основана на всех 10 точках данных.

С помощью ?cov2cor:

cov2cor scales a covariance matrix into the corresponding correlation matrix efficiently.

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

> c <- cov(x, use = "complete.obs")
> cov2cor(c)
           a1         a2        a3
a1  1.0000000 -0.2230619 0.2580796
a2 -0.2230619  1.0000000 0.6152059
a3  0.2580796  0.6152059 1.0000000
> cor(x, use = "complete.obs")
           a1         a2        a3
a1  1.0000000 -0.2230619 0.2580796
a2 -0.2230619  1.0000000 0.6152059
a3  0.2580796  0.6152059 1.0000000

Поэтому вам нужно спросить себя, какая информация вам нужна. Если вы используете попарно полные наблюдения, вы не должны использовать cov2cor по моему мнению.

Answer 2

Сводка

Оказывается, разница в том, как вычисляются отклонения. Если у нас есть только переменные x и y, где x имеет несколько NA, то функция cov будет вычислять y_mean и var_y, используя все наблюдения в y. Функция cor, с другой стороны, будет вычислять y_mean и var_y, используя только значения в y, которые соответствуют не пропущенным значениям в x.

Для меня "pairwise.complete.obs" должен означать, что должна быть вся статистика, которая может быть вычислена с учетом всех наблюдений, и, таким образом, функция cor некорректна для парных удалений. Возможно, вопрос мнения, и, возможно, не имеет большого значения при большем количестве наблюдений, но явно нежелательно, чтобы эти две функции молча давали разные результаты.

Мотивация

Следующий вопрос StackoverFlow объясняет, как смотреть на код C ++: Как просмотреть исходный код функции? и вот ссылка на коддля функции C_cov: https://svn.r -project.org / R / trunk / src / library / stats / src / cov.c

Ниже приведен R-код, демонстрирующийразличия в вычислениях. Сначала автоматический расчет:

dat <- data.frame(x = rnorm(10), y = rnorm(10))
dat$x[c(1,3)] <- NA

d <- cov(dat, use = "pairwise.complete.obs")
cov_version <- cov2cor(d)
cor_version <- cor(dat, use = "pairwise.complete.obs")

Вычислить вручную

#Computing the covariance is the same for both cor and cov functions
n <- length(dat$x[!is.na(dat$x)]) #n of non-missing values
x_bar <- mean(dat$x, na.rm = TRUE) 
y_bar <- mean(dat$y[!is.na(dat$x)]) #mean of y values where x is not NA
n1 <- n -1

s_x <- dat$x[!is.na(dat$x)] - x_bar
s_y <- dat$y[!is.na(dat$x)] - y_bar
ss_xy <- sum(s_x*s_y) #sums of squares (x, y)
cov_xy <- ss_xy / n1 #same results as cov(dat, use = "pairwise.complete.obs")
all.equal(cov_xy, d[1,2])
[1] TRUE

#The cor approach to computing variances
ss_x <- sum(s_x*s_x)
ss_y <- sum(s_y*s_y)
var_x <- ss_x / n1
var_y <- ss_y / n1

cor_xy <- cov_xy / (sqrt(var_x) *sqrt(var_y)) #same as cor(dat, use = "pairwise.complete.obs")
all.equal(cor_xy, cor_version[1,2])
[1] TRUE

#The cov approach to computing variances, different for the variable without missing values
y_bar2 <- mean(dat$y) #all values of y
s_y2 <- dat$y - y_bar2 
ss_y2 <- sum(s_y2*s_y2)
var_y2 <- ss_y2 / (length(dat$y) - 1) #n comes from all values of y
cov_cor <- cov_xy / (sqrt(var_x) * sqrt(var_y2))

all.equal(cov_cor, cov_version[1,2])
[1] TRUE