Hierarchical correction of p-values via an ultrametric tree running Ornstein-Uhlenbeck process

Statistical testing is classically used as an exploratory tool to search for association between a phenotype and many possible explanatory variables. This approach often leads to multiple testing under dependence. We assume a hierarchical structure between tests via an Ornstein-Uhlenbeck process on a tree. The process correlation structure is used for smoothing the p-values. We design a penalized estimation of the mean of the Ornstein-Uhlenbeck process for p-value computation. The performances of the algorithm are assessed via simulations. Its ability to discover new associations is demonstrated on a metagenomic dataset. The corresponding R package is available from https://github.com/abichat/zazou.


Introduction
In many fields, statistical testing is classically used as an exploratory tool to look for the association between a variable of interest and many possible explanatory variables.For example, in transcriptomics, the link between a phenotype and the expression of tens of thousands of genes is tested (McLachlan et al., 2005), in Genome Wide Association Studies (GWAS) the association between millions of markers and a phenotype is tested (Bush and Moore, 2012), in functional Magnetic Resonance Imaging (fMRI), the goal is to identify voxels that are significantly activated in two different conditions (Cremers et al., 2017).This problem of multiple comparisons dates back to the work of Tukey (Tukey, 1953).It has since been the subject of abundant literature and aims at controlling a probability of error of some sort.Most of the literature focus on the control of the Familiy Wise Error Rate (FWER) (Bland and Altman, 1995), being the probability of at least one false discovery among detections, or of the False Discovery Rate (FDR) (Benjamini and Hochberg, 1995), defined as the expected proportion of false positives among detections.Most of the correction procedures for controlling FWER or FDR, such as the popular Benjamini-Hochberg (BH) procedure, rely on independence, or some form of weak dependence, among the hypotheses, which is rarely observed in practice.Multiple testing under dependence is a difficult problem occurring in many fields.In transcriptomics, differential analysis has to deal with gene expressions that are often highly correlated.When performing GWAS, the linkage desiquilibrium imposes a strong spatial dependence between markers, and in Functional Magnetic Resonance Imaging (fMRI), two spatially close voxels have often comparable activation.The control of the FDR remains valid under arbitrary dependency structures by replacing the BH procedure with the more conservative BY procedure of Benjamini and Yekutieli (2001).However, based on results obtained from simulated datasets, it is obvious that there is a substantial loss of power when the real dependency structure is not taken into account, as discussed in depth in Blanchard et al. (2020).An alternative approach for dealing with multiple testing is to reduce the number of tests by aggregating certain hypotheses.Aggregation strategies vary and can be based on a priori knowledge (e.g.metabolic pathways, functional modules of genes) or on clustering algorithms (Renaux et al., 2020;Sankaran and Holmes, 2014).This article aims to take into account the dependencies between variables in order to offer a powerful statistical procedure of multiple testing.A hierarchical dependency structure between variables is assumed to be known up to certain constants.This assumption is common in our motivating example of microbiome studies (Huang et al., 2021;Matsen IV and Evans, 2013;Sankaran and Holmes, 2014;Silverman et al., 2017;Xiao et al., 2017), where the phylogeny is a natural hierarchical structure encoding similarities between variables (or namely species in that context).The hypotheses tested can then be organized in a tree structure which captures correlations at different scales of observation.This type of hierarchical structure is observable in transcriptomics differential analysis, where gene expressions can easily be represented by a hierarchy based on gene expression correlation.In GWAS and fMRI, spatial dependence also proves to be very suitable for hierarchical modeling (Ambroise et al., 2019;Eickhoff et al., 2015;Sesia et al., 2020).We propose to model the hierarchical structure of the multiple tests through an Ornstein-Uhlenbeck process on a tree.The process correlation structure is used for smoothing the p-values, after conversion to z-scores, similarly to the algorithm proposed in Xiao et al. (2017) but with an explicit underlying model.We then consider a three stage approach for our differential analysis procedure.The first stage reframes the initial problem as a linear regression problem that preserves the hierarchical structure.This linear problem is ill defined (p ∼ 2n) and we therefore resort to an 1 penalized estimation of the mean of the Ornstein-Uhlenbeck process.The second stage produces asymptotically valid p-values.The output of 1 penalized estimation produces are indeed biased and offer no theoretical guarantees about their asymptotic distribution; we therefore correct them using a debiasing procedure (Javanmard andMontanari, 2013, 2014;Zhang and Zhang, 2014) to compute valid p-values.The third and final stage controls the FDR of the overall procedure, using the tuning strategy of Javanmard et al. (2019).The selection strength of the Ornstein-Uhlenbeck process and the penalty parameter are hyperparameters of our model, whose selection is achieved via a Bayesian Information Criterion (BIC).We provide some background on hierarchical procedures in Section 2, introduce the model and statistical procedure in Section 3 and detail the computational steps in Section 4. The performances of the algorithm are assessed via simulations in Section 5.The use of the proposed model is illustrated in Section 6, where we demonstrate its ability to discover novel associations in a metagenomic dataset.

