[Vegan-commits] r2550 - pkg/vegan/R

[noreply at r-forge.r-project.org]

Author: jarioksa
Date: 2013-07-10 13:48:39 +0200 (Wed, 10 Jul 2013)
New Revision: 2550

Added:
   pkg/vegan/R/disp_weight.R
Log:
dispersion weighting - first commit

Added: pkg/vegan/R/disp_weight.R
===================================================================
--- pkg/vegan/R/disp_weight.R	                        (rev 0)
+++ pkg/vegan/R/disp_weight.R	2013-07-10 11:48:39 UTC (rev 2550)
@@ -0,0 +1,59 @@
+`disp_weight` <-
+    function(comm, group, nperm)
+{
+    # number of replicates per group
+    nrep <- tabulate(group) 
+    # workhorse
+    dfun <- function(comm, group, nperm, nrep) {
+        ## Calc Dispersion
+        # group means
+        means <-  tapply(comm, group, mean)
+        # omit groups with mean == 0
+        if(any(means == 0)) {
+            comm <- comm[means[group] != 0]
+            group <- group[means[group] != 0, drop = TRUE]
+            nrep <- nrep[means != 0]
+            means <- means[means != 0]
+        }
+        # group variances
+        vars <-  tapply(comm, group, var)
+        # dispersion
+        d <- vars / means
+        # average dispersion 
+        d_hat <- sum(d * (nrep - 1)) / sum(nrep - 1)
+        
+        ## Test
+        # original chisq values
+        chi_o <- sum((comm - means[group])^2 / means[group])
+        
+        # permutations
+        # calculate chisq
+        pfun <- function(comm, group){
+            means <-  tapply(comm, group, mean)
+            if(any(means == 0)) {
+                comm <- comm[means[group] != 0]
+                group <- group[means[group] != 0, drop = TRUE]
+                means <- means[means != 0]
+            }
+            chi <- sum((comm - means[group])^2 / means[group])
+            return(chi)
+        }
+        # sum of individuals per group
+        sums <- tapply(comm, group, sum)
+        # realocate randomly individuals to replications
+        perms <- vector('list', length(sums))
+        for(i in seq_along(sums)) {
+            perms[[i]] <- rmultinom(n = nperm, size = sums[i], prob = rep(1, nrep[i]) / nrep[i])
+        }
+        perms <- t(do.call(rbind, perms))
+        chi_p <- apply(perms, 1, pfun, sort(group))
+        p <- (sum(chi_p >= chi_o) + 1) / (nperm + 1)
+        out <- list(D = d_hat, p = p, weights = ifelse(p < 0.05, 1/d_hat, 1))
+        return(out)
+    }
+    # apply workhorse to every species
+    out <- apply(comm, 2, dfun, group, nperm, nrep)
+    # format output
+    out <- do.call(rbind.data.frame, out)
+    return(out)
+}
\ No newline at end of file