The Power of Microbiome Studies: Some Considerations on Which Alpha and Beta Metrics to Use and How to Report Results

1For a few decades now, researchers have left culture-based methods and used molecular technologies, and more recently mostly sequencing-based approaches, to characterize microbial communities within a certain environment, referred to as the microbiome. In humans and animals, the microbiome has an important role in health and disease. For example, animals raised without or fewer microbes showed an underdeveloped immune system and are more susceptible to diseases. Microbiome studies have as goal to investigate, characterize, and understand the compositional and functional variability of microbiomes. The question “What is different between different groups of interest?” can be translated into a hypothesis-testing procedure.

2Hypothesis testing rests on the definition and choice of four parameters: (i) the effect size, i.e., the quantification of the outcome of interest (in the simple case, the difference between two groups); (ii) the sample size n, i.e., the number of samples (to be) collected; (iii) the power of tests 1 - β, i.e., the probability of the test of rejecting the null hypothesis when actually false; and (iv) the confidence level α, i.e., the probability of rejecting the null hypothesis when actually true.

3It is necessary to perform power analysis before performing experiments. This is well acknowledged in all fields of research; however, microbiome studies are challenged with conflicting results. Underpowering and the failure to correct for false positives are among the causes underlying the lack of reproducibility of many biological findings.

4The power of a test is linked to the probability β of accepting the null hypothesis when actually false (false-negative error or Type II error), and α describes the false-positive error or Type I error. Once acceptable error rates α (usually 0.05 or 0.01) and β (usually 0.2, although context-dependent) and the effect that one is interested to assess statistically are chosen, it is possible, at least in principle, to determine the optimal sample size, i.e., the number of samples that one needs to collect/analyze to obtain, with probability 1 −β, a statistically significant result with confidence α.

5Given the nature of microbiome data, it is possible to quantify differences between groups at two levels: the alpha (within-sample) and beta (between-sample) diversity (Figure 1). Alpha diversity metrics summarize the structure of a microbial community with respect to its richness (number of taxonomic groups), evenness (distribution of abundances of the groups), or both. Commonly used alpha metrics are phylogenetic diversity (PD), observed number of amplicon sequence variants (ASVs), Chao1, Simpson’s, and Shannon’s indices. Beta diversity metrics summarize which samples differ from one another by considering sequence abundances or considering only the presence–absence of sequences. Commonly used beta metrics are the Bray–Curtis (BC) dissimilarity, Jaccard, unweighted UniFrac (UF), and weighted UniFrac.

6The choice of the diversity metrics affects the subsequent statistical testing and, as a result, how, and to which extent, power analysis can be performed.

7With the use of an alpha diversity metric, a single diversity value is obtained for each sample containing measurements of m taxa; thus, the problem of assessing differences between two (or more) groups can be addressed with a univariate test, like t-test, ANOVA, or a nonparametric test. The use of a beta diversity metric implies that all samples are to be considered simultaneously, and several methods to compare groups of samples measured on m > 1 have been proposed like analysis of similarities (ANOSIM) and permutational multivariate ANOVA (PERMANOVA) to replace classical multivariate tests like the Hotelling T2 or multivariate ANOVA, which are in general not applicable. This happens because basic assumptions are not met, such as independence of the sample units, the multivariate normality of errors, homogeneity of variance–covariance matrices among the groups, or because the number of variables is larger than the number of samples, making it impossible to apply the test.

8While sample size and Type I and Type II errors are well-defined concepts, the definition of effect size depends on the outcome quantity of interest and how this quantity is mathematically defined. A fundamental step when performing power analysis is then to define the effect size: for a simple two-sample t-test, the effect size can be expressed as Cohen’s δ.

9where μ1 and μ2 are the population means of the two groups to be compared and σ2 is the pooled variance. Since μ1 and μ2 are population parameters that are inaccessible and on which we want to perform inference, an a priori estimation, or educated guess, is necessary. This can be accomplished by taking the sample means m1 and m2 and pooled variance s2 from a pilot study or existing data to obtain estimates of the population parameters. A critical aspect that is not sufficiently acknowledged is that the effect size from Equation (1) is sensitive to the particular diversity metric used to the point that sample size calculations can be severely affected.

