alpha value structSSI hierarchical FDR

[csmiguel]

I am applying your workflow to analyze my metagenomics data. Thank you for publishing this versatile and powerful guide. I want to determine potential taxa associated to the rhizosphere of a plant. For this purpose I am following your recommendations on the Hierarchical multiple testing with structSSI. I also believe this is probably the best approach since it profits from the topology of the phylogenetic tree to carry out hypothesis testing across all depths in the tree, thus maximizing the information to be mined from the data (compared to other methods which that do not consider hierarchy such as implemented in DESeq2).

I would much appreciate your opinion when setting alpha value in the function hFDR.adjust. The algorithm stops whenever p-value from a node (the probability of finding differences in that node) falls above alpha, so shallower nodes/tips get no adjusted p-values (“NA”; the hypothesis is not rejected). I noticed you set in your example an alpha=.75. If I set this same value I only get 3 nodes evaluated because the p-values for my deepest nodes is very high. To overcome this and get adjusted p-values for all nodes/tips in the tree I set alpha to 1. To my understanding this will decrease the power of my analysis because the number of hypothesis evaluated increases. So:

1. Setting an alpha=1 is equivalent to doing non-hierarchical Benjamini-Hochberg procedure (such as implemented in DeSeq2)? would it be OK to interpret ajusted p-values using this alpha?
2. Is there trade-off between power (alpha set) and significant results? After checking the ASVs with most significant values their assigned taxonomy is very coherent with the environments of the metagenomic samples.

Thanks again

[krisrs1128]

Hi @sundamiguel. I think it's worth clarifying the role of alpha: it’s the False Discovery Rate (FDR) control rate within a single family of hypothesis. An FDR control rate of 1 is not useful -- all probabilities are less than 1. Admittedly, in our paper, the alpha = 0.75 is not a very meaningful bound, and it was chosen that high for the sake of illustration. As an aside, the hierarchical testing procedure comes with a guarantee on the tree-wide FDR of alpha * (#discoveries + #families tested) / (#discoveries + 1). [see table 1 at https://www.jstatsoft.org/article/view/v059i13]

With this context, hopefully the answers to your questions make more sense,

1. The adjusted p-values when you set alpha = 1 are not equivalent to doing a non-tree structured BH procedure. It’s possible to construct situations where the signal is so weak at a node that the procedure doesn’t descend down it’s subtree even when alpha = 1, it’s just that you wouldn’t expect this to happen in most situations.
2. There is a trade-off, but I think you might be confusing the direction of the association. Larger alpha give you more power, at the cost of a larger number of false positives. Intuitively, the lower your bar for calling a result significant, the higher power you'll have, at the cost of more potential false discoveries.

Besides trying to fiddle with alpha, you might want to consider two alternatives for dealing with the fact that the tree-wide p-values at higher internal nodes are not significant. Right now, the treePValues computes p-values by averaging data across descendant nodes and performing a t or F test. This is powerful in the case where there are many weak signals spread throughout the tree, but is weak in the case that there are a few strong signals (needles in the haystack). A better choice in this case is to use a Bonferroni adjusted p-value -- take the minimum p-value across descendants and multiple by the number of descendants. There are some more recent proposals for leveraging tree structure when testing trees. You might want to check out the approaches here [https://arxiv.org/abs/1705.07529] and here [https://arxiv.org/abs/1709.10250]