Simper analysis with Morisita-Horn
Hai Vesna, your code looks reasonable. However it will fail, when you try running the permutations: 'Error: object 'gpa' not found' This comes from your second inserted code block (should be 'ga' ect). However in the current version it is not intended that the 'normal' user knows about the permutations. And one must modify the function (change the second line) in order to get the permutations running. So you can skip your second code block when you don't want to run the permutations (and leave the function as is). Also have a look at the development-version of vegan at github: https://github.com/jarioksa/vegan/blob/master/R/simper.R for a slightly changed version. Hope this helps, Eduard
On 12/03/2012 10:18 AM, Vesna Gagic wrote:
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 R-sig-ecology at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology