RV or not RV ?

From: Stephane DRAY (dray@biomserv.univ-lyon1.fr)
Date: Wed Sep 29 2004 - 23:43:04 MEST

  • Next message: Pierre BADY: "Re: RV or not RV ?.... that is the question ?"

    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/
    --------------------------------------------------------------------------------------------------



    This archive was generated by hypermail 2b30 : Wed Sep 29 2004 - 23:44:43 MEST