10The aim of this study is to show how, and to what extent, different diversity metrics influence the sample size needed to assess the statistical significance of dissimilarities between different microbial communities. Both simulated and empirical 16S rRNA gene amplicon sequence datasets are used to perform retrospective calculations of the empirical power for microbiota studies. A broad selection of alpha and beta diversity metrics was used to compare the different microbiota datasets.

11This study generated insight into the sensitivity and bias of certain statistical methods used in microbial ecology on microbiota research outcomes. We conclude with some recommendations for the reporting of power analysis and sample size calculations for microbiome studies.

12To support the choice of alpha and beta diversity metrics to consider in our comparison, we performed a literature search on PubMed1 with a query: (microbiota [Title] OR microbiome [Title]) NOT Review [Publication Type] 2020/01:2020/02 [Date of Publication]). This strategy aimed to include a broad scope of microbiota studies. We limited our search to studies published in English with free full text.

13Richness is the number of taxa, most often defined as an operational taxonomic unit (OTU) or ASVs observed, where s is the number of observed taxa, calculated as

14PD is a phylogenetically weighted measure of richness. Although the name suggests diversity, it does not take into account the abundance of taxa. The PD is defined as the sum of the lengths of all those branches on the tree that span the members of the set, given the phylogenetic tree spanning s taxa:

15where s is the number of observed taxa and bi is the length of the ith branch in the tree. Index i runs on all branches.

16The Chao1 index is an abundance-based nonparametric estimator of taxa richness. It is defined as

17where s is the number of observed taxa, and F1 and F2 are the number of OTU/ASV with only one sequence (i.e., “singletons”) and two sequences (i.e., “doubletons”). This metric assumes the number of organisms identified for a taxon to follow a Poisson distribution. The definition rests on the concept that rare taxa bring most information about the number of missing taxa. This index gives more weight to the low-abundance taxa, and only the singletons and doubletons are used to estimate the number of missing taxa. This index is particularly useful for datasets skewed toward low-abundance taxa. However, singletons and doubletons are often removed from 16S rRNA amplicon sequence datasets because of the difficulty in robustly differentiating singleton errors from real singleton sequences.

18Shannon’s index H is an estimator of taxa diversity, combining richness and evenness). It is defined as

19where s is the number of OTU/ASV and pi is the proportion of the community represented by the ith OTU/ASV. Basically, this index is the entropy associated with a given sample and quantifies the uncertainty in predicting the taxa identity of an individual selected at random from the sample. Shannon’s index uses the relative abundances of different taxa; thus, diversity depends on both taxa richness and evenness with which organisms are distributed among the different taxa. This index places a greater weight on taxa richness.

20Simpson’s index D is an estimator of taxa diversity, combining richness and evenness. It is defined as

21where s is the number of OTU/ASV and pi is the proportion of the community represented by the ith OTU/ASV. This index considers taxa evenness more than taxa richness in its measurement; it indicates the taxa dominance and gives the probability of two individuals that belong to the same taxa being randomly chosen. It varies from 0 to 1, and the index increases as the diversity decreases.

22The BC index measures the compositional dissimilarity between the microbial communities of two samples i and j based on counts on each sample. It is defined as

23where Cij is the sum of the smallest values for only those taxa in common between the sample i and j, and Si and Sj are the total number of taxa counted in sample i and j, respectively. This index ranges between 0 (the two samples share all taxa) and 1 (the two samples do not share any taxa). It gives more weight to common taxa. The BC dissimilarity is computed pairwise between all samples.

24The Jaccard distance J between two samples i and j is defined as J = 1 - J(i, j), where J(i, j) is the Jaccard index, which is defined as

25which is the ratio between the number of members that are common between the two samples and the number of members that are distinct; it is a measure of similarity for the two communities and ranges between 0 (the communities are different) and 1 (the two communities are identical).

26UF and weighted UniFrac distances between two samples take into account the phylogenetic tree and thus phylogenetic distances between community members. In UF, the distance is calculated as the fraction of the branch length, and in weighted UniFrac, branch lengths are weighted by the relative abundance of sequences. The sum of unshared branch lengths is divided by the sum of all tree branch lengths, which results in the fraction of total unshared branch lengths that is defined as

