Exécution:
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")
vs exécution
c <- cov(x, use = "pairwise.complete.obs") cov2cor(c)
donne des résultats différents. Quelqu'un sait pourquoi et lequel donne des résultats corrects? Les deux fonctions appellent du code C ++ que je n'ai pas compris comment analyser.
Données reproductibles:
cor(x, use = "pairwise.complete.obs")`
3 Réponses :
Je ne suis pas tout à fait sûr, mais je suppose que cela dépend du nombre de points de données que vous utilisez pour le calcul. Vos données:
> 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
Ici vous voyez, que a1 a 2 valeurs manquantes. Ainsi, la corrélation entre a1 et a3, et a1 et a2 n'est basée que sur 8 points de données et la corrélation entre a2 et a3 est basée sur les 10 points de données.
De l'aide de ? Cov2cor :
cov2cor scales a covariance matrix into the corresponding correlation matrix efficiently.
Lorsque cov2cor met à l'échelle votre matrice de covariance, il suppose le même nombre de points de données. Lorsque vous n'utilisez que des observations complètes, le résultat est le même:
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
Vous devez donc vous demander de quelles informations vous avez besoin. Si vous utilisez des observations complètes par paires, vous ne devriez pas utiliser cov2cor à mon avis.
Hmm, pourriez-vous clarifier ce que vous entendez par le fait que cov2cor suppose le même nombre de points de données? Peut-être que quelque chose me manque, mais à ma connaissance, cov2cor ne fait pas une telle hypothèse, il met plutôt à l'échelle les covariances simplement en utilisant la diagonale de la matrice de covariance calculée elle-même (selon le cov2cor code).
C'est ce que je voulais dire. Mais la covariance de a1 avec les deux, a2 et a3, n'est calculée qu'avec 8 points de données. La covariance entre a2 et a3 est quant à elle calculée par 10 points de données. Lorsque cov2cor met ensuite à l'échelle la matrice de covariance, il ne connaît pas les valeurs manquantes dans l'ensemble de données d'origine. C'est pourquoi le résultat est le même avec l'option use = complete.obs
J'ai maintenant réussi, avec quelques essais et erreurs et en regardant le code C ++, pour comprendre le problème. Vous aviez raison en ce sens que cela a à voir avec le nombre d'observations, mais le problème ne vient pas des hypothèses de la fonction cov2cor mais c'est plutôt une différence dans la façon dont le cor et les fonctions cov calculent les variances. Je laisserai une réponse complète dans un instant
Il s'avère que la différence réside dans la façon dont les variances sont calculées. Si nous n'avons que les variables x et y où x a des NA, alors la fonction cov calculera le y_mean et var_y en utilisant toutes les observations de y . La fonction cor , quant à elle, calculera y_mean et var_y en utilisant uniquement les valeurs de y qui correspondent à valeurs non manquantes dans x .
Puisque toutes les variables impliquées dans une analyse doivent avoir le même nombre d'observations, cela n'a pas de sens lors du calcul des statistiques de corrélation que l'une des entrées (covariances) se compose de moins d'observations que l'autre entrée ( variances). Peut-être une question d'opinion, et pourrait ne pas avoir tant d'importance tant que la proportion manquante est faible, mais clairement indésirable que ces deux fonctions fournissent silencieusement des résultats différents.
La question StackoverFlow suivante explique comment regarder le code C ++: Comment puis-je voir le code source d'une fonction? et voici un lien vers le code de la fonction C_cov : https: //svn.r-project. org / R / trunk / src / library / stats / src / cov.c
Ci-dessous un R-code qui montre les différences de calcul. Tout d'abord le calcul automatique:
#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
Calculez manuellement
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")
@Jaesoc a donné une excellente explication, voici une intuition légèrement différente:
Intuitivement, lorsque vous avez des valeurs manquantes, vous vous retrouvez avec un nombre différent d'observations pour chaque paire, vous devez donc stocker différentes versions de la variance pour chaque variable (disons var (y) avec tous les obs, var (y) avec le même obs que z, etc.). Maintenant, le problème est que lorsque vous calculez la matrice de covariance, vous ne pouvez stocker qu'une seule version de la variance pour chaque variable (c'est-à-dire les éléments diagonaux). Donc, faire cov2cor () ou même calculer manuellement à partir de la matrice de covariance ne vous donnera pas le résultat correct.
Notez que dans le cas extrême où vous avez moins d'observations que de variables, cov2cor () peut même vous donner des corrélations supérieures à | 1 |!
set.seed(123) dat <- data.frame(x = rnorm(4), y = rnorm(4), z=rnorm(4))#, zz=rnorm(4)) dat[1,3] <- NA dat #> x y z #> 1 -0.56047565 0.1292877 NA #> 2 -0.23017749 1.7150650 -0.4456620 #> 3 1.55870831 0.4609162 1.2240818 #> 4 0.07050839 -1.2650612 0.3598138 ## Count number of pairwise observations: crossprod(as.matrix(!is.na(dat))) #> x y z #> x 4 4 3 #> y 4 4 3 #> z 3 3 3 ## Compare results cov2cor(cov(dat, use = "pairwise.complete.obs")) #> x y z #> x 1.00000000 -0.01630914 0.9632927 #> y -0.01630914 1.00000000 -0.4893276 #> z 0.96329271 -0.48932759 1.0000000 cor(dat, use = "pairwise.complete.obs") #> x y z #> x 1.00000000 -0.01630914 0.9408492 #> y -0.01630914 1.00000000 -0.4005502 #> z 0.94084924 -0.40055017 1.0000000 ## Extreme case: N<p dat_S <- dat[1:3,] cov2cor(cov(dat_S, use = "pairwise.complete.obs")) #> x y z #> x 1.0000000 -0.1777287 1.109409 #> y -0.1777287 1.0000000 -1.060258 #> z 1.1094092 -1.0602577 1.000000 cor(dat_S, use = "pairwise.complete.obs") #> x y z #> x 1.0000000 -0.1777287 1 #> y -0.1777287 1.0000000 -1 #> z 1.0000000 -1.0000000 1
Créé le 26/08/2020 par le package reprex (v0.3.0)