[R] Simper analysis with Morisita-Horn

[VG]

Dear ecology fellows,

I tried to implement Morisita-Horn distance (instead of Bray that is in the
current version)  in the code for the Simper analysis in vegan. I would be
very grateful if someone can check if the code  is right.

function (comm, group, ...) 
{
    if (any(rowSums(comm, na.rm = TRUE) == 0)) 
        warning("you have empty rows: results may be meaningless")
    permutations <- 0
    trace <- FALSE
    comm <- as.matrix(comm)
    comp <- t(combn(unique(as.character(group)), 2))
    outlist <- NULL
    P <- ncol(comm)
    nobs <- nrow(comm)
    if (length(permutations) == 1) {
        perm <- shuffleSet(nobs, permutations, ...)
    }
    else {
        perm <- permutations
    }
    if (ncol(perm) != nobs) 
        stop(gettextf("'permutations' have %d columns, but data have %d
rows", 
            ncol(perm), nobs))
    nperm <- nrow(perm)
    if (nperm > 0) 
        perm.contr <- matrix(nrow = P, ncol = nperm)
    for (i in 1:nrow(comp)) {
        group.a <- comm[group == comp[i, 1], ]
        group.b <- comm[group == comp[i, 2], ]
        n.a <- nrow(group.a)
        n.b <- nrow(group.b)
        contr <- matrix(ncol = P, nrow = n.a * n.b)
        for (j in 1:n.b) {
            for (k in 1:n.a) {

 ########### Morisita-Horn 
               
    aN <- sum(group.a[k, ])
    bN <- sum(group.b[j, ])
    da <- sum(group.a[k, ]^2)/aN^2
    db <- sum(group.b[j, ]^2)/bN^2
    top <- (group.a[k, ] * group.b[j, ])
    contr[(j - 1) * n.a + k, ] <- 2 * top/((da + db) * aN * bN)
#############
            }
        }
        average <- colMeans(contr)
        if (nperm > 0) {
            if (trace) 
                cat("Permuting", paste(comp[i, 1], comp[i, 2], 
                  sep = "_"), "\n")
            contrp <- matrix(ncol = P, nrow = n.a * n.b)
            for (p in 1:nperm) {
                groupp <- group[perm[p, ]]
                ga <- comm[groupp == comp[i, 1], ]
                gb <- comm[groupp == comp[i, 2], ]
                for (j in 1:n.b) {
                  for (k in 1:n.a) {
 
 ########### Morisita-Horn 

    aNp <- sum(gpa[k, ])
    bNp <- sum(gpb[j, ])
    dap <- sum(gpa[k, ]^2)/aNp^2
    dbp <- sum(gpb[j, ]^2)/bNp^2
    topp <- (gpa[k, ] * gpb[j, ])
    contrp[(j - 1) * n.a + k, ] <- 2 * topp/((dap + dbp) * aNp * bNp)
#############
                  }
                }
                perm.contr[, p] <- colMeans(contrp)
            }
            p <- (apply(apply(perm.contr, 2, function(x) x >= 
                average), 1, sum) + 1)/(nperm + 1)
        }
        else {
            p <- NULL
        }
        overall <- sum(average)
        sdi <- apply(contr, 2, sd)
        ratio <- average/sdi
        ava <- colMeans(group.a)
        avb <- colMeans(group.b)
        ord <- order(average, decreasing = TRUE)
        cusum <- cumsum(average[ord]/overall)
        out <- list(species = colnames(comm), average = average, 
            overall = overall, sd = sdi, ratio = ratio, ava = ava, 
            avb = avb, ord = ord, cusum = cusum, p = p)
        outlist[[paste(comp[i, 1], "_", comp[i, 2], sep = "")]] <- out
    }
    attr(outlist, "permutations") <- nperm
    class(outlist) <- "simper"
    outlist
}




--
View this message in context: http://r.789695.n4.nabble.com/Simper-analysis-with-Morisita-Horn-tp4651283.html
Sent from the R help mailing list archive at Nabble.com.