Examples of multiple testing strategies
A classic example in genomics consists in grouping the markers according to whether they belong to the same genes (aggregation by a prior).The genes can then be grouped according to their similarity, computed for example from expression profiles.Kim et al. (2010) have, for example, proposed a hierarchical testing strategy controlling the FWER in a hierarchical manner, by testing clusters of genes, then individual genes associated with a phenotype with the goal of finding genomic regions associated with a specific type of cancer.This type of top-down approach uses the concept of sequential rejection principle (Goeman and Finos, 2012;Meinshausen, 2008;Renaux et al., 2020).fMRI is another domain where tests are aggregated: neighboring voxels that are highly correlated are aggregated into a single voxel cluster.Benjamini and Heller (2007) propose an adaptation of the False Discovery Rate (FDR) to allow for cluster-level multiple testing for fMRI data.Ad hoc aggregating methods for multiple testing also exist in Metagenomics.LEfSe (Segata et al., 2011) performs a bottom up approach where a factorial Kruskal-Wallis rank sum test is applied to each feature with respect to a class factor, followed by a pairwise Wilcoxon test, and a linear discriminant analysis.MiLineage (Tang et al., 2017) performs multivariate tests concerning multiple taxa in a lineage to test the association of lineages to a phenotypic outcome.

Independence assumption
The assumption of independence of tests is convenient as it enables for both exact analyses and simple error bounds for classical procedures (Benjamini and Hochberg, 1995, e.g.).It is however unrealistic in practice.In many fields, including all the previous examples, measurements typically exhibit strong correlations.Some correction procedures, like the one proposed by Benjamini and Yekutieli (2001), make few assumptions while guaranteeing control of the FDR.Those general guarantees come with a high cost in terms of statistical power: the nominal FDR typically is much smaller that the target, resulting in many FN.Permutation procedures are an appealing alernative that can automatically adapt to the dependence structure of the p-values (Tusher et al., 2001) but may fail when confronted to unbalanced design or correlated data.Knowledge of the correlation structure can be leveraged to increase the power while still controling the FDR below a given target.Several approaches have been developed along those lines when the tests are organized along a hierarchical structure, typically encoded in a tree.

