Re: ACPVI et ACP des valeurs predites

From: Daniel Chessel (chessel@biomserv.univ-lyon1.fr)
Date: Mon Nov 08 1999 - 12:49:46 MET


Pierre WAVRESKY reprend la fiche 3.2 et 3.5 avec quelques difficultés.

Il y a sans doute des imprécisions qui méritent d'être corrigées.

>Utilisateur recent d'ADE, j ai eu quelques difficultes avec son
>utilisation. J ai essaye de suivre la fiche 3.2 "regression lineaire" (pp.1
>a 10), puis la fiche 3.5 "ordination sous contraintes" (a partir de la page
>28, c est a dire quand on fait une ACPVI a partir des tableaux "forme" et
>"Mil").

Dans la fiche 3.2, pour obtenir le fichier Forme, il faut :
1) extraire les variables 2a7 de TLog dans Provi et la colonne 1 dans
Taille (FilesUtil: Row-Col Selection)
2) ajouter à la colonne 1 de Taille une colonne de 1 ce qui donne 1ntaille
(FilesUtil: Add column 1n)
3) définir le sous-espace de projection associé à 1ntaille
(Projectors: Table->Orthonormal Basis)
Ceci donne 1ntaille.@ob
4) projeter provi sur cette base avec le titre Reg
Projectors: Table Projection
Orthonormal basis: E:\Ade4\TRUITE\1ntaille.@ob
It has 33 rows and 2 columns
Dependant variable file: E:\Ade4\TRUITE\Provi
It has 33 rows and 6 columns
|----|----------|----------|----------| |------|------|
| |Subspace A| A Orthogo| Total | | A+ | A- |
|----|----------|----------|----------| |------|------|
| 1|5.9870e-02|1.2362e-02|7.2231e-02| | 8288| 1711|
| 2|8.4546e-03|4.7545e-03|1.3209e-02| | 6400| 3599|
| 3|1.9323e-03|1.8089e-03|3.7412e-03| | 5164| 4835|
| 4|9.8621e+00|4.3938e-01|1.0301e+01| | 9573| 426|
| 5|7.3828e-03|2.0414e-03|9.4243e-03| | 7833| 2166|
| 6|2.7732e+01|4.2685e-01|2.8158e+01| | 9848| 151|
|----|----------|----------|----------| |------|------|
File Reg.mod contains the predicted variables
It has 33 rows and 6 columns

File Reg.coe contains the canonical weights
coefficients of linear combination of explanatory variables
It has 2 rows (explanatory v.) and 6 columns (dependent v.)

File Reg.res contains the (observed - predicted) values
It has 33 rows and 6 columns

5) Renommer Reg.res en Forme
Forme est listé dans la fiche 3.2 p. 10

L'exercice permet de voir que Projectors fait aussi bien de la régression
par l'origine que de la régression ordinaire.

> Comme je ne suis pas arrive a definir le sous espace de projection
>(fiche 3.5, a partir p.28),

Il suffit de faire
1) Normaliser Mil en Mil0 par Bin->Bin: Centring

Le fichier sorti de la carte de données s'appelle Habitat (d'accord, c'est
idiot) et il a été renommé Mil
X input file: E:\Ade4\TRUITE\Mil
--- Number of rows: 33, columns: 12
Y input file: Mil0
--- Number of rows: 33, columns: 12
Uniform row weights
1 = D-centring | 2 = D-standardization | 3 = D-normalization
--- Selected option: 3

2) Définir le sous-espace de projection associé à Mil0 par Projectors:
Table->Orthonormal Basis

Orthonormalization: subspace generated by quantitative variables
------------------------------------------
Explanatory variable file: E:\Ade4\TRUITE\Mil0
It has 33 rows and 12 columns
------------------------------------------
Orthonormal basis: E:\Ade4\TRUITE\Mil0.@ob
It has 33 rows and 12 columns
Row weight file: E:\Ade4\TRUITE\Mil0.@pl
Uniform row weight = 0.030303
Coordinates of the vectors of the orthonormal basis
in the initial basis in : E:\Ade4\TRUITE\Mil0.@co
File E:\Ade4\TRUITE\Mil0.@co has 12 rows and 12 columns
------------------------------------------

