p-value at gamma scale

[rueuntal]

Previous discussions at #62 .

[FelixMay]

Hi @rueuntal , well, somehow I was confused and thought two much about confidence intervals for the single values instead of a confidence interval for the difference. Thanks for your detailed explanation. I now think that makes sense, but I will discuss it with Jon on Tuesday and then start implementing the p-value at the gamma scale.

[rueuntal]

Hi @FelixMay - sounds great, thanks!

[FelixMay]

Sorry for the delay, but now I started implementing the significance test at the gamma scale. As a first attempt I calculate bootstrapped confidence intervals for each index and each group. I did this by using resampling with replacement within each group. The drawback is that it does not work for species richness, because due to the nature of richness the richness of any resample WITH replacement is lower than in the real data. Thus, the observed value is larger than the "upper" confidence bound. As second step I plan to implement a confidence interval for the difference between control and treatment as suggested by @rueuntal above. This is straightforward when there are just two groups, but not when there are three. So the user need to define a reference group that is compared to everything else I guess? At the moment the reference group is set in the function get_delta_stats. One option is that I add a corresponding argument also in get_mob_stats. However, as this introduces redundancy, I wonder of the reference group should already be defined in make_mob_in? What do you think @rueuntal and @dmcglinn ?

[FelixMay]

In addition I used the package plotrix for plotting the confidence intervals. If we decide to adopt this, the package needs to be added to the dependencies in the DESCRIPTION file.

[dmcglinn]

The PR #140 implements @FelixMay 's solution. Just wanted to note that to provide a link between these comments threads.

[dmcglinn]

Thanks for this contribution @FelixMay! The graphics look cool. I see what you mean about the interval for S only going down based on the resample.

Based on your code it looks like you carried out the bootstrap procedure by randomly drawing plots with replacement. Did you and Jon ever consider randomly drawing individuals with replacement as an alternative? Also I know in #134 you and @rueuntal discussed whether bootstrap or sample permutation was the better approach. In the end did you decide to code up the bootstrap because it provides a better visual sense of the uncertainty in each group with the CI's or was there a another reason?

With respect to the specification of a reference group. I would prefer you add it as an argument to get_mob_stats rather than make_mob_in despite the redundancy because I think user's should have the flexibility of using a single mob_in data object in different kinds of comparisons just by tweaking the argument to the analysis function rather than having to go and recreate the data input.

[FelixMay]

Hi @dmcglinn, alright I will add the reference group argument to get_mob_stats. So far I did not think about the resampling of individuals, but I will discuss it with Jon. With the definition of a reference group I can also implement the approach discussed with @rueuntal , which will provide a p-value for the difference between the reference group and the rest of the data, when I got it correctly. How to proceed? Will you first merge #140 and then I implement all the other changes or should I keep committing to #140? (sorry I am still not always confident what is the best way to do this stuff in GitHub)

[dmcglinn]

I'm working on a code review of #140 right now which I should finish in a few minutes. Once you get that and respond to my suggestions with comments and possibly new commits then we'll merge the PR. I'll fix the conflict when I merge. You'll want to do a pull after #140 is merged because it will differ from your local version due to the conflict in the DESCRIPTION file.

[rueuntal]

Hi @FelixMay @dmcglinn - thanks of working this out! Seems that you've already got a solution. But if my input is needed please let me know.

Felix - My apologies for the delay. I'll be working full-time in a data science workshop for another four weeks, and will probably only have a few hours every week to work on anything else. But I'll try my best to get back to you guys if needed.

[FelixMay]

Hi @rueuntal, @dmcglinn , @ngotelli - or dear stats wizards In the most recent mobr implemention there are now p-values for treatment comparisons at the plot and at the group level. However, these are sometimes zero when the observed test statistics are more extreme than any of the the permutation-based values (see mob_boxplots.R lines 123 and 457). Jon pointed to the issue that p-value of zero are strange, especially when they are based on a finite number of permutations.

So my first question: Is the implementation correct? For the F-test I tried to follow Nick's advice as far as I remember and for the group-level test I considered the Legendre $ Legendre Numerical Ecology book.

When everything is correct, do we somehow have to substitute the p = 0 values by something like p <= 1/npermutations ?

[ngotelli]

Hi @FelixMay

Yes, this is the correct way to handle a calculated p value of zero based on a randomization test: if the observed value is more extreme than any of the simulated values, we convert it to a tail probability, with p <= 1/(# of permutations)

Thanks!