Re: ktab dans ade4 sous R

From: Daniel Chessel (chessel@biomserv.univ-lyon1.fr)
Date: Wed Sep 03 2003 - 16:05:18 MEST


At 17:50 01/09/2003 +0200, you wrote:
>Bonjour,
>
>Essayant pour la première fois un ktab avec le module ade4 de R, je
>n'arrive pas a créer mon fichier ktab.
>
>J'ai un data frame constitué comme suit :
>
>période localite espece nombre
>juin sta1 sp1 201
>juin sta1 sp2 35
>...
>
>Il y a 743 lignes, 5 périodes, 13 localités et 71 espèces qui ne sont
>bien sûr pas présentes partout ni à toutes les périodes.
>
>J'ai voulu utiliser la fonction ktab.data.frame mais manifestement le
>data frame d'entrée doit avoir une autre structure que le mien.
>Quelqu'un peut-il m'aider ou m'orienter vers une doc qui puisse
>m'éclairer ?

J'ai placé un data.frame de ce type à
http://pbil.univ-lyon1.fr/R/donnees/totor.txt
faire

> download.file ("http://pbil.univ-lyon1.fr/R/donnees/totor.txt","toto")
> toto = read.table("toto",h=T)
> summary(toto)
 sai sta esp ab
 sp:22 Min. :1.00 Brh :10 Min. : 2.0
 su:24 1st Qu.:1.00 Bsp : 8 1st Qu.: 4.0
         Median :3.00 Hla : 7 Median : 5.0
         Mean :3.07 Eig : 4 Mean : 5.5
         3rd Qu.:4.75 Rhi : 4 3rd Qu.: 7.0
         Max. :5.00 Bni : 2 Max. :11.0
                        (Other):11

> toto$sta = as.factor(toto$sta)
> summary(toto)
 sai sta esp ab
 sp:22 1:14 Brh :10 Min. : 2.0
 su:24 2: 3 Bsp : 8 1st Qu.: 4.0
         3: 7 Hla : 7 Median : 5.0
         4:10 Eig : 4 Mean : 5.5
         5:12 Rhi : 4 3rd Qu.: 7.0
                Bni : 2 Max. :11.0
                (Other):11

2 saisons (facteur sai), 5 stations (facteur sta), 11 espèces (facteur esp), abondance (numeric ab)
les absences ne sont pas indiquées.
Lorsque les blocs du k-tableau sont définies par une variable qualitative (facteur)
la stratégie générale est celle du passage intra-classe à k tableaux.

1) Il faut d'abord faire un tableau :

> toto$sai:toto$sta
 [1] sp:1 sp:1 sp:1 sp:1 sp:1 sp:1 sp:1 sp:2 sp:2 sp:3 sp:3 sp:3 sp:3 sp:4 sp:4
[16] sp:4 sp:4 sp:5 sp:5 sp:5 sp:5 sp:5 su:1 su:1 su:1 su:1 su:1 su:1 su:1 su:2
[31] su:3 su:3 su:3 su:4 su:4 su:4 su:4 su:4 su:4 su:5 su:5 su:5 su:5 su:5 su:5
[46] su:5
Levels: sp:1 sp:2 sp:3 sp:4 sp:5 su:1 su:2 su:3 su:4 su:5

énumérer les couples espace-temps
> nomrel=as.character(levels(toto$sai:toto$sta))
énumérer les espèces
> nomesp=as.character(levels(toto$esp))
préparer une matrice, en faire un data.frame
> tab=matrix(0,length(nomrel),length(nomesp))
> tab=data.frame(tab)
> row.names(tab)=nomrel
> names(tab)=nomesp
> tab
     Bni Bpu Brh Bsp Cae Cen Ecd Eda Eig Hla Par Rhi