Hierarchical testing
The Hierarchical FDR (hFDR) introduced by Yekutieli (2008) and implemented in the R package structSSI (Sankaran and Holmes, 2014), proposes a topdown algorithm to sequentially reject hypotheses organized in a tree.The same approach is used in (Renaux et al., 2020) to select a group of variables arranged in a clustering tree.However, this approach suffers from some limitations, as shown in (Bichat et al., 2020;Huang et al., 2021).First, the algorithm in its vanilla formulation commonly fails to move down on the tree because of failure to reject the topmost node.Second, it only controls for an a posteriori FDR level, which is a complex function of the (user-defined) a priori FDR level and the structure of rejected nodes.This makes it difficult to calibrate the a priori FDR that would achieve a target a posteriori FDR and thus to compare it to other correction methods.Finally, it does not produce a corrected p-value, or q-value, per leaf, but only a reject / no reject decision and was shown in (Bichat et al., 2020) to perform no better than BH in many instances.Given all these drawbacks, we did not include the hFDR in our benchmark and use BH as a baseline instead.
StructFDR (Xiao et al., 2017) was developed for metagenomics Differential Abundance Testing (DAT) and relies on z-scores / p-values smoothing followed by permutation correction.Given any taxa-wise DAT procedure, p-values p are first computed for all m taxa (i.e.leaves of the tree) and then transformed to z-scores z.The tree is used to compute a distance matrix (D i,j ) and then turned into a correlation matrix C ρ = (exp (−2ρD i,j )) between taxa using a Gaussian kernel.The z-scores are then smoothed using the following hierarchical model: where µ captures the effect size of each taxa and z is a noisy observation of µ.The maximum a posteriori estimator µ * of µ is given by The FDR is controlled by means of a resampling procedure to estimate the distribution of µ * under H 0 and estimate adjusted p-values q sf .This method is implemented in the StructFDR package (Chen, 2018).
TreeclimbR (Huang et al., 2021) is a bottom-up approach also developed for metagenomics DAT but with a broader scope.It relies on aggregating abundances at each node of the tree (understood as a cluster of taxa) and performing a test to compute one p-value per node (compared one test per leaf for StructFDR).The main idea is then to use those p-values to compute a score for node i where B(i) is the set of descendants of node i, p k and s k ∈ {−1, +1} are the p-value of the node k and the sign of the associated effect, and t is a tuning parameter.A node i will be considered as candidate if U i (t) 1 and p i < α.This ensure that all descendants are (i) significant at level t with (ii) effects of coherent sign.At the end, multiplicity correction is only done on nodes (including leaves) that do not descend from another candidate.

Models and algorithms
Our correction methods assumes that p-values, or rather z-scores, evolve according to an Ornstein-Uhlenbeck process on a tree.We thus use the corresponding correlation structure to decorrelate the z-scores and, in turn, the p-values.This is similar in spirit to the smoothing algorithm of Xiao et al. (2017) but we derive our procedure from first principles and explicit assumptions.We first remind a few properties of Ornstein-Uhlenbeck processes before proceeding to our model and procedure.