27defined n as the total number of branches in the tree; bi as the length of branch i; Ai and Bi as the numbers of sequences that descend from branch i in communities A and B, respectively; and AT and BT as the total numbers of sequences in communities A and B, respectively. In order to control for unequal sampling effort, Ai and Bi are divided by AT and BT.

28This dataset contains 16S rRNA gene amplicon sequence data obtained from a broiler chicken experiment. The dataset is described in detail in; briefly, chickens were raised under three different housing conditions with the same medium-chain fatty acid (MCFA) feed intervention. Between those housing conditions, bird management was kept as similar as possible. At the hatchery, the chicks were randomly allocated to three different experimental facilities. Dataset A contains 70 broilers from a grow-out feed trial facility, dataset B contains 70 broilers raised at a floor stable, and dataset C contains 70 broilers raised in isolators. A feed intervention was used as a tool to create differences in cecal microbiota between broilers within the same housing condition.

29This dataset was obtained from the Human Microbiome Project (HMP) phase I. It contains 16S rRNA gene amplicon sequence data of 169 stool samples, 150 oral samples, 86 vaginal samples, and 69 skin samples. The bodyside microbiomes are all diverse in terms of community membership.

30In all datasets, ASVs were defined as unique sequences. All data were analyzed using NG-Tax. Taxonomy was assigned using the SILVA 128 16S rRNA gene reference. An overview of dataset characteristics (sample size, the number of ASVs, and mean values of different alpha and beta diversity metrics) is shown in Table 1.

31We simulated two different scenarios described by Simulated datasets 1 and 2. An overview of the characteristics of the simulated datasets is shown in Table 2.

32Simulated dataset 1 is built starting from 1,995 microbial features observed in 169 stool samples from the HMP, indicated as X1 in the following. Datasets (named X2 for the sake of simplicity) were created where 2, 5, 10, 25, 50, and 75% of bacterial features were randomly removed. This simulation generates data under an exist–non-exist binary scenario where bacterial features are either present or absent in X1 and X2.

33Simulated dataset 2 is a case–control scenario (X1 controls and X2 cases, with 1,995 features and 169 samples each) where one-fourth of bacterial features in X2 are 1, 2, 5, 10, 15, and 20% more abundant than in X1. Both simulated datasets have the same phylogenetic structure as the HMP dataset. This simulation thus generates data under a differential abundant scenario where bacterial features are present in the different compositions in X1 and X2.

34Differences between groups using alpha diversity as determined by using PD, richness (defined as observed), Chao1, Simpson, and Shannon were assessed using the Kruskal–Wallis test. A significance threshold α = 0.01 was used in all calculations.

35The differences between groups using beta diversity as determined by using the BC, Jaccard, UF, and weighted UniFrac were assessed using the PERMANOVA. PERMANOVA is a robust approach to compare groups of samples measured on m > 1 variables. It constructs ANOVA-like test statistics from a matrix of resemblances (distances, dissimilarities, or similarities) calculated among the sample units and assesses the significance of the observed differences using random permutations of observations among the groups. The null hypothesis H0 tested by PERMANOVA is that the centroids of the groups (in the space of the chosen resemblance measure) are the same for all groups. This test assumes that samples are exchangeable under the null hypothesis, are independent, and have similar multivariate dispersion. The PERMANOVA test statistic is a pseudo ANOVA F-ratio:

36where SSB is the total sum of squares of the (diss)similarities between groups, SSW is the total sum of squares of the (diss)similarities within groups, g is the number of groups, and n is the total number of samples.

37The significance of the F-statistics is calculated by means of permutations (k = 9,999). The distribution of F under the null hypothesis is generated by permuting g times the sample group labels and recalculating F on the permuted data. Significance is expressed as a p-value calculated as the fraction of permuted F-statistics, which are equal to or greater than the pseudo F-ratio observed on the original data.