sp:1 0 0 0 0 0 0 0 0 0 0 0 0
sp:2 0 0 0 0 0 0 0 0 0 0 0 0
sp:3 0 0 0 0 0 0 0 0 0 0 0 0
sp:4 0 0 0 0 0 0 0 0 0 0 0 0
sp:5 0 0 0 0 0 0 0 0 0 0 0 0
su:1 0 0 0 0 0 0 0 0 0 0 0 0
su:2 0 0 0 0 0 0 0 0 0 0 0 0
su:3 0 0 0 0 0 0 0 0 0 0 0 0
su:4 0 0 0 0 0 0 0 0 0 0 0 0
su:5 0 0 0 0 0 0 0 0 0 0 0 0

Le remplir :
> nomrow=as.character(toto$sai:toto$sta)
> nomcol=as.character(toto$esp)
> for (i in 1:nrow(toto)) tab[nomrow[i],nomcol[i]]=toto$ab[i]
> tab
     Bni Bpu Brh Bsp Cae Cen Ecd Eda Eig Hla Par Rhi
sp:1 9 0 10 7 0 0 0 4 0 9 4 5
sp:2 0 0 8 0 0 0 0 0 0 4 0 0
sp:3 0 0 5 5 0 0 0 0 0 5 0 2
sp:4 0 0 6 3 0 0 0 0 0 6 0 3
sp:5 0 0 6 5 0 0 5 0 4 4 0 0
su:1 0 10 10 7 0 0 0 6 2 7 0 2
su:2 0 0 9 0 0 0 0 0 0 0 0 0
su:3 0 0 8 6 0 2 0 0 0 0 0 0
su:4 0 0 11 7 5 2 0 0 5 2 0 0
su:5 2 3 9 6 2 0 4 0 7 0 0 0

2) Il faut faire un facteur qui définit les blocs :
> substr(rownames(tab),1,2)
 [1] "sp" "sp" "sp" "sp" "sp" "su" "su" "su" "su" "su"
> sai = as.factor(substr(rownames(tab),1,2))
> sai
 [1] sp sp sp sp sp su su su su su
Levels: sp su

3A) On pense qu'un tableau élémentaire relève d'une ACP centrée (par espèce) :
> pca1=dudi.pca(tab,scale=FALSE,scann=FALSE)
> wit1=within(pca1,sai,scann=FALSE)
> ktab1=ktab.within(wit1)
> ktab1
class: ktab

tab number: 2
  data.frame nrow ncol
1 sp 12 5
2 su 12 5
...
On obtient un k-tableaux de 2 tableaux, espèces en lignes, sites en colonnes, centrés par ligne
Le second est :
$su
    su:1 su:2 su:3 su:4 su:5
Bni -0.4 -0.4 -0.4 -0.4 1.6
Bpu 7.4 -2.6 -2.6 -2.6 0.4
Brh 0.6 -0.4 -1.4 1.6 -0.4
Bsp 1.8 -5.2 0.8 1.8 0.8
Cae -1.4 -1.4 -1.4 3.6 0.6
Cen -0.8 -0.8 1.2 1.2 -0.8
Ecd -0.8 -0.8 -0.8 -0.8 3.2
Eda 4.8 -1.2 -1.2 -1.2 -1.2
Eig -0.8 -2.8 -2.8 2.2 4.2
Hla 5.2 -1.8 -1.8 0.2 -1.8
Par 0.0 0.0 0.0 0.0 0.0
Rhi 1.6 -0.4 -0.4 -0.4 -0.4

3B) On pense qu'un tableau élémentaire relève d'une AFC :

> coa1=dudi.coa(tab,scann=F)
> coa1=dudi.coa(tab,scann=F)
> ktab2=ktab.within(wit2)
> ktab2
class: ktab

tab number: 2
  data.frame nrow ncol
