Merci beaucoup de la rapidite et de l'efficacite de votre reponse !!
Bien a vous.
Gael.
-----Message d'origine-----
De : owner-adelist@biomserv.univ-lyon1.fr [mailto:owner-adelist@biomserv.univ-lyon1.fr] De la part
de Daniel Chessel
Envoyé : jeudi 6 novembre 2003 18:07
À : adelist@biomserv.univ-lyon1.fr
Objet : RE : dudi.pca et variance biaisee
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