38The experimental and simulated data were used to generate K random datasets of different sizes to take into account both data generation variability and the calculation of the empirical power. More specifically, from N1 × m and N2 × m data matrices X1 and X2 (either experimental or from Simulations 1 and 2), we randomly sampled with replacements K = 1,000 n1 × m and n2 × m datasets for different sample sizes n1 and n2. For the sake of simplicity, we consider n1 = n2 = n, and we varied n between 5 and 50 or 100 in steps of 5. We used K = 1,000 for the analysis of univariate analysis (alpha) and K = 100 for the analysis of multivariate analysis (beta), both simulated and experimental data. In total, more than 40,000 randomly generated datasets were analyzed.

39A statistical test to assess the difference between X1 and X2 (as quantified by any of the alpha and beta metrics) at significance level α = 0.01 under the assumption of the null hypothesis being false was applied on the K randomly generated datasets for different sample sizes n. The empirical power of the test is defined as the empirical probability EPr of H0 being rejected, calculated as

40where # indicates the number of times that H0 is rejected.

41All statistical analyses were performed in R version 4.0.2 (R Foundation for Statistical Computing, Austria; ), using the following packages: Phyloseq, Microbiome, and Vegan. PERMANOVA was performed using the adonis function for the Vegan package. Other power calculations were performed using the G*power software using the “Means: Difference between two independent means (two groups)” as Statistical test and “a priori” and “post hoc” option for the Type of Power analysis. Differentially abundant microbiota profiles were simulated with the microbiomeDASim R package using the gen_norm_microbiome function. The R scripts can be found on the Github page: https://github.com/mibwurrepo/KersSaccenti-Power.

42We begin with a motivational example to show how the choice of the diversity metrics affects the power of a microbiome study and how the same study may be underpowered if a different metric is used.

43Let us suppose we want to plan an experiment to assess whether gut and oral microbial communities are different. A very simple and basic study design would be to collect n1 = n2 gut and oral samples and compare the alpha diversity between the two conditions (gut vs. oral) using a two-sample Kruskal–Wallis t-test.

44We can base our estimation d of a very similar effect size δ on data from HMP (Table 1). Using four different alpha metrics and Equation (1), we obtained d = 1.27 (PD), d = 0.3621 (Shannon), d = 0.58 (Chao1), and d = 0 (Simpson). These values are markedly different: fixing the power to 0.8 (β = 0.2) and confidence α = 0.05 will lead to dramatically different required total sample size (Figure 2A). This clearly indicates that microbiome studies may be severely underpowered depending on which alpha metric was used to compare two (or more) groups.

45We also explored the achievable power by fixing the sample size (n = n1 + n2 = 50 + 50 = 100) and using different effect sizes (Figure 2B). Consistently with what is observed in Figure 2A, results vary strongly, providing a clear indication of the risk of underpowering when Shannon’s diversity is used.

46Note that with the use of beta diversity metrics, performing a priori power analysis becomes much more complicated. The classical tools for power analysis cannot be applied since the statistical tools are not parametric: solutions have been proposed in the literature; see, for instance.

47Our literature search returned 632 papers matching the search criteria. We selected randomly 100 papers, and of those, the materials and methods or full text were investigated to obtain an overview of the most frequently used alpha and beta diversity metrics and the sample sizes used. Of the 100 full-text papers, 92% of the papers contained alpha metrics, and 83% of the papers contained beta metrics.

48In 58% of the papers, more than one alpha metric was used. In 21% of the papers, more than one beta metric was used. An overview of the frequency of the different metrics showed that Shannon’s index and the BC dissimilarity are the most common metrics (Table 3). There was a wide variance in the used sample size, defined as the smallest number per group: 46% of the papers had a sample size of ≤ 10 samples, 34% of the papers used between 11 and 50 samples, 7% of the papers used between 51 and 100 samples, 10% of the papers used between 101 and 1,000 samples, and three papers used > 1,000 samples.

49This (small) literature offers an indication of what the most used alpha and beta metrics are for the analysis of microbiome data. The results are not surprising and agree with the author’s experience. We can note that none of the papers screened mentioned the use of Hill’s numbers, a mathematically unified family of diversity indices (differing among themselves only by a parameter) that incorporate species richness and species relative abundances. The use of Hill’s number has found consensus in ecology, and they have also been used for the analysis of microbiome data. However, their use seems to be not widespread, and their utility is not fully acknowledged: in this study, we will focus on the more commonly used metrics.

