[R-sig-eco] Simper analysis with Morisita-Horn

[Eduard Szöcs]

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
>