Ornstein-Uhlenbeck process on a tree
An Ornstein-Uhlenbeck (OU) process (W t ) with optimal value (also called drift) β ou , selection strengh (also called mean reversion parameter) α ou and variance of the white noise σ 2 ou , is a Gaussian process that satisfies the stochastic differential equation: The important properties of OU processes are bounded variance and convergence to a stationary distribution centered on the optimal value β ou , namely ou /2α ou when t → ∞.Thanks to those properties, OU processes have become a popular model applied in various subfields of biology, ranging from evolution of continuous traits, such as body mass (Freckleton et al., 2003), fitness (Lande, 1976) or CpG enrichment in viral sequences (MacLean et al., 2021) to animal movement (Dunn and Gipson, 1977) and epidemiology (Nåsell, 1999).They naturally emerge as the continuous limit of broad range of discrete-time evolution models (Lande, 1976).Ornstein-Uhlenbeck processes can be readily adapted to tree-like structures as illustrated in Fig. 1.Formally, we consider a rooted ultrametric tree T with m leaves and n branches (n = 2m − 1 for binary trees).The internal nodes are labeled N 1 (the root) to N n−m and the leaves T 1 to T m .Let i be a node, W i the value of the trait at that node and denote pa(i) its unique parent.By convention, we set t N1 = 0 and assume W N1 = 0.The branch leading to i from pa(i) is denoted b i and has length l i = t i − t pa(i) where t i is the time elapsed between the root and node i.Since the tree is ultrametric, t i = h for all i ∈ {T 1 , . . ., T m }.For any pair of nodes (i, j), let t ij be the time elapsed between the root and the most recent common ancestor of i and j and denote d ij = t i − t j − 2t ij the distance in the tree between nodes i and j.The distribution of the trait at node i is given by: where λ i = exp(−α ou l i ) and β ou,i is the optimal value on branch i. Remark that the process mean value does not immediately shift to β ou,i but lags behind it with a shrinkage parameter controlled by 1 − λ i .If β ou,i = 0 for all i, straightforward computations show that W = (W T1 , . . ., W Tm ) is a gaussian vector with distribution When, the optimal value can shift on a branch (e.g. the branch b N4 leading to N 4 in Fig. 1), the mean vector of W is a slightly more complex and depends on both the tree topology and the location and magnitude of the shifts.Denote U the m × (n + m) incidence matrix of T with rows labeled by leaves (i ∈ {T 1 , . . ., T m }) and columns labeled by inner nodes and leaves (j ∈ {N 1 , . . ., N n−m , T 1 , . . ., T m }), with entries defined as U ij = 1 if and only if leaf i is in the subtree rooted at node j.Intuitively, column U .j encodes all leaves descending from node j and row U i. encodes all ancestors of leaf i. Denote ∆ the dimension n column vector with entries defined as ∆ i = β ou,i − β ou,pa(i) where i ∈ {N 1 , . . ., N n−m , T 1 , . . ., T m }.Non-zero entries of ∆ correspond to shifts location, nodes for which the optimal value β ou,i differ from its parent's and their values to shifts magnitude (see Figure 2 for an example).Finally let Λ be the n diagonal matrix with diagonal entries Bastide et al. (2017) for detailed derivations) show that W is a gaussian vector with joint distribution: (2) When T is known, the matrix T = U Λ is completely specified up to parameter α ou .The shifted Ornstein-Uhlenbeck model, with parameters α ou , σ 2 ou and shift vector ∆, has been used (Bastide et al., 2017;Khabbazian et al., 2016) to find Figure 1: (A) Phylogenetic tree with 5 leaves and 4 internal nodes (root N 1 included).A shift occurs on the branch leading to N 4 .(B) Ornstein-Uhlenbeck process with shifts on the tree defined in the left panel.At each node, the process spawns two independent process with the same initial value.The shifts on the optimal value on the branch leading to N 4 results in a different mean value for N 4 and all its offsprings (T 1 and T 2 ).adaptive events, modeled as non zero values in ∆, in the evolution of continuous traits of interest (turtle shell size, great monkey brain shape, etc).In this work, we apply the same mathematical framework to the joint distribution of p-values transformed to z-scores.

Procedure
We show here how to use the previously described Ornstein-Uhlenbeck process to incorporate the tree structure T in the correction of the p-values vector p.
Framework.Noting m 1 i (resp.m 2 i ) the median count (or relative abundance) of taxon i under condition 1 (resp.condition 2), we want to test i and assume that we have a testing procedure that outputs p-values, e.g. the Wilcoxon-Mann-Whitney test (Mann and Whitney, 1947;Wilcoxon, 1992).We first convert the p-values to z-scores using the quantile function Φ −1 of the standard gaussian: Provided the use of a correct statistical test, we known that p i ∼ U([0, 1]) under H i0 , so that z i ∼ N (0, 1).We also know that p i U([0, 1]) and thus z i N (0, 1) under H i1 .We could also test we only require the procedure to output p-values that satisfy the previous distributional assumptions for these H i0 and H i1 .Note that, even if the test statistic is itself a z-score before being transformed to a p-value, the z-score z i may differ from the raw test statistic z i because of the intermediate Figure 2: Incidence matrix U , shift vector ∆ and mean vector µ associated with Fig. 1.Λ N4 = 1−e αou(h−t N 3 ) is the shrinkage parameter from equation (1).
p-value p i .Indeed when considering the simple case of testing equality of means in two samples of size n, with gaussian distributions and known variance σ, the relation between z i and z i = √ n( m1 i − m2 i )/2σ is given by: After transformation, the test can be thus always be reframed as one-sided on z i : We make two assumptions regarding the distribution of z.
Assumption (A1) is very classic when working with z-scores (McLachlan and Peel, 2000): finding the alternative hypotheses is equivalent to finding the negative entries of µ.Assumption (A2) allows us to specify the joint distribution of z as: where Σ is fully specified by the parameters σ ou and α ou .Note that the diagonal coefficients of Σ are all equal to σ 2 ou /2α ou (1 − 2e −2αouh ).As they correspond to marginal variances, this forces the equality σ 2 ou = (1 − 2e −2αouh )/2α ou so that Σ depends only on α ou , i.e.Σ = Σ(α ou ).Finally, the decompositon µ = T ∆, where T acts as a phylogenetic design matrix, ensures that alternative hypotheses are likely to form clades, i.e. groups of leaves obtained by cutting a single branch in the tree.
This framework allows us to use T as a prior structure in the mean vector µ and variance matrix Σ and to recast the hypothesis testing problem as a regression problem.