50As shown in the motivational example, the choice of a particular alpha (and beta) diversity metric determines the number of samples required to achieve a predetermined power. Based on this observation, we examined two simulated datasets using both alpha and beta diversity metrics to understand the relationship between the sample size, the observed power, and the diversity metrics, together with two experimental datasets (chicken and HMP datasets).

51As representatives of testing procedures using alpha and beta diversity measures to compare two groups, the Kruskal–Wallis test (for alpha metrics) and PERMANOVA (for beta metrics), respectively, were selected. The Kruskal–Wallis test is the nonparametric choice for comparing two groups when the normality assumption does not hold. When comparing two (or more) groups using beta diversity metrics, PERMANOVA and ANOSIM are popular choices. The two approaches are equally popular (359 hits on PubMed for PERMANOVA and 341 for ANOSIM); however, showed that while both approaches are sensitive to unbalanced designs and differences in variance within groups, PERMANOVA is a more robust approach: on this ground, we based our choice for PERMANOVA.

52While in this study we focus only on inferential approaches for the analysis of microbiome data, as carried on using univariate and multivariate tests, we should comment that sample size is also relevant for exploratory approaches like principal component analysis (PCA), principal coordinate analysis (PCoA), and multidimensional scaling (MDS). These approaches are not inferential (as long one does not consider the inferential setting for dimension assessment;, ), but the number of samples affects the reliability and stability of loadings and stability; thus, asking what the (minimal) sample size is to obtain stable and reproducible component loading estimations is relevant. Very little is known on the topic in classical PCA and factor analysis setting from a theoretical point of view, which cannot be directly extended to PCoA and MDS, which are the most commonly used approaches in microbiome analysis.

53For the simulated datasets, the effect size is known a priori, and it is expressed as the % of differentially abundant or present/absent microbial features (ASV) (Figure 3). The achievable power for Simulated dataset 1 is shown as a function of the sample size (n) for different percentages of present/absent ASV. If 2%, 5% of the ASVs are deleted from the dataset, none of the alpha diversity metrics was able to capture the difference between datasets X1 and X2, irrespective of the sample size used (Figures 3A,B). When more than 10% of ASVs were removed from dataset X2 (Figures 3C–E), all measures were somehow able to capture the difference, but the resulting actual power was very different. Overall, Chao1 and observed diversity allowed higher power with the lower sample size (are more sensitive to observe differences), especially in the medium range of differences (10–25%, Figures 3C,D), whereas differences are minimal for > 25%. Note that in contrast with the motivation example, here, the PD was not the metric resulting in the smallest sample size.

54The same approach was used across different beta diversity metrics (Figure 4). The Jaccard diversity metric was the most sensitive and the weighted UniFrac was the least sensitive to observe the differences in the presence/absence between the datasets (Figures 4B–F). When more than 10% of ASVs were removed from dataset X2, no difference between datasets was observed by the UniFrac metric (Figure 4C), while with 25% removed, the BC and UF showed a comparable power and sample size (Figure 4D). In this simulated dataset, the weighted UniFrac distance needed the highest sample size to observe the difference (Figures 4D–F).

55The achievable power for Simulated dataset 2 is also shown as a function of the sample size for different percentages of differentially abundant ASVs (Figure 5). If ≤ 5% of the ASVs were differentially abundant in dataset X2 as compared with X1, Simpson’s metric needed the lowest sample size (is most sensitive) to observe differences between the data (Figures 5A–C). However, if 10% of the ASVs were differentially abundant, the PD and Chao1 were more sensitive and Simpson’s and Shannon’s metrics less sensitive (Figure 5D). With 15% of the ASVs differentially abundant, no differences were observed with Simpson’s metrics (Figure 5E).

56The same approach was used across different beta diversity metrics (Figure 6). The BC distance was the most sensitive to observe differences, whereas UF needed the largest sample size. If 2% of the ASVs were differentially abundant, the power of the beta metrics was totally different; for example, a sample size of 15 would result in the power of 100 for the BC, 50 for weighted UniFrac, 40 for the Jaccard distance, and just 10 for UF (Figure 6B). However, if 10% of the ASVs were differentially abundant, all metrics would result in a power higher than 0.80 (Figures 6D–F).