3) Faire l'ACP normée de Forme (PCA: Correlation matrix PCA)
Classical Principal Component Analysis (Hotelling 1933)
Input file: E:\Ade4\TRUITE\Forme
File :E:\Ade4\TRUITE\Forme.cnta
|Col.| Mini | Maxi |
|----|----------|----------|
| 1|-1.684e+00| 3.957e+00|
...
| 6|-1.831e+00| 1.874e+00|
...
     means and variances:
     Col.: 1 | Mean: -7.7892e-09 | Variance: 1.2362e-02
     Col.: 2 | Mean: -2.8786e-09 | Variance: 4.7545e-03
     Col.: 3 | Mean: -1.6086e-09 | Variance: 1.8089e-03
     Col.: 4 | Mean: -9.0310e-08 | Variance: 4.3939e-01
     Col.: 5 | Mean: -2.3424e-09 | Variance: 2.0414e-03
     Col.: 6 | Mean: -1.5895e-07 | Variance: 4.2685e-01
...
----------------------- Correlation matrix -------------------
[ 1] 1000
[ 2] 653 1000
[ 3] 314 -249 1000
[ 4] -91 -144 450 1000
[ 5] 137 -231 197 -229 1000
[ 6] -280 -253 -69 -47 607 1000
--------------------------------------------------------------
...
-----------------------
Total inertia: 6
-----------------------
Num. Eigenval. R.Iner. R.Sum |Num. Eigenval. R.Iner. R.Sum |
01 +1.9687E+00 +0.3281 +0.3281 |02 +1.5558E+00 +0.2593 +0.5874 |
03 +1.4786E+00 +0.2464 +0.8339 |04 +6.7654E-01 +0.1128 +0.9466 |
05 +2.2203E-01 +0.0370 +0.9836 |06 +9.8266E-02 +0.0164 +1.0000 |
...

4) Faire l'ACPVI de Forme.cnta sur Mil0.@ob (Projectors: PCA on
Instrumental Variables)

Instrumental variables
------------- input -----------------
Orthonormal basis: E:\Ade4\TRUITE\Mil0.@ob
It has 33 rows and 12 columns
Dependent variable file: E:\Ade4\TRUITE\Forme.cnta
It has 33 rows and 6 columns
------------- output ---------------------
Projected variable file: FM.ivta
It has 33 rows and 6 columns
Inertia: 2.9849e+00
...
Num. Eigenval. R.Iner. R.Sum |Num. Eigenval. R.Iner. R.Sum |
01 +1.1925E+00 +0.3995 +0.3995 |02 +9.0626E-01 +0.3036 +0.7031 |
03 +5.1868E-01 +0.1738 +0.8769 |04 +2.4315E-01 +0.0815 +0.9584 |
05 +1.0420E-01 +0.0349 +0.9933 |06 +2.0089E-02 +0.0067 +1.0000 |

> j ai applique la remarque du bas de la page 33
>de la fiche 3.5 : "on obtient exactement le meme resultat qu avec une ACP
>centree ordinaire du tableau des previsions des variables de Forme par
>celle de Mil".
>J ai donc fait des regressions (sous SAS) des variables de "Forme" en
>fonction des variables de "Mil". J obtiens quasiment la meme chose que la
>fiche 3.5, p.30 : les R2 sont presque (pourquoi presque?) identiques a la
>colonne "Subspace A". Ensuite j ai fait une ACP (avec SAS ou ADE, ca donne
>la meme chose) centree, puis centree reduite des valeurs ajustees (y
>chapeau) de ces regressions.
>Et les valeurs propres de l ACP centree reduite sont les suivantes :
>2.26 1.80 1.16 0.52 0.22 0.04. (avec une somme de 6, pas de 2.98,
>comme dans la fiche).
>le nuage des individus (et le cercle de correlation : biom1+ est legerement
>du cote negatif de l axe 2 au lieu de l etre legerement du cote positif)
>est quasiment identique a celui montre p.34.
>Ou est mon erreur?