Parameter Estimation
Estimation of μ.Assume first that Σ, or equivalently α ou , is known.Our main goal is to estimate the negative components of µ.
To leverage the known tree structure, we use the decomposition µ = T ∆ and estimate µ by means of ∆.Since ∆ has dimension n compared to dimension m for µ, we force ∆ to be sparse using a constrained lasso penalty (Tibshirani, 1996) : ∆ = argmin where Intuitively, the decomposition together with the 1 penalty works as a nested group lasso penalty for the components of µ, where the groups correspond to clades of T , while the constraint T ∆ ∈ R m − forces components of µ to be non positive.For compacity, we define the feasible set Finally, we use the Cholesky decomposition Σ −1 = R T R to simplify the problem into the very well studied optimisation problem: with y = Rz ∈ R m and X = RT ∈ R m×n .Note that y is a whitened version of z, with independent components and spherical covariance matrix.This is a lasso problem with a convex feasability constraint on ∆.The optimisation algorithm used to solve this problem is detailed in Section 4.
Estimation of Σ and tuning of λ.Remember first that Σ is completely determined by α ou because of the link between α ou and σ 2 ou .There are no closedform expression for the maximum likelihood estimator of α ou .We therefore resort to numerical optimisation.To tune the parameter λ, we test several values to estimate models with different sparsity levels and select the best one using a modified BIC criterion: (6) where ∆ α,λ is the solution of problem (4) for Σ(α) and λ.In practice, α and λ vary in a bidimensional grid and we select the values that minimize the objective.We use a modified BIC, where log(log m) log m replaces log m, to account for the fact that m scales like n as suggested in Fan and Tang (2013).

Confidence intervals
Lasso procedures are known to produce biased estimators and do not return confidence intervals for the point estimate μi .Instead of simply returning all negative components of μ = T ∆, we first debias the estimates and construct confidence intervals for the components of ∆, and in turn of μ, using the debiasing procedure of Javanmard andMontanari (2013, 2014); Zhang and Zhang (2014).
Debiasing.All debiasing procedures assume a model Y ∼ N m X∆, σ 2 I m and require both an initial estimator ∆(init) of ∆ and σ of σ.We use the scaled lasso (Sun and Zhang, 2012) with the same negativity constraint as in (4): Problem ( 7) can be solved efficiently by iterating between updates of (i) σ using the closed-form expression σ = y − X ∆ 2 / √ m and (ii) of ∆ by solving the constrained lasso problem (5) with tuning parameter λ scaled = λmσ.Debiasing is achieved by the corrected update: where the s j form a score-system (SS).Intuitively, s j should form a relaxed orthogonalization of x j against other column-vectors of X.The s j are used to decorrelate the estimators.We used the strategy of Zhang and Zhang (2014) and take the residuals of a lasso regression of x j against X −j .We also considered the alternative debiasing strategy of Javanmard andMontanari (2013, 2014), which is based on a pseudo-inverse of Σ = X T X m .Their debiased estimate is again a simple update of the initial scaled lasso estimator: but the decorrelation matrix S is computed in a so-called colwise inverse approach (CI), by inverting Σ in a columnwise fashion.Column s j is solution of the optimization problem : where e j is the j th canonical vector and γ ≥ 0 is a slack hyperparameter.If γ is too small, the problem is not feasible (unless Σ is non singular).If γ is too large, the unique solution is s j = 0.
Confidence Interval.Zhang and Zhang (2014) showed that asymptotically ∆ ∼ N (∆, V ) with the covariance matrix V defined by Similarly, the columnwise-inverse estimator of Javanmard and Montanari (2013) has asymptotic distribution N (∆, V ) with variance matrix V = S ΣS T /m.For both procedures, the bilateral confidence interval at level α for ∆j is Note that the estimator of the i th component of µ can be written μi = t T i.