1 sp 12 5
2 su 12 5
...
On obtient un k-tableaux de 2 tableaux, espèces en lignes, sites en colonnes, collection d'AFC sur modèles
voir Dufour, A.-B., S. Pavoine, and D. Chessel. 2003. Introduction of correspondence analysis in multiway methods of simultaneaous ordination. Abstract 25 in M. Greenacre and J. Blasius, editors. International Conference on Correspondence Analysis and Related Methods (CARME 2003), Barcelona
http://pbil.univ-lyon1.fr/R/articles/arti111.pdf
On a ainsi l'AFM pour AFC : Bécue-Bertaut, M., and J. Pagès. 2003. A principal axes method for comparing contingency tables : MFACT. Computational Statistics & Data Analysis In press.

4) Autre façon de procéder à partir des données brutes
Les points 3) sont toujours possibles. Si toutes les stations sont présentes à chaque date,
les tableaux peuvent être juxtaposés par les lignes dans un tableau
> tab1=cbind.data.frame(tab[1:5,],tab[6:10,])
> tab1
     Bni Bpu Brh Bsp Cae Cen Ecd Eda Eig Hla Par Rhi Bni Bpu Brh Bsp Cae Cen
sp:1 9 0 10 7 0 0 0 4 0 9 4 5 0 10 10 7 0 0
sp:2 0 0 8 0 0 0 0 0 0 4 0 0 0 0 9 0 0 0
sp:3 0 0 5 5 0 0 0 0 0 5 0 2 0 0 8 6 0 2
sp:4 0 0 6 3 0 0 0 0 0 6 0 3 0 0 11 7 5 2
sp:5 0 0 6 5 0 0 5 0 4 4 0 0 2 3 9 6 2 0
     Ecd Eda Eig Hla Par Rhi
sp:1 0 6 2 7 0 2
sp:2 0 0 0 0 0 0
sp:3 0 0 0 0 0 0
sp:4 0 0 5 2 0 0
sp:5 4 0 7 0 0 0
On choisit un centrage : ici par espèce et par saison pour retrouver le précédent :
> tab10=data.frame(scalewt(tab1,scale=F))
Enfin, faire un k-tableaux en passant les nombres de colonnes par paquets :
> ktab3=ktab.data.frame(tab10,c(12,12))
> ktab3
class: ktab

tab number: 2
  data.frame nrow ncol
1 Ana1 5 12
2 Ana2 5 12
On obtient un k-tableaux de 2 tableaux, espèces en colonnes, sites en lignes, centrés par colonne
Le second est :

> ktab3$Ana2
      Bni Bpu Brh Bsp Cae Cen Ecd Eda Eig Hla Par Rhi
sp:1 -0.4 7.4 0.6 1.8 -1.4 -0.8 -0.8 4.8 -0.8 5.2 0 1.6
sp:2 -0.4 -2.6 -0.4 -5.2 -1.4 -0.8 -0.8 -1.2 -2.8 -1.8 0 -0.4
sp:3 -0.4 -2.6 -1.4 0.8 -1.4 1.2 -0.8 -1.2 -2.8 -1.8 0 -0.4
sp:4 -0.4 -2.6 1.6 1.8 3.6 1.2 -0.8 -1.2 2.2 0.2 0 -0.4
sp:5 1.6 0.4 -0.4 0.8 0.6 -0.8 3.2 -1.2 4.2 -1.8 0 -0.4

Ici les stratégies sont toutes possibles parce que l'appariement est possible par site et par espèce.
Attention :
> row.names(ktab1)
 [1] "Bni" "Bpu" "Brh" "Bsp" "Cae" "Cen" "Ecd" "Eda" "Eig" "Hla" "Par" "Rhi"
> col.names(ktab1)
 [1] "sp:1" "sp:2" "sp:3" "sp:4" "sp:5" "su:1" "su:2" "su:3" "su:4" "su:5"

> row.names(ktab3)
[1] "sp:1" "sp:2" "sp:3" "sp:4" "sp:5"
> col.names(ktab3)
 [1] "Bni" "Bpu" "Brh" "Bsp" "Cae" "Cen" "Ecd" "Eda" "Eig" "Hla" "Par" "Rhi"
[13] "Bni" "Bpu" "Brh" "Bsp" "Cae" "Cen" "Ecd" "Eda" "Eig" "Hla" "Par" "Rhi"