Je ne sais pas s'il y a une erreur mais l'exercice vaut la peine. On va le
faire dans S-PLUS (chacun ses manies)
J'ai une petite pile qui s'appelle Splus.sta qui transforme un dossier de
travail d'ADE en dossier de travail S-PLUS
Si il y a des amateurs, on la met sur le serveur

Working data will be in E:\ADE4\Truite\_Data
Importer Forme.txt
> Forme
           V1 V2 V3 V4 V5 V6
 1 0.0339550 -0.043108 -0.005779 -0.38660 -0.010573 -0.42023
 2 0.0251380 -0.006679 0.001024 0.30403 0.040768 1.12990
 3 -0.0205320 0.003173 -0.018778 -0.89336 -0.006634 0.12759

Importer Habitat.txt

> Habitat
     V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12
 1 1017 12.0 169 1.80 5.6 0.23 0.29 0.55 16 139 482 53.0
 2 1010 12.2 162 2.02 8.4 0.20 0.30 0.41 7 104 1281 23.5
 3 970 12.6 155 4.00 7.0 0.30 0.32 0.48 19 120 1350 25.0
 4 830 13.0 165 0.50 5.0 0.30 0.26 0.73 22 727 2575 26.0
 5 800 11.0 93 2.60 14.5 0.32 0.52 1.45 5 321 8205 36.0
...

Pour centrer :
> Mil0_Habitat-matrix(colMeans(Habitat),33,12,byrow=T)
Pour normaliser :
> Mil0_Mil0/sqrt(matrix(colVars(Habitat, unbiased=F),33,12,byrow=T))
> Mil0[5,]
      V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11
 -0.7841 -0.9475 -0.5254 -0.03 1.942 0.6534 2.739 3.199 -1.082 -0.226 3.319
     V12
 -0.1882

Dans Mil0 de ADE-4 on a :
   5 | -0.7841 -0.9475 -0.5254 -0.0300 1.9420 0.6534 2.7390 3.1986
-1.0816 -0.2260 3.3192 -0.1882

> Form0_Forme-matrix(colMeans(Forme),33,6,byrow=T)
> Form0_Form0/sqrt(matrix(colVars(Forme, unbiased=F),33,6,byrow=T))
> Form0[5,]
      V1 V2 V3 V4 V5 V6
 -0.6527 -0.1631 -0.4208 -0.4733 -0.5469 -0.07868
Dans Forme.cnta de ADE-4 on a :
    5 | -0.6527 -0.1631 -0.4208 -0.4733 -0.5469 -0.0787

> Mil0_as.data.frame(Mil0)
> Form0_as.data.frame(Form0)
provi_function(x)
{
        auxi <- cbind.data.frame(x, Mil0)
        predict(lm(X1 ~ V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 + V10 + V11 +
V12, data = auxi))
}

Pour faire les régressions colonne par colonne :
> res_apply(Form0,2,provi)
> res
         V1 V2 V3 V4 V5 V6
 1 0.09311 -0.43525 0.0747983 -0.09552 -0.13966 -0.07874
 2 0.07645 -0.41074 0.4871985 0.05363 0.52001 1.02416
 3 -0.74453 -0.78738 -0.5337949 -0.68343 0.01575 0.53576
 4 0.41889 -0.02212 0.2204857 -0.09913 0.07343 -0.22355
 5 -0.60266 -0.18636 -0.4100500 -0.30374 -0.02488 0.51912

Dans ADE-4
Binary input file: E:\Ade4\TRUITE\FM.ivta - 33 rows, 6 cols.
   1 | 0.0931 -0.4353 0.0748 -0.0955 -0.1397 -0.0788
   2 | 0.0765 -0.4107 0.4872 0.0536 0.5200 1.0241
   3 | -0.7445 -0.7874 -0.5338 -0.6834 0.0158 0.5357
   4 | 0.4189 -0.0221 0.2205 -0.0991 0.0734 -0.2236
   5 | -0.6027 -0.1864 -0.4101 -0.3037 -0.0249 0.5191

Pour les valeurs propres de l'ACP centrée du tableau des prédictions
> (princomp(res,cor=F)$sdev)^2
 Comp. 1 Comp. 2 Comp. 3 Comp. 4 Comp. 5 Comp. 6
   1.193 0.9063 0.5187 0.2431 0.1042 0.02009