57Shannon’s index was the most sensitive alpha metric and Chao1 and PD were the less sensitive metrics to observe a difference between the groups in dataset B (Figure 7A). In dataset B, Shannon’s index was also the most sensitive alpha metric, but Simpson’s index was the least sensitive metric (Figure 7B). UF distance was the most sensitive beta diversity metric to observe a difference between groups in dataset A (Figure 8A). The Jaccard distance was the only metric that showed that dataset C needed the smallest sample size, indicating that in dataset C, specific ASVs are differentially present between the groups (Figure 8). Weighted UniFrac was more sensitive than UF to observe a difference between the groups based on their microbial communities (Figure 8). In general, the alpha diversity measures were less sensitive to observe differences between the broilers than the beta diversity (Figures 7A,B, 8).

58Although no difference in alpha diversity was observed between broilers fed with or without MCFAs raised in housing condition 1, the average daily gain and the average daily feed intake were lower in MCFA broilers. Therefore, the difference only observed based on the beta diversity might already be biologically relevant and hence sufficient to draw conclusions in this case. For this dataset, we observed that Shannon is the most sensitive alpha diversity metric to observe differences between groups, resulting in the lowest needed sample size. The sensitivity of the beta diversity, however, was different per dataset. Based on this retrospective power calculation, two conclusions can be drawn on this study design. First, not enough chickens were sampled to observe a difference in the alpha diversity between broilers fed with or without MCFA raised in dataset A. Second, 15 chicken samples instead of 35 samples per group would have resulted in the same conclusion.

59The samples in this dataset were collected from different body sites and are known to have a very distinct origin and therefore expected to be different in microbial composition. The comparison between different body sites showed a wide variation in sample size across different alpha diversities (Figures 7C–F). The difference in sample size was small in the comparison between skin and oral microbiome samples, a total of 10 samples (threshold power 80, 1-β) (Figure 7C). In the skin vs. gut microbiome samples, Simpson’s and Shannon’s alpha diversities did not differ, and the PD was the most sensitive to observe differences (Figure 7D). In contrast, when comparing the gut vs. the oral microbiome, the PD was the least sensitive to observe differences, whereas Shannon’s and Simpson’s metrics were different between the gut and oral samples (Figure 7E). In the skin vs. vaginal microbiome comparison, Simpson’s and Shannon’s alpha diversities were more sensitive than the observed/Chao1 and PD (Figure 7F). Based on the different beta diversity metrics, all comparisons between different body sites supported significant differences even when just five samples were compared (data not shown), due to the large difference between communities (Supplementary Figure 1). Therefore, the retrospective power calculations were not informative for this dataset.

60Figure 9A shows the distribution of the sample size (n1) of the datasets that were analyzed, using the Chao1 diversity measure in 28 of the 100 papers considered in the literature review. The distribution is highly skewed toward 0 with a median of 39 samples per group and a mode of 8 samples. Removing the two outlying studies with > 300 samples resulted in a median of 23 samples.

61Even considering that Chao1 was one of the best-performing measures in both simulated and experimental datasets, these numbers appear to be worryingly low: on experimental data, which have a complicated structure that is impossible to replicate in simulations, it is rarely possible to attain a power of 80% with less than 40 samples per group. Similar considerations hold when PERMANOVA is applied (see Figure 9B), with a median group size of 22 and mode 3.

62One of the studies we examined in the literary review reported that sample size and power analysis were performed: “Sample sizes were chosen on the basis of pilot experiments and on our experience with similar experiments.” This is commendable, but we believe that the way forward is to employ and report in full a standardized summary of sample size calculations performed. The software G*power generates a Protocol for power analysis. For instance, for a two-group comparison with a Mann–Whitney/Kruskal–Wallis test, a possible (modified) reporting is given in Table 4.

63Together with this, information should be provided on how effect size was determined, i.e., which pilot data were used and how the effect size was calculated.

64A similar reporting protocol could be devised if simulations are used in a PERMANOVA setting (Table 5). Since simulation and/or pilot data must be used in this case, details on the simulations or pilot data should be reported. For instance, using Chicken data 1 as the pilot, one could report the following protocol, taking 100 resamplings of size 6 to calculate the achievable power:

