Gael Millot (millot@ijm.jussieu.fr) me pose directement une question très intéressante :
>J'aurais une petite question sur la fonction dudi.pca dans le package
>ade4. Sauf si je me trompe, j'ai l'impression que quand on demande de
>travailler sur les données initiales centrées réduites, la fonction
>dudi.pca utilise la formule variance sans correction du biais (/n au
>lieu de /(n-1)). On peut le voir dans $tab. Du coup les valeurs propres
>et vecteurs propres différent un peu quand on utilise les données
>centrées uniquement (diagonalisation de la matrice des variances-covariances).
>
>Si je ne me trompe pas, ma question est pourquoi preferer ce type de
>variance dans l'ACP ? Serait-il possible de trouver une fonction
>similaire a dudi.pca mais avec utilisation des variances non biaisées ?
Il m'autorise à dupliquer sur le forum la réponse que je lui fait et qui peut concerner nombre d'utilisateurs.
La question du 1/n ou 1/(n-1) est toujours présente dans une ACP centrée.
Elle recouvre une question sévère qui a pour origine le fait que la définition de l'ACP est multiple.
On peut définir l'ACP comme portant sur une loi de probabilité X multivariée quelconque (c'est un problème de maths) de matrice de covariances C (théorique) qui définit les composantes principales comme les composantes de la variable Y=XU où U contient les vecteurs propres de C (Rao, C. R. 1973. Linear statistical inference and its applications, 2nd edition. Wiley, New York. p. 590).
On peut définir l'ACP comme portant sur un nuage de point quelconque pour lequel on cherche un plan le plus proche des données en diagonalisant la matrice de covariances en 1/n (EXACTEMENT dans Pearson, K. 1901. On lines and planes of closest fit to systems of points in space. Philosophical Magazine 2:559-572)
On peut encore définir l'ACP comme portant sur un échantillon d'une variable gaussienne multivariée. On est alors sur l'estimation au maximum de vraisemblance des composantes principales où on retrouve les estimations sans biais en 1/(n-1) (Anderson, T. W. 1971. An Introduction to Multivariate Statistical Analysis. John Wiley and Sons, New York. p. 460)
On peut encore définir l'ACP comme étant l'estimation des positions de n points dans un sous-espace de dimension k de Rp connus avec une erreur sphérique gaussienne(Besse, P., H. Caussinus, L. Ferre, and J. Fine. 1986. Some guidelines for principal component analysis. Pages 23-30 in I. A. f. S. Computing., editor. COMPSTAT 1986. Physica-Verlag, Heidelberg)
Si on prend 1/(n-1) on pense estimation et loi multivariée gaussienne dont on a un échantillon. Si on prend 1/n on pense axes d'inertie d'un nuage de points.
Pourquoi préférer l'un à l'autre ? dudi.pca est une analyse d'inertie comme dudi.coa, dudi.acm, dudi.mix, ... donc utilise 1/n. C'est tout.
Je propose alors ci-dessous un exercice qui montre que
1 princomp et dudi.pca font exactement la même ACP pondérée en 1/n comme indiquée dans la doc.
2 prcomp fait l'ACP en 1/(n-1)
3 qu'on peut faire une ACP en 1/(n-1) avec as.dudi, ce qui était la seconde partie de la question
as.dudi fait l'analyse d'un tableau X quelconque avec une pondération des colonnes quelconque, une pondération des lignes quelconque : c'est la partie commune à tous les dudi.
Cordialement
library(ade4)
data(euro123)
tab=euro123$in78
tab
pri sec ter
Belgium 0.032 0.359 0.609
Denmark 0.079 0.319 0.602
Spain 0.206 0.372 0.422
France 0.092 0.368 0.540
Greece 0.320 0.297 0.383
Ireland 0.206 0.320 0.474
Italy 0.155 0.381 0.464
Luxembourg 0.062 0.392 0.546
Netherlands 0.054 0.330 0.616
Portugal 0.313 0.348 0.339
Germany 0.061 0.444 0.495
United_Kingdom 0.028 0.390 0.582
pca1=dudi.pca(tab,scale=FALSE,scan=F)
princomp1=princomp(tab)
pca1$eig
[1] 0.017363 0.001983
princomp1$sdev^2
Comp.1 Comp.2 Comp.3
0.017363 0.001983 0.000000
pca1$c1
CS1 CS2
pri 0.7540 0.3133
sec -0.1057 -0.8096
ter -0.6483 0.4963
princomp1$loadings
Loadings:
Comp.1 Comp.2 Comp.3
pri 0.754 0.313 0.577
sec -0.106 -0.810 0.577
ter -0.648 0.496 0.577
pca1$li
Axis1 Axis2
Belgium -0.14358 0.019974
Denmark -0.09938 0.063610
Spain 0.10748 -0.028849
France -0.05456 -0.002761
Greece 0.22664 0.048234
Ireland 0.07926 0.039061
Italy 0.04084 -0.031268
Luxembourg -0.08360 -0.028613
Netherlands -0.12847 0.053820
Portugal 0.24450 -0.017089
Germany -0.05679 -0.096340
United_Kingdom -0.13237 -0.019779
princomp1$scores
Comp.1 Comp.2 Comp.3
Belgium -0.14358 0.019974 -4.398e-18
Denmark -0.09938 0.063610 -8.321e-18
Spain 0.10748 -0.028849 4.049e-17
France -0.05456 -0.002761 4.433e-17
Greece 0.22664 0.048234 4.077e-17
Ireland 0.07926 0.039061 1.032e-17
Italy 0.04084 -0.031268 6.403e-17
Luxembourg -0.08360 -0.028613 6.603e-17
Netherlands -0.12847 0.053820 -4.743e-19
Portugal 0.24450 -0.017089 6.402e-17
Germany -0.05679 -0.096340 6.372e-17
United_Kingdom -0.13237 -0.019779 7.237e-18
princomp et dudi.pca font exactement la même ACP pondérée en 1/n comme indiquée dans la doc.
tab0=scalewt(tab,scale=F)
tab0=as.data.frame(tab0)
tab0
pri sec ter
Belgium -0.102 -0.001 0.103
Denmark -0.055 -0.041 0.096
Spain 0.072 0.012 -0.084
France -0.042 0.008 0.034
Greece 0.186 -0.063 -0.123
Ireland 0.072 -0.040 -0.032
Italy 0.021 0.021 -0.042
Luxembourg -0.072 0.032 0.040
Netherlands -0.080 -0.030 0.110
Portugal 0.179 -0.012 -0.167
Germany -0.073 0.084 -0.011
United_Kingdom -0.106 0.030 0.076
cov(tab)
pri sec ter
pri 0.010981 -0.0020578 -0.0089229
sec -0.002058 0.0016295 0.0004284
ter -0.008923 0.0004284 0.0084945
t(tab0)%*%tab0/(nrow(tab0)-1)
pri sec ter
pri 0.010981 -0.0020578 -0.0089229
sec -0.002058 0.0016295 0.0004284
ter -0.008923 0.0004284 0.0084945
tab0.rw=rep(1/(nrow(tab0)-1),nrow(tab0))
tab0.rw
[1] 0.09091 0.09091 0.09091 0.09091 0.09091 0.09091 0.09091 0.09091 0.09091 0.09091 0.09091 0.09091
tab0.cw=rep(1,ncol(tab0))
tab0.cw
[1] 1 1 1
pca2=as.dudi(tab0,tab0.cw,tab0.rw,F,2,match.call(),"standard")
prcomp1=prcomp(tab,scal=F)
pca2$eig
[1] 0.018942 0.002163
prcomp1$sdev^2
[1] 1.894e-02 2.163e-03 1.288e-33
pca2$c1
CS1 CS2
pri 0.7540 0.3133
sec -0.1057 -0.8096
ter -0.6483 0.4963
prcomp1$rotation
PC1 PC2 PC3
pri 0.7540 0.3133 -0.5774
sec -0.1057 -0.8096 -0.5774
ter -0.6483 0.4963 -0.5774
pca2$li
Axis1 Axis2
Belgium -0.14358 0.019974
Denmark -0.09938 0.063610
Spain 0.10748 -0.028849
France -0.05456 -0.002761
Greece 0.22664 0.048234
Ireland 0.07926 0.039061
Italy 0.04084 -0.031268
Luxembourg -0.08360 -0.028613
Netherlands -0.12847 0.053820
Portugal 0.24450 -0.017089
Germany -0.05679 -0.096340
United_Kingdom -0.13237 -0.019779
prcomp1$x
PC1 PC2 PC3
Belgium -0.14358 0.019974 -2.053e-17
Denmark -0.09938 0.063610 -2.737e-17
Spain 0.10748 -0.028849 -1.012e-17
France -0.05456 -0.002761 -4.555e-17
Greece 0.22664 0.048234 -4.166e-17
Ireland 0.07926 0.039061 -2.386e-17
Italy 0.04084 -0.031268 -4.190e-17
Luxembourg -0.08360 -0.028613 -5.257e-17
Netherlands -0.12847 0.053820 -3.450e-17
Portugal 0.24450 -0.017089 -3.150e-17
Germany -0.05679 -0.096340 -1.026e-17
United_Kingdom -0.13237 -0.019779 -1.323e-17
Daniel Chessel - chessel@biomserv.univ-lyon1.fr
This archive was generated by hypermail 2b30 : Tue Sep 07 2004 - 13:45:25 MEST