Re: analyses des correspondances internes

From: Daniel Chessel (chessel@biomserv.univ-lyon1.fr)
Date: Fri Sep 24 2004 - 19:00:41 MEST

  • Next message: Stephane DRAY: "RV or not RV ?"

    At 13:17 24/09/2004 +0200, Marie Sémon wrote:
    >Bonjour,
    >
    >Est-ce qu'il y a une fonction dans ADE4 sous R permettant de produire directement
    >les 9 graphes des valeurs propres associés aux 9 analyses de l'analyse des
    >correspondances internes? (cf Figure 1 de l'article de Lobry et Chessel, 2003
    >J. Appl. Genet.)

    Non, mais en voilà une.

    ww est un objet de la classe witwit créé par witwit.coa
    exemple :
    data(ardeche)
    coa1 <- dudi.coa(ardeche$tab, scann = FALSE, nf = 4)
    ww <- witwit.coa(coa1, ardeche$row.blocks, ardeche$col.blocks, scann = FALSE)
    witwitsepan(ww,c(4,6))

    Détails dans http://pbil.univ-lyon1.fr/R/querep/qr7.pdf

    witwitsepan <- function (ww, mfrow = NULL, csub = 2, plot = TRUE) {
        if (!inherits(ww, "witwit")) stop ("witwit object expected")
        appel <- as.list(ww$call)
        rowblo <- eval(appel[[3]], sys.frame(0))
        colblo <- eval(appel[[4]], sys.frame(0))
        anal <- eval(appel[[2]], sys.frame(0))
        tab <- eval(as.list(anal$call)[[2]], sys.frame(0))

        rowfac=as.factor(rep(1:length(rowblo),rowblo))
        if (is.null(names(rowblo))) names(rowblo) <- as.character(1:length(rowblo))
        levels(rowfac)=names(rowblo)
        
        colfac=as.factor(rep(1:length(colblo),colblo))
        if (is.null(names(colblo))) names(colblo) <- as.character(1:length(colblo))
        levels(colfac)=names(colblo)
        
        listblocrow = split(tab,rowfac)
        listbloc = NULL
        lapply(listblocrow, function(x)
             listbloc <<- c(listbloc,split(as.data.frame(t(x)),col.fac)))

        fun1 <- function(x) {
            x <- data.frame(x)
            if (nrow(x) <2) return (NULL)
            if (ncol(x) <2) return (NULL)
            sumlig <- apply(x,1,sum)
            if (sum(sumlig>0)<2) return (NULL)
            sumcol <- apply(x,2,sum)
            if (sum(sumcol>0)<2) return (NULL)
            return(dudi.coa(x,scann=F)$eig)
        }
        
        names(listbloc) <- t(outer(names(rowblo),names(colblo),function(x,y) paste(x,y,sep="/")))

        result <- lapply(listbloc,fun1)
        if (!plot) return(result)
        
        opar <- par(ask = par("ask"), mfrow = par("mfrow"), mar = par("mar"))
        on.exit(par(opar))
        par(mar = c(0.6, 2.6, 0.6, 0.6))
        nbloc <- length(result)
        if (is.null(mfrow))
            mfrow <- n2mfrow(nbloc)
        par(mfrow = mfrow)
        if (nbloc > prod(mfrow))
            par(ask = TRUE)
        neig <- max(unlist(lapply(result,length)))
        maxeig <- max(unlist(result))
        for (ianal in 1:nbloc) {
            w <- result[[ianal]]
            su0 <- names(result)[ianal]
            scatterutil.eigen(w, xmax = neig, ymax = maxeig, wsel = 0,
                sub = su0, csub = csub, possub = "topright")
        }
        return(invisible(result))

    }

    Daniel Chessel - chessel@biomserv.univ-lyon1.fr



    This archive was generated by hypermail 2b30 : Fri Sep 24 2004 - 19:02:02 MEST