Merge branch 'morisita'

Author: jarioksa

R/vegdist.R

@@ -55,6 +55,9 @@
     if (method %in% c(7, 13, 15) && !isTRUE(all.equal(x, round(x))))
         warning("results may be meaningless with non-integer data in method ",
                 dQuote(inm))
+    if (method %in% 7 && round(max(x)) == 1) # morisita
+        warning("results may be meaningless with largest integer 1 in method ",
+                dQuote(inm))
     d <- .Call(do_vegdist, x, as.integer(method), PACKAGE = "vegan")
     d[d < ZAP] <- 0
     if (any(is.na(d)))

man/vegdist.Rd

@@ -226,7 +226,10 @@
   in the species present at site \eqn{j} that occur with only one
   individual in site \eqn{k} (Chao et al. 2005).
 
-  Morisita index can be used with genuine count data (integers) only. Its
+  Morisita index can be used with genuine count data (integers) only. It
+  is based on the idea of resampling without replacement and should not
+  be used with presence/absence data, but may give meaningless results
+  if compared sampling units (rows) have largest integer 1. Its
   Horn--Morisita variant is able to handle any abundance data.
 
   Mahalanobis distances are Euclidean distances of a matrix where

src/vegdist.c

@@ -352,8 +352,13 @@ static double veg_kulczynski(double *x, int nr, int nc, int i1, int i2)
      return dist;
 }
 
-/* Morisita index.  Can only be used with integer data, and may still
- * fail with unfortunate pairs of species occurring only once.
+/* Morisita index.  Can only be used with non-negative integer data
+ * (counts) and is hardly meaningful for presence/absence data or
+ * comparisons where largest integer is 1.  Function involves terms
+ * tlam1, tlam2 which for x=1 are x*(x-1) = 0 and if both are 0 this
+ * will give sim/0 which for sim > 0 is -Inf (made to 0), and for sim
+ * == 0 is NaN (made to 1). While this technically works, this is not
+ * very meaningful.
  */
 
 static double veg_morisita(double *x, int nr, int nc, int i1, int i2)
@@ -380,8 +385,11 @@ static double veg_morisita(double *x, int nr, int nc, int i1, int i2)
 	  i2 += nr;
      }
      if (count==0) return NA_REAL;
-     dist = 1 - 2*sim/(tlam1/t1/(t1-1) + tlam2/t2/(t2-1))/t1/t2;
-     if (dist < 0)
+     if (sim <= 0)
+	 dist = 1;
+     else
+	 dist = 1 - 2*sim/(tlam1/t1/(t1-1) + tlam2/t2/(t2-1))/t1/t2;
+     if (dist < 0) /* round-off error or sim > 0 && tlam1+tlam2 == 0 */
 	  dist = 0;
      return dist;
 }