Further problems with ChaoShared

[PhilippeJeanneret]

Dear Anne,
Further testing of 'ChaoShared' showed errors difficult to understand:
test <- data.frame(v1 = c(1,1,0,0), v2 = c(1,0,1,1))
ChaoShared(test, datatype="abundance", se=TRUE, nboot=200, conf=0.95)
gives outputs but this:
data.frame(v1 = c(1,1,0,1), v2 = c(1,0,1,0))
results in "Error in rmultinom(B, n2, p2) : NA in probability vector
Called from: rmultinom(B, n2, p2)"
Quite strange, isn't it? In both cases, both communities have 1 species shared (the first one), one community has 1 exclusive species and the other community has 2 exclusive species.
But more remarkably, by reversing v1 and v2, 'ChaoShared' gives an output:
test <- data.frame(v1 = c(1,0,1,1), v2 = c(1,1,0,0))

I hope you can help me to understand, respectively resolve the problem.
Thanks in advance for your answer.
Best regards,
Philippe

[PhilippeJeanneret]

Dear Anne,
Thank you for considering my mail.
I know that no statistics can be derived from communities with three or four species.
However, I was not confident in the code "ChaoShared" when seeing that the ranking of species occurences in two communities did have an impact on the result, see my previous comment.
I further checked with larger communities and was happy to see reasonable results.
The question remains about how many species must the communities have to give reasonable results in the sense that the ranking of the species has no impact on the functionning of "ChaoShared"?
Furthermore, I still remain surprised to not have "ChaoShared" working when one community is a subset of the second one:
test <- data.frame(v1 = c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1), v2 = c(0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1))
ChaoShared(test, datatype="abundance", se=TRUE, nboot=200, conf=0.95)
Any comment on that?
Thank you in advance.
Philippe

Error in rmultinom(B, n2, p2) : NA in probability vector
Called from: rmultinom(B, n2, p2)
Browse[1]> Q

[dsiddons]

The problem seems to stem from "SortDataFun"

specifically this line of code

sit1 <- sample(sit1[1:(D1-D12)], replace=F)

When the result of D1-D12 is 1, it gives three numbers (a known problem/result of doing this with a single value) rather than the result of "3".

try it:

x1 = c(1,0,1,0), x2 = c(1,1,0,1)
sit1 <- order(x1 > 0 & x2 == 0, decreasing=T)
set.seed(123)
sit1 <- sample(sit1[1:(D1-D12)], replace=F)  # random sample
special1 <- cbind(x1[sit1[1:(D1-D12)]], x2[sit1[1:(D1-D12)]])

the result for special one is "1,1" and not "1,0" as it should be.

a quick fix would be an if/else on the value of D1-D12 (and I suppose for D2-D12 for sit2). e.g.

if (D1-D12 == 1){
    special1 <- cbind(x1[sit1[1]], x2[sit1[1]])  
    } else {
    set.seed(123)
    sit1 <- sample(sit1[1:(D1-D12)], replace=F)  # random sample
    special1 <- cbind(x1[sit1[1:(D1-D12)]], x2[sit1[1:(D1-D12)]])
    }