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 : Fri Sep 01 2006 - 18:40:55 MEST