L'analyse triadique est possible sur le second mais pas sur le premier à cause des noms de colonnes
On a strictement l'identité avec :
col.names(ktab1) = c(c("s1","s2","s3","s4","s5"),c("s1","s2","s3","s4","s5"))
row.names(ktab3) = c("s1","s2","s3","s4","s5")
pta(ktab1,scan=F)$eig
pta(ktab3,scan=F)$eig
> col.names(ktab1) = c(c("s1","s2","s3","s4","s5"),c("s1","s2","s3","s4","s5"))
> row.names(ktab3) = c("s1","s2","s3","s4","s5")
> pta(ktab1,scan=F)$eig
[1] 40.711 17.426 6.626 3.357
> pta(ktab3,scan=F)$eig
[1] 40.711 17.426 6.626 3.357

5) autre façon en assemblant des ACP à pondération commune :
dudi1=dudi.pca(tab[1:5,],scal=F,scan=F)
dudi2=dudi.pca(tab[6:10,],scal=F,scan=F)
ktab4=ktab.list.dudi(list(dudi1,dudi2))
> ktab4
class: ktab

tab number: 2
  data.frame nrow ncol
1 Ana1 5 12
2 Ana2 5 12
ktab4$Ana2
      Bni Bpu Brh Bsp Cae Cen Ecd Eda Eig Hla Par Rhi
sp:1 -0.4 7.4 0.6 1.8 -1.4 -0.8 -0.8 4.8 -0.8 5.2 0 1.6
sp:2 -0.4 -2.6 -0.4 -5.2 -1.4 -0.8 -0.8 -1.2 -2.8 -1.8 0 -0.4
sp:3 -0.4 -2.6 -1.4 0.8 -1.4 1.2 -0.8 -1.2 -2.8 -1.8 0 -0.4
sp:4 -0.4 -2.6 1.6 1.8 3.6 1.2 -0.8 -1.2 2.2 0.2 0 -0.4
sp:5 1.6 0.4 -0.4 0.8 0.6 -0.8 3.2 -1.2 4.2 -1.8 0 -0.4

6) On peut aussi transposer des schémas de dualité avant de les assembler
> dudi1=t.dudi(dudi1)
> dudi2=t.dudi(dudi2)
> ktab5=ktab.list.dudi(list(dudi1,dudi2))
> ktab5
class: ktab

tab number: 2
  data.frame nrow ncol
1 Ana1 12 5
2 Ana2 12 5
> ktab5$Ana2
    su.1 su.2 su.3 su.4 su.5
Bni -0.4 -0.4 -0.4 -0.4 1.6
Bpu 7.4 -2.6 -2.6 -2.6 0.4
Brh 0.6 -0.4 -1.4 1.6 -0.4
Bsp 1.8 -5.2 0.8 1.8 0.8
Cae -1.4 -1.4 -1.4 3.6 0.6
Cen -0.8 -0.8 1.2 1.2 -0.8
Ecd -0.8 -0.8 -0.8 -0.8 3.2
Eda 4.8 -1.2 -1.2 -1.2 -1.2
Eig -0.8 -2.8 -2.8 2.2 4.2
Hla 5.2 -1.8 -1.8 0.2 -1.8
Par 0.0 0.0 0.0 0.0 0.0
Rhi 1.6 -0.4 -0.4 -0.4 -0.4

On obtient des ktab avec les fonctions
"[.ktab" extraction de tableaux comme pour les lignes ou les colonnes dans un tableau
"c.ktab" combinaison de ktab associés par les lignes
"ktab.data.frame" à partir de blocs colonnes par passage des nombres de colonnes par blocs
"ktab.list.df" à partir d'une liste de data.frame
"ktab.list.dudi" à partir d'une liste de schémas de dualité
"ktab.within" à partir d'une intra-classe

Bon courage



This archive was generated by hypermail 2b30 : Tue Sep 07 2004 - 13:45:24 MEST