Pourquoi les fonctions R 'princomp' et 'prcomp' donnent-elles des valeurs propres différentes?

22

Vous pouvez utiliser le jeu de données de décathlon {FactoMineR} pour le reproduire. La question est de savoir pourquoi les valeurs propres calculées diffèrent de celles de la matrice de covariance.

Voici les valeurs propres en utilisant princomp:

> library(FactoMineR);data(decathlon)
> pr <- princomp(decathlon[1:10], cor=F)
> pr$sd^2
      Comp.1       Comp.2       Comp.3       Comp.4       Comp.5       Comp.6 
1.348073e+02 2.293556e+01 9.747263e+00 1.117215e+00 3.477705e-01 1.326819e-01 
      Comp.7       Comp.8       Comp.9      Comp.10 
6.208630e-02 4.938498e-02 2.504308e-02 4.908785e-03 

Et la même chose en utilisant PCA:

> res<-PCA(decathlon[1:10], scale.unit=FALSE, ncp=5, graph = FALSE)
> res$eig
          eigenvalue percentage of variance cumulative percentage of variance
comp 1  1.348073e+02           79.659589641                          79.65959
comp 2  2.293556e+01           13.552956464                          93.21255
comp 3  9.747263e+00            5.759799777                          98.97235
comp 4  1.117215e+00            0.660178830                          99.63252
comp 5  3.477705e-01            0.205502637                          99.83803
comp 6  1.326819e-01            0.078403653                          99.91643
comp 7  6.208630e-02            0.036687700                          99.95312
comp 8  4.938498e-02            0.029182305                          99.98230
comp 9  2.504308e-02            0.014798320                          99.99710
comp 10 4.908785e-03            0.002900673                         100.00000

Pouvez-vous m'expliquer pourquoi les valeurs propres directement calculées diffèrent de celles-ci? (les vecteurs propres sont les mêmes):

> eigen(cov(decathlon[1:10]))$values
 [1] 1.381775e+02 2.350895e+01 9.990945e+00 1.145146e+00 3.564647e-01
 [6] 1.359989e-01 6.363846e-02 5.061961e-02 2.566916e-02 5.031505e-03

De plus, la prcompméthode alternative donne les mêmes valeurs propres que le calcul direct:

> prc <- prcomp(decathlon[1:10])
> prc$sd^2
 [1] 1.381775e+02 2.350895e+01 9.990945e+00 1.145146e+00 3.564647e-01
 [6] 1.359989e-01 6.363846e-02 5.061961e-02 2.566916e-02 5.031505e-03

Pourquoi faire PCA/ princompet prcompdonner différentes valeurs propres?

George Dontas
la source
L'ACP vous donnera des résultats différents selon que vous utilisez la matrice de covariance ou la matrice de corrélation.
charles.y.zheng
7
Les différences semblent relativement faibles, mais probablement trop importantes pour être de simples problèmes numériques. Serait-ce la différence entre la normalisation par ou n - 1 , par exemple, lors du calcul d'une estimation de la covariance avant le calcul de la décomposition SVD ou valeurs propres? nn-1
Cardinal
7
@cardinal Nice guess! Notez que les deux séquences différentes de valeurs propres ont des rapports successifs identiques. Ainsi, un ensemble est un multiple constant de l'autre. Le multiple est 1,025 = 41/40 ( exactement ). Je ne sais pas d'où cela vient. Peut-être que l'ensemble de données comprend 41 éléments et que l'OP ne révèle que les 10 premiers?
whuber
7
@cardinal Indeed: Page d'aide pour princomp: "Notez que le calcul par défaut utilise le diviseur N pour la matrice de covariance." Page d'aide pour prcomp: "Contrairement à princomp, les variances sont calculées avec le diviseur N-1 habituel."
caracal
2
@caracal, vous devez copier votre commentaire dans une réponse (et peut-être le faire CW) afin qu'il puisse être accepté et la question peut être marquée comme résolue.
cardinal du

Réponses:

16

princompNprcompcovN-1N

Ceci est mentionné dans la section Détails de help(princomp):

Notez que le calcul par défaut utilise le diviseur 'N' pour la matrice de covariance.

et la section Détails de help(prcomp):

À la différence princomp, les variances sont calculées avec le diviseur habituel N - 1.

princompNn.obscv

else if (is.null(covmat)) {
    dn <- dim(z)
    if (dn[1L] < dn[2L]) 
        stop("'princomp' can only be used with more units than variables")
    covmat <- cov.wt(z)
    n.obs <- covmat$n.obs
    cv <- covmat$cov * (1 - 1/n.obs)
    cen <- covmat$center
}

Vous pouvez éviter cette multiplication en spécifiant l' covmatargument au lieu de l' xargument.

princomp(covmat = cov(iris[,1:4]))$sd^2

Mise à jour concernant les scores PCA:

cor = TRUEprincompprincompzN

princomp(scale(data))$scoresprincomp(data, cor = TRUE)$scores(N-1)/N

Joshua Ulrich
la source
1
Vous pouvez envisager de remplacer «deviné» par «déjà confirmé» (voir le flux de commentaires ci-dessus.) Vous pouvez également envisager de modifier votre réponse pour en faire CW. À votre santé.
cardinal du
@cardinal Je n'ai pas vu ces commentaires. Je n'ai vu que ceux qui avaient été votés. Merci. Pouvez-vous également expliquer la justification de la réponse CW? Quelles sont les règles / directives pour cela?
Joshua Ulrich
Quelqu'un peut-il deviner pourquoi le code ne cv <- cov.wt(z, method="ML")rend pas simplement inutiles les 2 lignes suivantes?
caracal
2
@Joshua: Ma suggestion concernant la réponse CW était due au fait que la réponse est apparue via un flux de commentaires et a été générée par une discussion "communautaire". Puisqu'il a été résolu dans les commentaires, je pense qu'il est plus logique de le reformuler comme une réponse, marquée comme CW pour indiquer cette collaboration, et cela permet à la réponse d'être acceptée et à la question d'être marquée comme résolue. (Sinon, il sera automatiquement sauvegardé par le logiciel après un certain temps.)
Cardinal
2
@amoeba, il aurait été utile de le mentionner dans votre commentaire d'édition. "Ajout de 860 caractères au corps" à une réponse de ~ 450 caractères n'aide personne à évaluer si la modification est raisonnable.
Joshua Ulrich