Bonjour,
la précision de l'article non plus.
J'ai pensé que le mieux était de retourner dans mes archives aux
données utilisées. J'ai réactiver d'anciennes piles hypercard sous un
LC630 (quelle jeunesse)... si si cela marche toujours. J'ai pu
retrouvé toutes les données de l'ensemble du fascicule spécial sauf
celles de Friday 87. Je suis alors retourné aux données disponibles
dans ADE.data du metacard windows, j'y ai ajouté les 2 variables
manquantes (en tenant compte des "astuces de décimales") ce qui m'a
mis la puce à l'oreille sur les transformation de variables.
Effectivement, on ne retrouve pas le résultat de la figure 3. Par
contre si on transforme toutes les variables en log naturel alors on
retrouve le résultat à une inversion de signe sur l'axe 2 près. Ceci
tout aussi bien avec la version metacard ADE-4 qu'avec la version
ade4 sous R. Ouf.
Alors par rapport à la reproductibilité des résultats, le problème
est qu'il n'est écrit nulle part dans l'article de 1994 que cette
opération de transformation en log a été effectué. Bref j'ai honte...
(mais j'ai transpiré aussi). J'en profite pendant qu'on y est pour
signaler un autre abus à la fig.2 du même article ou on parle de
"maximal correlation" et "maximal standard deviation". En fait une
"maximal covariance" signifie une optimisation de la correlation et
des variances. C'est la covariance qui est maximale.
Bien à vous,
Sylvain
Le 1 Sep 2006 à 18:39, Tito de Morais Luis a écrit :
> Le mercredi 30 août 2006 à 10:09 +0200, Daniel Chessel a écrit :
>> Oui, ben j'ai encore sauté une info.
>> La précision n'est plus ce qu'elle était.
>
> (...)
>
> Bon, j'ai donc essayé de refaire les calculs de l'article (en partie).
> Mais c'est bien vrai que la précision n'est plus ce qu'elle était, à
> moins que je ne me sois trompé, ce qui est loin d'être à exclure ;-(
>
> Pour les analyses séparés, j'obtiens des résultats (à la rotation
> près)
> voisins mais pas identiques. Pour la co-inertie, cela se gâte un peu.
> J'ai dû me tromper quelque part...
>
> Cela dit, sauf si j'ai vraiment fait une grosse bourde, que je vous
> remercie de relever, ce n'est pas la peine de passer trop de temps la
> dessus. Même si la valeur pédagogique de la chose s'en trouve
> diminuée,
> cela peut convenir ainsi.
>
> Merci déjà pour l'aide apportée.
>
> L. Tito
>
> Voici toujours mon code commenté :
>
> library(ade4)
> data(friday87)
> # Dans les données friday87 fournies avec ade4 il manque la turbidité
> # et la profondeur, qu'il faut donc rajouter
> Dept <- c(7.0, 4.5, 2.5, 6.5, 3.0, 3.0, 2.0, 4.5, 3.0, 3.0, 5.0,
> 4.5, 4.5,11.0, 7.0,2.5)
> Turb <- c(0.8, 3.6, 0.6, 1.0, 3.0, 15.0, 10.0, 15.0, 15.0, 20.0,
> 5.5, 10.0, 60.0, 2.5, 5.7, 3.0)
> friday87$mil <- cbind(friday87$mil,turb,depth)
> # Pour avoir les bons labels du milieu comme dans l'article
> names(friday87$mil) <- c
> ("Parea","Varea","pH","Cond","BOD","THar","Alkal","PO4","NO3","Turb","
> Dept")
> #
> # la pca ("standardized") sur le milieu
> dudi1 <- dudi.pca(friday87$mil, scale = TRUE, center=TRUE, scan =
> FALSE, nf = 2)
> #
> # pour avoir la figure 3 de l'article
> #(sauf 3a qui ne peut être faite faute des données géographiques)
> barplot(dudi1$eig) #3b
> s.corcircle(dudi1$co) #3c
> s.label(dudi1$li) #3d
> #
> # la pca ("centred") sur la faune
> dudi2 <- dudi.pca(friday87$fau, center=TRUE, scan = FALSE, nf = 3)
> #
> # pour avoir la figure 4 de l'article
> barplot(dudi2$eig) #4a
> par(mfrow=c(2,5))
> j=0 ; k=0 ; for(i in as.numeric(friday87$fau.blo)) {j=i+j ; k=k+1 ;
> s.class(dudi2$co[(j-i+1):j,1:2],fac=as.factor(rep(1,i)), cellipse
> =0, label = friday87$tab.names[k],grid = F )} #4b
> # comme le code précédent ne place pas les étoiles à l'origine
> comme dans la fig. 4B du papier
> # j'essaye le suivant mais cela ne fonctionne pas. Je dois me
> tromper quelque part...
> j=0 ; k=0 ; for(i in as.numeric(friday87$fau.blo)) {j=i+j ; k=k+1 ;
> s.class(dudi2$co[(j-i+1):j,1:2],fac=as.factor(rep(1,i)), cellipse
> =0, label = friday87$tab.names[k],grid = F, origin=c(mean(dudi2$co
> [(j-i+1):j,1]) , mean(dudi2$co[(j-i+1):j,2])) )}
> # Par ailleurs l'allure des figures obtenues ci-dessus est un peu
> différente de celle du papier...
> par(mfrow=c(1,1))
> s.value(dudi2$li[,1:2], dudi2$li[,3], csi = 0.75, cleg =
> 0.75,method='greylevel') #4c
> s.label(dudi2$li[,1:2],add.plot=T)
> # à une rotation près c'est proche...
> #
> # coinertie
> # à partir d'ici j'ai des différences assez importantes avec le papier
> #
> coin1 <- coinertia(dudi1,dudi2, scan = FALSE, nf = 2)
> #
> # pour avoir la figure 5
> par(mfrow = c(1,2))
> s.corcircle(coin1$aX) #5a
> s.corcircle(coin1$aY) # 5b
> par(mfrow = c(1,1))
> #
> # Pour obtenir les valeurs du tableau 1 (mais elles sont
> différentes, sauf erreur de ma part)
> summary(coin1)
> #
> # Pour ce qui est des tests...
> randtest.coinertia(coin1, nrepet =100)
> plot(randtest.coinertia(coin1, nrepet =100))
> # Je n'ai pas réussi à faire les tests de la figure 6...
> #
> # figure 7
> barplot(coin1$eig) #7a
> s.arrow(coin1$co) #7b
> par(mfrow = c(1,2))
> s.label(coin1$lX) #7c
> s.label(coin1$lY) #7d
> par(mfrow = c(1,1))
> #
> par(mfrow=c(2,5))
> j=0 ; k=0 ; for(i in as.numeric(friday87$fau.blo)) {j=i+j ; k=k+1 ;
> s.class(coin1$li[(j-i+1):j,1:2],fac=as.factor(rep(1,i)), cellipse
> =0, label = friday87$tab.names[k],grid = F )} #7e
> #
> # Figure 8a
> s.label(cbind(coin1$lX[,1],coin1$lY[,1]), label = row.names(coin1$lX))
> abline(lm(coin1$lY[,1]~coin1$lX[,1]))
> cor(coin1$lY[,1], coin1$lX[,1])
> text(4,-4,labels=paste("R = ", round(cor(coin1$lY[,1], coin1$lX[,
> 1]),digits=2)))
> #
> # Figure 8b
> coi.mil <- coin1$lX
> coi.fau <- coin1$lY
> names(coi.mil) <- c("axe1", "axe2")
> names(coi.fau) <- c("axe1", "axe2")
> s.match(coi.mil, coi.fau,clab=0)
> s.label(coi.mil, clab = 1, add.plot = TRUE)
>
>
>
> --
> Luis Tito de Morais
> IRD - UR070 (RAP)
> Centre IRD de Belair
> Route des hydrocarbures
> BP 1386, Dakar
> Sénégal
> Téléphone : +221 849 36 58
> Fax : +221 832 16 75
> Courriel : tito@ird.sn
> Site web : http://www.ird.sn/activites/rap
> http://www.mpl.ird.fr/suds-en-ligne/ecosys/estuarien/
> estuarien.html
>
This archive was generated by hypermail 2b30 : Wed Sep 06 2006 - 09:04:30 MEST