A comparer avec :
Num. Eigenval. R.Iner. R.Sum |Num. Eigenval. R.Iner. R.Sum |
01 +1.1925E+00 +0.3995 +0.3995 |02 +9.0626E-01 +0.3036 +0.7031 |
03 +5.1868E-01 +0.1738 +0.8769 |04 +2.4315E-01 +0.0815 +0.9584 |
05 +1.0420E-01 +0.0349 +0.9933 |06 +2.0089E-02 +0.0067 +1.0000 |

C'est nickel.

Pour les valeurs propres de l'ACP normée du tableau des prédictions
> (princomp(res,cor=T)$sdev)^2
 Comp. 1 Comp. 2 Comp. 3 Comp. 4 Comp. 5 Comp. 6
   2.261 1.8 1.161 0.5253 0.2155 0.03813
à comparer avec :
>2.26 1.80 1.16 0.52 0.22 0.04. (avec une somme de 6, pas de 2.98,
>comme dans la fiche).

Moralité :
Une ACPVI est bien une ACP centrée (et non normée) des prédictions séparées
des variables dépendantes par les explicatives. En utilisant
colVars(Habitat, unbiased=F), c'est-à-dire les variances avec 1/n on
reproduit strictement les calculs d'ADE-4 dans S-PLUS.

On retrouve cette nuance dans les R2 :
Dans S-PLUS
> provi <- function(x)
{
        auxi <- cbind.data.frame(x, Mil0)
        lm0 <- lm(X1 ~ V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 + V10 + V11 +
                V12, data = auxi)
        var(lm0$fitted.values)
}
> apply(Form0, 2, provi)
     V1 V2 V3 V4 V5 V6
 0.5343 0.5236 0.5871 0.4152 0.3987 0.6193

> provi <- function(x)
{
        auxi <- cbind.data.frame(x, Mil0)
        lm0 <- lm(X1 ~ V1 + V2 + V3 + V4 + V5 + V6 + V7 + V8 + V9 + V10 + V11 +
                V12, data = auxi)
        var(lm0$fitted.values, unbiased = F)
}
> apply(Form0, 2, provi)
     V1 V2 V3 V4 V5 V6
 0.5181 0.5077 0.5693 0.4026 0.3866 0.6005
Dans ADE-4

Projected inertia on a subspace
------------------------------------------
Orthonormal basis: E:\Ade4\TRUITE\Mil0.@ob
It has 33 rows and 12 columns
Dependent variable file: E:\Ade4\TRUITE\Forme.cnta
It has 33 rows and 6 columns
|---|----------|----------|----------| |-----|-----|
| |Subspace A| A Orthogo| Total | | A+ | A- |
|---|----------|----------|----------| |-----|-----|
| 1|5.1814e-01|4.8186e-01|1.0000e+00| | 5181| 4818|
| 2|5.0770e-01|4.9230e-01|1.0000e+00| | 5077| 4922|
| 3|5.6930e-01|4.3070e-01|1.0000e+00| | 5692| 4307|
| 4|4.0265e-01|5.9735e-01|1.0000e+00| | 4026| 5973|
| 5|3.8660e-01|6.1340e-01|1.0000e+00| | 3865| 6134|
| 6|6.0053e-01|3.9947e-01|1.0000e+00| | 6005| 3994|
|---|----------|----------|----------| |-----|-----|
|Tot|2.9849e+00|3.0151e+00|6.0000e+00| | 4974| 5025|
|---|----------|----------|----------| |-----|-----|

Le "pourquoi presque ?" vient souvent de ces fameuses variances en 1/n ou
1/(n-1)
Bon, tout ça est cohérent, mais c'est quand même plus simple de cliquer sur
le bouton

Cordialement
Daniel Chessel
Universite Lyon 1 - Biométrie et Biologie Evolutive - Bât 741
69622 Villeurbanne CEDEX
Tel : 04 72 44 82 77 - (33) 4 72 44 82 77



This archive was generated by hypermail 2b30 : Sat Feb 10 2001 - 10:36:02 MET