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
}

Regards,
Vesna
	[[alternative HTML version deleted]]

