Richness estimates confidence intervals very high.

[mestaki]

Hey folks, I've got a simple nested 2x3 design of mouse gut tissue which I have denoised with dada2 (with pooled=True) and then phyloseq+breakawy to estimate richness:

plot(ps.richness.breakaway, ps, color="group") 
# `group` here is the merged version 2x3=6 factors. Couldn't figure out plotting the 2x3 design.

First question, I assumed the error bars are 95% CI but I actually couldn't find this in the packaged DivNet vignette. Breakaway didn't come packaged with a vignette. Is this correct? Would be good to know to describe in a publication. Second, as you can see the CI are rather highly variable, in some instances reaching almost 50k which biologically doesn't make sense to me..is this typical? I thought those instances may be the result of low sequencing depths, but I couldn't confirm this pattern given the sequencing depth of the samples.

 DM01  DM02  DM03  DM04  DM05  DM06  DM07  DM08  DM09  DM10  DM11  DM12  DM13  DM14  DM15  DM16 
 6715   145 11006  3560 11668  3981  6603 10553  2738 11676 20038 13292   537  4628 46729 23484 
 DM18  DM19  DM20  DM21  DM22  DM23  DM24  DM26  DM27  DM28  DM29  DM30  DM32  DM33  DM37  DM38 
18304  5885 21195 19050 14348 13088 19430 19489  9656 18852  9507  6645 12403  7445 19437 24984 
 DM39  DM40  DM42  DM43  DM44  DM45  DM46  DM47  DM48  DM49  DM53  DM55 
 7130 13516 15407 11437 42219 26092  6607  7494 12115  7186  4701  2766 ```

Any thoughts on how to explain this and more importantly if these values are to be trusted when this high?

Thanks!

[adw96]

Hey Bod @mestaki ! So sorry again for the delay -- but it's updates time for breakaway 😺 so expect some answers to your questions and issues soon.

Indeed these are 95% confidence intervals. Because breakaway is a wrapper for multiple species richness estimates, the level of the set is actually listed within those functions. e.g. In the documentation for kemp is says ci: An asymmetric 95\% confidence interval for diversity. Sorry that's not easier to find.

To answer your second question, the upper endpoints of species richness estimates are indeed typically high. Here are couple of reasons for this

• Species richness estimation is an extrapolation problem, where we extrapolate off the range of the data (in this case, the frequency counts). Therefore there is a lot of uncertainty, and a high variance/standard error in the estimates.
• The species richness estimates that breakaway implements are all asymptotically normally distributed, but in finite samples we typically observe a distribution more like a log normal. For this reason, the default confidence intervals that we show are asymmetric and based off the log normal distribution. Therefore, I would argue that although it's unappealing to have so much error, this is a correct mathematical description of how uncertain each individual species richness estimate is.

There are a couple of reasons why this isn't for cause for concern:

• Each individual samples' species richness doesn't matter too much. If you want to do inference on whether diversity changes, you're going to be using all of the samples together. The overall uncertainty in the inference problem is much less than the uncertainty in the individual estimation problems.
• If you use betta to do inference, we don't use the asymmetric intervals -- we use symmetric intervals, which are considerably narrower.

To answer your final questions -- yes this is typical, it is reliable, and it is correct 😸 That's why you collected lots of samples!

Let me know if you have any questions about how to run betta for inference on diversity.

Amy