65The aim of this study was to assess how, and to what extent, different diversity metrics and compositions of the microbiota influence the needed sample size to observe statistically significant dissimilar groups. Based on our literature survey, we observed that Shannon’s and Bray–Curtis metrics are the most published metrics. This might be because they are often the most sensitive metrics to observe differences between groups, resulting in a lower sample size. Our results are in line with those of a previous literature that showed that the choice of distance metric may significantly influence the observed results.

66A well-known phenomenon that can hamper progress in every research field concerns publication biases in reporting mainly positive findings. In microbiota research, this might even occur rather unintentionally, by using certain alpha and beta diversity metrics, but it might also be that researchers selectively report only results for the metric that shows significance even when other metrics had been assessed during the analyses. Our results lead to the speculation that many microbiome studies may be underpowered or, conversely, only reporting evidence of very large effects that can be assessed to be statistically significant also with a small sample size. However, since effect size and test statistics are not reported, it is impossible to judge the quality of the results. This also hampers the use of published studies as pilot studies to perform power analysis and sample size calculations, as long as data are not de novo reanalyzed.

67None of the 100 microbiome studies that we have considered reported the effect size. A collaborative project aiming to investigate the reproducibility of 100 high-profile psychological studies reported that the average effect size observed in the replication studies was approximately half the magnitude of those given in the original studies, leading to a replication success of only 36%. The lack of reported effects makes it impossible to analyze retrospectively microbiome studies and to perform a meta-analysis and, more importantly, makes it impossible to check the consistency of the statistical analysis or detect errors.

68On the basis of this, reporting of effects and test statistics should be made compulsory in microbiome studies. For the highly used Kruskal–Wallis test, the H test statistic is given by (in absence of ties)

69where N is the total number of samples, ni is the number of samples in group i, and Ri is the average rank of observations in the ith group. Note that H (or alternative formulas) obtainable from most software packages.

70For the Kruskal–Wallis, the most common effect is the η2, which is defined as

71where H is the value obtained in the Kruskal–Wallis test and k is the number of groups. For instance, for the comparison of the two feedings (Feed A and B; see Table 1) from the chicken dataset using the observed alpha diversity, one could report the following:

72Feed A (n1 = 35) and Feed B (n2 = 35) samples were compared with Kruskal–Wallis test using the Chao1 metric: H (df) = H (1) = 14.68, p-value 0.0001, η2 = 0.66, δ = 0.58,

73where df indicates the degrees of freedom.

74Note here that for a two-group comparison, the Kruskal–Wallis test is equivalent to the Wilcoxon–Mann–Whitney (WMW) test. For the WMW test, Cohen’s δ effect size definition (Equation 1) also applies. This greatly simplifies power analysis and sample size calculations: we advise to also report δ when two groups are considered.

75In addition, performing power analysis for a Kruskal–Wallis is not a simple matter and requires the use of a rather advanced statistical machinery; for instance, such calculations are not included in G*power, which is the most complete software for power analysis. The Kruskal–Wallis is the nonparametric counterpart of one-way ANOVA and as such is used in situations where there are more than two groups. However, whereas power analysis and sample size calculation for a one-way ANOVA with more than two groups are “easily” accessible within R or other software packages, this is not the case in the Kruskal–Wallis testing. We could locate an R package “MultNonParam” that performs power analysis for the Kruskal–Wallis test with more than two groups; however, it requires the specification of the offsets for the various populations, under the alternative hypothesis. Relating such tools to determine the effect size observed in microbiome data is a matter we believe to be worthy of exploration and brings us back to the problem that statistics and effect size are not easily available for microbiome studies.

76The principle of reporting the effect size should also apply when testing is performed using beta diversity metrics, in which case the PERMANOVA pseudo F-statistics (see Equation 10) and the effect size should be reported. Typical effect measures in ANOVA are Choen’s f2, η2, and ϖ2:

77Feed A (n1 = 35) and Feed B (n2 = 35) samples were compared with PERMANOVA test using the Unweighted Unifrac metric: F(dfB, dfW) = F(1, 68) = 6.27, p-value = 0.0001, f2 = 0.092, 1,000 permutations,