∆ with t T
i. the i th row of T .Its unilateral confidence intervals at level α is thus given by −∞, μi + t T i. V t i. φ −1 (1 − α) .We can thus simply check whether 0 falls in the interval to test H i0 : {µ i = 0} versus H i1 : {µ i < 0} at level α or compute the p-value of the one-sided test as:

FDR control
The debiasing procedure achieves marginally consistent interval estimation of the shifts ∆ but additional care is required to control the FDR when testing all components of µ simultaneously.We use the procedure proposed in Javanmard et al. (2019), which is specific to debiased lasso estimators, and relies on the 1/2 .Briefly, for FDR control at a given level α, let t max = √ 2 log m − 2 log log m and set: where R(t) = m i=1 1 {ti≤−t} is the total number of rejections at threshold t, or t = √ 2 log m if the previous expression is empty.Applying the procedure from Javanmard et al. ( 2019) strictly would replace 2m with m in the numerator, as we're considering one-sided tests instead of two-sided ones for µ i .However, numerical analysis showed that the extra 2 led to better control of the FDR and we thus kept it.Hypothesis H i0 is rejected if t i ≤ −t or in term of q-values if Since t itself depends on α, the corrected p-values depend on α, unlike in the standard BH procedure, where they only depend on the order statistics.

Algorithm
The algorithm 1 summarises our procedure.We call it zazou for "z-scores az Ornstein-Uhlenbeck".
Algorithm 1 Zazou procedure 1: Compute the vector p of raw p-values 2: Transform it to the vector z of raw z-scores 3: for values of α and λ varying in a grid do 4: Compute Σ, R, y and X

Sign-constrained lasso
Our inference procedure is based on very standard estimates but requires to solve the following constrained lasso problem: For arbitrary vector y and matrices X and T .This a convex problem as both the objective function and feasibility set are convex.We therefore adapt the shooting algorithm (Fu, 1998), an iterative algorithm used to solve the standard lasso by looping over coordinates and solving simpler unidimensional problem, to our constrained problem.Let X −j (resp.∆ −j ) be the matrix X (resp.vector ∆) deprived of its j th column (resp.j th coordinate).We can isolate ∆ j in (5) and decompose the objective as We can likewise decompose T ∆ = u j + v j ∆ j where u j = T −j ∆ −j ∈ R m and v j = t j .When updating ∆ j , we can thus consider the simpler univariate problem in θ: Let I + = {i : v i > 0} and I − = {i : v i < 0} and denote θ max = min I+ {−u i /v i } and θ min = max I− {−u i /v i } with the usual conventions that max(∅) = −∞ and min(∅) = +∞.Problem ( 13) is feasible only if (i) θ min ≤ θ max and (ii) for all i, v i = 0 ⇒ u i ≤ 0, in which case the feasible region is [θ min , θ max ].Computing the subgradient ∂h(θ) of h and looking for values θ such that 0 ∈ ∂h(θ) leads to the usual shrinked estimates: By convexity of h, the solution of ( 13) can be found by projecting the previous unconstrained minimum to the feasibility set.If problem ( 13) is feasible, its solution is thus given by where P I : u → max(θ min , min(u, θ max )) is the projection of u on the segment I = [θ min , θ max ].
5 Synthetic Data

