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]]
_______________________________________________
R-sig-ecology mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
