Implanting first pass on no w/in plot agg about mobr HOT 7 CLOSED

Here is a better proof of concept and a function to carry out the computation of the individual swapped community by site matrix. The code also performs a quick test at the bottom that checks that species sums have not changed

S = 6
N = 100
nplots = 10
rand_comm = matrix(rpois(S*nplots, 1),
                   ncol=S, nrow=nplots)

intra_sp_swap = function(comm, groups) {
  group_levels = unique(groups)
  comm_group_noagg = matrix(NA, ncol=ncol(comm), nrow=nrow(comm))
  for(i in seq_along(group_levels)) {
    row_indices = groups == group_levels[i]
    comm_group = comm[row_indices, ]
    sp_sums = colSums(comm_group)
    tmp_comm = sapply(sp_sums, function(x) 
                      table(c(sample(1:nrow(comm_group), x,
                                    replace=T),
                              1:nrow(comm_group))) - 1)
    comm_group_noagg[row_indices, ] = tmp_comm
  }  
  comm_group_noagg
}

# test species sums for equality
replicate(20, all.equal(colSums(rand_comm), 
          colSums(intra_sp_swap(rand_comm, rep(1:2, each=5)))))

from mobr.

here is the best form of the two functions needed for incoperating this into our code base. There are two functions: 1) computes shuffled comm matrix, 2) computes averaged sample based rarefaction curve. These two functions can just be called before normal sample based rarefaction for both the observed and null perms. it is very fast luckily

intra_sp_swap = function(comm, groups) {
    group_levels = unique(groups)
    comm_group_noagg = matrix(NA, ncol=ncol(comm), nrow=nrow(comm))
    for(i in seq_along(group_levels)) {
        row_indices = groups == group_levels[i]
        comm_group = comm[row_indices, ]
        sp_sums = colSums(comm_group)
        tmp_comm = sapply(sp_sums, function(x) 
                          table(c(sample(1:nrow(comm_group), x,
                                        replace=T),
                                 1:nrow(comm_group))) - 1)
        comm_group_noagg[row_indices, ] = tmp_comm
    }  
    comm_group_noagg
}

avg_swap_rare = function(comm, groups, nperm=1000, effort=NULL){
    S = replicate(nperm, 
                  rarefaction(intra_sp_swap(comm, groups),
                              'samp', effort))
    Savg = apply(S, 1, mean)
    Savg
}

from mobr.

Bummer minor bug in code not quite as fast as I thought maybe we could turn down the default 1000 replicates to 500 or so.

from mobr.

I should note though that the functions I provide above appear to be working

from mobr.

Looks great, thanks!

from mobr.

Forget about what I said this morning - let's keep the functions as they are for now. It's actually easier to fit them into our original framework this way.

from mobr.

Yea with my functions you can supply them just a single group

Dan

On Jul 15, 2016, at 11:07 AM, Xiao Xiao [email protected] wrote:

Forget about what I said this morning - let's keep the functions as they are for now. It's actually easier to fit them into our original framework this way.

—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub, or mute the thread.

from mobr.