Metagenomics
Metagenomics data are made up of three components.The first component is the count or abundance matrix X = (x ij ), with 1 ≤ i ≤ m and 1 ≤ j ≤ p, which represents the quantity of taxa i in sample j.The second component is a set of sample covariates, such as disease status, environmental conditions, group, etc.The final component is a phylogenetic tree which captures the shared evolutionary history of all taxa.When performing DAT, we are interested in taxa whose abundance is significantly associated to a covariate.Most DAT procedures proceed with univariate tests (one test per species) followed by a correction procedure.In the synthetic datasets, we consider discrete covariates only.Dozens of full-fledged testing pipelines are published each year, including some designed with omics data in mind.Since our goal is this study is to compare correction procedures rather than full testing procedures, we use Wilcoxon or Kruskall-Wallis tests, which are classical and widespread non parametric tests in metagenomics.

Simulations
Simulation scheme.We use the following simulation scheme: 1. start with a homogeneous dataset, 2. assign each sample to group A or B at random 3. select differentially abundant taxa in a phylogenetically consistent manner (diffentially abundant taxa) 4. apply a fold-change to the observed abundance of diffentially abundant taxa in group B.
This non-parametric simulation scheme was previously used in Bichat et al. (2020).We considered two variants for step 3, respectively called positive and negative.In the negative variant, differentially abundant taxa were selected randomly across the tree, so that the phylogeny is not informative.In the positive variant, taxa are instead selected in a phylogenetically consistent manner.Formally, the phylogeny was first used to compute the cophenetic (Sneath et al., 1973) distance matrix between taxa.A partioning around medoids algorithm was then used to create cluster of related species.One or more clusters were then picked at random and all species in those clusters were selected as differentially abundant.
For each fold-change (fc ∈ {3, 5, 10}), 500 simulated datasets were created, with a proportion of differentially abundant species ranging from 3 % to 35 %.For each simulation, we corrected p-values using no correction (Raw), BH procedure (BH), BY procedure (BY), StructFDR (TF) or our procedure with either score system (SS) or colwise inverse debiasing (CI), targeting in all instances a 5% FDR level.We compared the 6 procedures in terms of True Positive Rate (TPR), nominal FDR and AUC (Area Under the Curve).
Positive simulations.The results of positive simulations (i.e.where the phylogeny is informative) are shown in Figure 3.All correction methods have controlled the FDR at the target rate or below when the fold change is larger than 5.For smaller fold changes, both SS and CI variations of zazou exhibit nominal FDR slightly above the target level (up to 9% in the worst case).In all settings, BY had the lowest TPR, whereas TF was comparable to vanilla BH, in line with results of Bichat et al. (2020).Finally, zazou (both SS and CI variations) had the best overall TPR, with largest gains observed in the lowest fold-change setting.The higher than intended FDR of zazou methods suggests that the problem of finding an adequate threshold for p ss i is not completely solved by Javanmard et al. ( 2019) procedure.To assess the performance of zazou in a thresholdindependent manner, we also compared the AUC of all procedures.Fig. 4 shows that zazou (both variants) has higher AUC than all other methods.As reported previously, TF and BH are at the same level and BY has the lowest ROC curve.Focus on the beginning of left hand side side of the curve shows that zazou is more efficient starting from the first discoveries.Negative simulations.The negative simulations are designed to assess the robustness of our algorithm with respect to uninformative phylogenies, or equivalently mispecified hierarchies.Fig. 5 shows that, as expected, standard BH outperforms competing methods (in terms of AUC) when the tree is mispecified.Forcing an inadequate tree structure results in AUC losses ranging from 15 to 20 percentage points compared to no structure.The puzzling lack of AUC loss for the TF procedure is explained by an implementation trick: TF always performs BH correction in parallel to its hierarchical procedure and falls back to BH when the hierarchical procedure detects much fewer species than BH (Bichat et al., 2020;Xiao et al., 2017).

