Re: RV or not RV ?.... that is the question ?

From: Pierre BADY (pierre.bady@univ-lyon1.fr)
Date: Thu Sep 30 2004 - 01:29:02 MEST

  • Next message: Daniel Chessel: "Re: RV or not RV ?.... that is the question ?"

    Bonjour,

    Il me semble que la différence entre les deux matrices de RV se justifie par les
    « types de tableaux » utilisés :

    - La pta (ou statis sur X) établie sa matrice de RV sur la base des tableaux.
    - Statis (statis sur les WD, par défaut dans R) établie sa matrice de RV à
    partir des opérateurs d’inertie (Escoufier 1973, Lavit 1988).

    Escoufier, Y. (1973) Le traitement des variables vectorielles. Biometrics, 29,
    750-760.
    Lavit, C. (1988) Analyse conjointe de tableaux quantitatifs Masson, Paris.

    Si on reprend la fonction myRV2 et qu’on considère les deux premiers tableaux de
    kta2 (cf data meaudret), on obtient :
    #-----------------------------------------------
    myRV2=function(x,y) {
    tr=function(x) sum(diag(x))
    phi1=x
    phi2=y
    res=tr(t(phi1)%*%phi2)/sqrt(tr(t(phi1)%*%phi1)*tr(t(phi2)%*%phi2))
    return(res)
    }

    data(meaudret)
    wit1 <- within.pca(meaudret$mil, meaudret$plan$dat, scan = FALSE, scal =
    "partial")
    kta1 <- ktab.within(wit1, colnames = rep(c("S1","S2","S3","S4","S5"), 4))
    kta2 <- t(kta1)

    t1 <- as.matrix(kta2[[1]])
    t2 <- as.matrix(kta2[[2]])

    # calcul du RV sur les tableaux (aux poids près).
    pta1 <- pta(kta2, scann = FALSE)
    pta1$RV

    # spring summer autumn winter
    # spring 1.0000000 0.6934558 0.7886185 0.2834592
    # summer 0.6934558 1.0000000 0.7671756 0.5340456
    # autumn 0.7886185 0.7671756 1.0000000 0.4794976
    # winter 0.2834592 0.5340456 0.4794976 1.0000000

    myRV2(t1,t2)
    # [1] 0.6934558

    # calcul du RV sur les opérateurs (aux poids près).
    # WD = X Q t(X) D

    statis1 <- statis(kta2, scann = FALSE)
    statis1$RV

    # spring summer autumn winter
    # spring 1.0000000 0.6770721 0.8735325 0.3733631
    # summer 0.6770721 1.0000000 0.8856941 0.8571096
    # autumn 0.8735325 0.8856941 1.0000000 0.6355403
    # winter 0.3733631 0.8571096 0.6355403 1.0000000

    w1 <- t2%*%t(t2)
    w2 <- t1%*%t(t1)
    myRV2(w1,w2)
    # [1] 0.6934558

    #-----------------------------------------------

    Donc au final, on peut considérer que les RV sont calculés de la même façon dans
    les deux analyses. La différence est due à la considération des tableaux dans un
    cas et des opérateurs d&#8217;inertie dans l&#8217;autre.

    En espérant que cela aidera.

    P.BADY

    PS : bonne nuit ... ;o))

    Quoting Stephane DRAY <dray@biomserv.univ-lyon1.fr>:

    > Bonjour,
    > j'ai trouve quelque chose de bizarre, il me semble que c'est un bug...
    > ou
    > alors, j'aurais besoin d'un explication.
    > Une petite manip dans R:
    >
    > > installed.packages()[2,]
    > Package LibPath
    > Version
    > "ade4" "C:/Rdev/R-1.9.1/library"
    > "1.2-2"
    > Priority Bundle
    > Depends
    > NA NA
    > NA
    >
    >
    > data(meaudret)
    > wit1 <- within.pca(meaudret$mil, meaudret$plan$dat, scan =
    > FALSE,
    > scal = "partial")
    > kta1 <- ktab.within(wit1, colnames =
    > rep(c("S1","S2","S3","S4","S5"), 4))
    > kta2 <- t(kta1)
    > pta1 <- pta(kta2, scann = FALSE)
    > stati1 <- statis(kta2, scann = FALSE)
    >
    > > pta1$RV
    > spring summer autumn winter
    > spring 1.0000000 0.6934558 0.7886185 0.2834592
    > summer 0.6934558 1.0000000 0.7671756 0.5340456
    > autumn 0.7886185 0.7671756 1.0000000 0.4794976
    > winter 0.2834592 0.5340456 0.4794976 1.0000000
    > > stati1$RV
    > spring summer autumn winter
    > spring 1.0000000 0.6770721 0.8735325 0.3733631
    > summer 0.6770721 1.0000000 0.8856941 0.8571096
    > autumn 0.8735325 0.8856941 1.0000000 0.6355403
    > winter 0.3733631 0.8571096 0.6355403 1.0000000
    >
    > Arghh... Les RV ne sont pas les memes.
    >
    > Je me fais une petite fonction, pour calculer le RV (d'apres l'article
    > de
    > Escoufier, je reprends la notation de Heo et Gabriel):
    >
    > myRV=function(x,y) {
    > tr=function(x) sum(diag(x))
    > phi1=x%*%t(x)
    > phi2=y%*%t(y)
    > res=tr(t(phi1)%*%phi2)/sqrt(tr(phi1%*%phi1)*tr(phi2%*%phi2))
    > return(res)}
    >
    > Et je vois ce que cela donne:
    > > myRV(as.matrix(kta2[[1]]),as.matrix(kta2[[2]]))
    > [1] 0.6770721
    >
    > La fonction statis semble donc donner les bons resultats. Il y aurait
    > donc
    > un probleme dans pta ??
    > J'ai ete voir dans le code, et je refait une nouvelle fonction:
    >
    > myRV2=function(x,y) {
    > tr=function(x) sum(diag(x))
    > phi1=x
    > phi2=y
    > res=tr(t(phi1)%*%phi2)/sqrt(tr(t(phi1)%*%phi1)*tr(t(phi2)%*%phi2))
    > return(res)}
    >
    > > myRV2(as.matrix(kta2[[1]]),as.matrix(kta2[[2]]))
    > [1] 0.6934558
    >
    > J'obtiens les memes resultats que pta !
    >
    > D'une question de logiciel, on passe a une question theorique.
    > Dans la these de Laurence Blanc (avec reference a Escoufier) , on
    > definit:
    >
    > COVV(X1,X2)=trace(X1tX2)
    > VAV(X1)=trace(X1tX1)
    >
    > RV1=COVV(X1,X2)/sqrt(VAV(X1)*VAV(X2))
    >
    > C'est ce qu'utilise pta. On retrouve la meme definition dans
    > http://pbil.univ-lyon1.fr/R/cours/bsc.pdf sans la reference a Escoufier
    >
    > (detail important je pense !).
    >
    > Par contre la definition original du RV est :
    >
    > RV2=COVV(X1,X2)/sqrt(VAV(X1)*VAV(X2))
    >
    > avec :
    >
    > COVV(X1,X2)=trace(X1tX2X2tX1)
    > VAV(X1)=trace(X1tX1X1tX1)
    >
    > Et c'est ce qu'utilise statis.
    >
    > Je n'ai pas les references sur les methodes sous la main. Mais s'agit-il
    >
    > d'un bug (j'en doute a la fin de ce mail) ou plutot d'un choix delibere
    > de
    > nommer deux mesures differentes par la meme appelation ?
    >
    > le RV2 est compris entre 0 et 1 alors que RV1 peut etre negatif.
    > Si on considere deux variables, RV1 est egal a cor(x,y) alors que RV2
    > vaut
    > son carre.
    > Le terme de correlation vectorielle, serait donc peut etre plus
    > approprie
    > au RV1 ?
    >
    >
    > Merci d'avance !
    >
    >
    >
    > Stéphane DRAY
    >
    --------------------------------------------------------------------------------
    ------------------
    >
    >
    > Département des Sciences Biologiques
    > Université de Montréal, C.P. 6128, succursale centre-ville
    > Montréal, Québec H3C 3J7, Canada
    >
    > Tel : (514) 343-6111 poste 1233 Fax : (514) 343-2293
    > E-mail : stephane.dray@umontreal.ca
    >
    --------------------------------------------------------------------------------
    ------------------
    >
    >
    > Web
    > http://www.steph280.freesurf.fr/
    >
    --------------------------------------------------------------------------------
    ------------------
    >
    >
    >
    >
    >
    >

    ---------------------------------------------------------
    Pierre BADY <°))))><
    Université Claude Bernard Lyon 1
    UMR CNRS 5023, LEHF
    bat Alphonse Forel
    43 boulevard du 11 novembre 1918
    F-69622 VILLEURBANNE CEDEX
    FRANCE
    TEL : +33 (0)4 72 44 62 34
    FAX : +33 (0)4 72 43 28 92
    MEL : pierre.bady@univ-lyon1.fr
    http://limnologie.univ-lyon1.fr



    This archive was generated by hypermail 2b30 : Thu Sep 30 2004 - 01:29:34 MEST