78where dfB and dfW indicate the between-groups and within-groups degrees of freedom, respectively. These notations follow the guidelines of the American Psychological Association, which provides standardized formats for the reporting of statistical analysis for statistical procedures.

79For PERMANOVA, the matter complicates considerably: to estimate statistical power and calculate sample size, one must quantify the expected within-group variance and the effect to be expected when comparing two or more groups. A package like micropower in principle allows estimation of PERMANOVA effects quantified by the ϖ2 value (limited to the UniFrac measure): unfortunately, to the best of our knowledge, the package seems not to be maintained and lacks a proper manual. The original paper presents a table with effects calculated from different studies that could be used as a guide; however, this metric is not standard. It can be calculated from the PERMANOVA table as

80For the chicken data, we observed ϖ2 values in the range 0.04–0.1, depending on the beta metrics, and these are consistent with those reported in Table 1 from.

81A more common effect measure is Cohen’s f2-value; i.e., the between-group to within-group ratio can be easily obtained by the ANOVA table provided by software like the R package Vegan by taking the ratio between the Treatment sum of squares and the Residual sum of squares. The f2-value is used for power calculation in the ANOVA setting; however, it should not be used to perform sample size calculation for PERMANOVA, not even to obtain a rough indication, since the corresponding F-statistics do not follow an F distribution. For instance, when comparing Feed A and B from the first chicken data with PERMANOVA, we can derive a Cohen’s f2 = 0.38; if this value is used to perform power calculation for an ANOVA with two groups with power 80% at α = 0.01, we obtain that 42 samples per group are needed. However, comparing the results in Figure 8, we see that a 100% power can be obtained with 25 samples per group, regardless of which measure is used: this is a clear indication that power analysis for PERMANOVA can be obtained only by means of simulations. In this light, Figure 8 can be viewed as a priori power calculations using pilot data.

82Furthermore, there is a convention, more or less widely accepted, of classifying Cohen’s δ effect, which we have used in the power calculation for the Kruskal–Wallis, into trivial (δ < 0.2), small (δ = 0.2), medium (δ = 0.5), and large (δ > 0.8). However, this classification is based on what is observed in psychology and does not apply automatically to other fields of research. In microbiome studies, the effects may be in the same order of magnitude that may be considered large or very large using the standard convention.

83We have seen that different alpha and beta diversity metrics lead to different study power: on the basis of this observation, one could be naturally tempted to try all possible metrics until one or more are found that give a statistically significant test result, i.e., p-value < α. This way of proceeding is one of the many forms of the so-called p-value hacking (p-hacking). p-Hacking (also called data dredging, significance chasing, significance questing, or selective inference) is the improper use of data (like adding or removing observations) or statistical procedures (like applying many different tests) until a configuration is found that produces a statistically significant result at the desired confidence level. p-Hacking is an illegitimate practice that promotes unreproducible results, polluting literature and adding to publication bias.

84To this end, in our opinion, the only way to protect ourselves from (the temptation of) p-hacking would be to publish, and we stress here the word publish, a statistical plan before experiments are initiated: this practice is customary for clinical trials where a statistical plan describing the endpoints and the corresponding statistical analyses must be disclosed before the start of the study and must be adhered to if results are going to be published. This is the only guarantee that data analysis is not manipulated toward artificially inflated significant results. We appreciate that clinical trials are inherently different from microbiome (and other omics) studies, which are often exploratory in nature, but as far as statistics are concerned, they are the prey of the same traps and pitfalls. It is obvious that such a change in the approach to microbiome studies requires the concerted cooperation of researchers, journal editors, reviewers, and publishers.

85HMP data used in this study are available at https://github.com/mibwurrepo/Microbial-bioinformatics-introductory-course-Material-2018/tree/master/input_data; http://doi.org/10.5281/zenodo.1436630. The poultry data are available at https://www.ncbi.nlm.nih.gov/bioproject/ with accession number PRJNA553870. Simulated datasets are available at https://github.com/mibwurrepo/KersSaccenti-Power.

86JK and ES wrote the manuscript. Both authors read and approved the final manuscript.

87The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

88All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.