Application
We use our zazou procedure on a gut microbiota dataset from the Fiji Islands (Brito et al., 2016;Pasolli et al., 2017) to identify species that are differentially abundant between adults and children.The data sets consists in the abundances of p = 387 species among n = 146 islanders, split into 112 adults and 34 children.
To mimick the simulation study, we used Wilcoxon tests for the univariate tests.Without correction, 21 species were detected as differentially abundant at the 5% level.None of them remained significant after correction by BH, BY, TreeFDR or treeclimbR.By contrast, zazou detected differentially abundant species with both desparsification methods: 17 for SS and 6 for CI.Fig. 6 shows that they are not a strict subset of the 21 detected with no correction.Smoothing salvages some species that are closely related to one of the 21 without being significant on their own (red box in the figure).It also illustrate some numerical problems associated with colwise-inverse debiasing, which is highly sensitive to the choice of the slack hyperparameter γ.The window of relevant values for γ is narrow and too large or too small values γ respectively lead to no correction or a faulty p-value correction.

Conclusion
In this work, we introduced zazou, a new method for correcting p values in a hierarchical context.zazou is based on recasting the testing problem as a regression problem, under the framework of stochastic processes on an ultrametric tree, and using the tree topology as a regularization parameter.It outperforms competing methods, hierarchical (TreeFDR, TreeclimbR) or not (BH, BY) in terms of AUC but this does not translate immediately to superior results in terms of FDR and TPR.The threshold for rejecting hypotheses is turned out to be quite difficult to calibrate while controling the FDR and warrants further work.There are several other parts of the procedure that are not as powerful as expected.First, the BIC step used to select λ and in turn the number of shifts tends to choose models with very few shifts, and sometimes even none.In such instances, the relevance of the debiasing step is limited.Second, the correction procedure proposed by Javanmard et al. (2019) is too conservative for our purpose.It was indeed developed to control both the FDR and the directional FDR (i.e.proportion of Type S errors, where the effect size have the wrong sign, in the discoveries) whereas we only need to control the former.For both these steps, specific developments taking into account the sign constraint on μ and the structure of the topology matrix of tree T could lead to better performances for zazou.Antoine Bichat, Jonathan Plassais, Christophe Ambroise, and Mahendra Mariadassou.Incorporating phylogenetic information in microbiome differen-

Figure 3 :
Figure 3: Boxplots and average (red point) TPR and FDR across positive simulation settings.Each facet corresponds to a different fold-change (fc) and each boxplot is computed over 500 simulation replicates.All corrections control the FDR at the target level or slightly above but zazou (SS and CI) achieve higher TPR, especially for small fold changes.

Figure 4 :
Figure 4: AUC boxplots (top) and average ROC curves (bottom) across positive simulations settings.Facets correspond to fold-changes (fc).ROC curves are computed for each simulation and linearly interpolated over a fixed grid before being averaged.Each boxplot and each curve are computed over 500 replicates.In all settings, SS/CI have the highest AUC / ROC curve, followed by BH/TF while BY has the lowest values.

Figure 5 :
Figure 5: AUC boxplots (computed over 500 replicates) in negative simulations.BH outperforms SS and CI, highlighting the cost of imposing a mispecified hierarchical structure.

Figure 6 :
Figure 6: Phylogeny of the 387 species from the Fidji dataset with associated zscores (inner circle), evidence (middle circle) and detection status (outer circle) under different correction procedures.Species detected by zazou are generally close-by on the tree and often, but not always, detected by raw p-values.The red strip highlight the smoothing property of the procedure in a subtree where individual species are not detected when using independant univariate tests but are detected when accounting for the hierarchical structure.