Je Remercie

Anne-Laure Boulesteix

Elisabeth Gassiat

Jean-Philippe Vert D'

Je Remercie Également

Philippe Berthet

Béatrice Laurent-Bonneau

Franck Picard

Patricia Raynaud-Bouret

Bertrand Thirion D'avoir

Je Souhaite

Emmanuel Barillot

Stéphane Boucheron D'avoir

Julien Chiquet

Cyril Dalmasso

Rigaill Guillem

Clément, Fanny, François Cathy

Laurent Jean-Marc Gersende

Nicole Pierre

Morgane Alia Dehman

Benjamin Sadacca

Nathanaël Guillermo Durand

Christophe Ambroise

Catherine Matias

Fabien Reyal

Etienne Roquain

Nathalie Vialaneix

Marie Chavent

moi et je suis honoré que des scientifiques de votre stature jugent et apprécient mes travaux

Keywords: Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction. Cambridge University scientific challenges .

A new class of "high throughput" biomedical devices (...) routinely produce hypothesis-testing data for thousands of cases at once. This is not at all the situation envisioned in the classical frequentist testing theory of Neyman, Pearson, and Fisher.

-Bradley Efron,

Chapter 1

Introduction and overview

This chapter provides a brief account of my scientific path, followed by context and motivation for the research described in this manuscript. After a description of the recent transition from hypothesis-driven to data-driven research, we introduce specific types of genomic data and associated biomedical questions, which will serve as a basis for the leading applications in this document. Next, we highlight some challenges and opportunities raised by the analysis of such of high-throughput genomic data. Finally, we give an overview of the contributions described in this document. 1.1 Foreword: scientific path My first professional experience was a one-year internship at Crédit Lyonnais in the Groupe de Recherche Opérationnelle (2000)(2001), where I enjoyed being involved in methodological research projects while being in contact with their application [R2]. However, this was also the exciting period where the human genome was first sequenced and I decided to dive into the domain of statistics applied to genomics. For my masters thesis project in 2003 I had the opportunity to work with bacterial sequences at the Statistique et Génome lab headed by Bernard Prum in Évry [R1]. After graduating from ENSAE (2003) I obtained a 18 month position as a research engineer in biostatistics in the Bioinformatics group at Institut Curie, which had just been created by Emmanuel Barillot. This experience convinced me of the richness of the statistical questions raised by genomic data analysis. I decided to start a PhD in Statistics (2004Statistics ( -2008; co-supervised by Stéphane Boucheron at Université Paris 7 and Emmanuel Barillot at Institut Curie) in order to find an appropriate balance between theory and application. I was lucky enough to maintain this balance by being hired as a post doc in Terry Speed's group at the Statistics Department of UC Berkeley (2008Berkeley ( -2010)).

When I got back to France I had the opportunity to join the Statistique et Génome lab in Évry for a one-year postdoc. In 2011, I obtained a tenured researcher position at CNRS, and stayed in Évry for five more years. I moved to Toulouse in 2016 where I joined the Institut de Mathématiques de Toulouse.

From hypothesis-driven to data-driven research

The number and size of available data sets of different types is both a consequence and a cause of scientific and technological breakthroughs. For example, the completion of the human genome sequencing in 2003 by an international consortium triggered the development of high-throughput molecular profiling technologies: microarrays, followed by massively parallel sequencing of group of cells ("bulk" sequencing), and now of individual cells ("singlecell" sequencing); see [START_REF] Goodwin | Coming of age: ten years of next-generation sequencing technologies[END_REF] for a review of sequencing technologies. An emblematic example of large-scale initiative to produce molecular data is the Cancer Genome Atlas (TCGA) project from the US National Cancer Institute. Quoting the web page of TCGA1 , The Cancer Genome Atlas (TCGA), a landmark cancer genomics program, molecularly characterized over 20,000 primary cancer and matched normal samples spanning 33 cancer types. (...) TCGA generated over 2.5 petabytes of genomic, epigenomic, transcriptomic, and proteomic data. The data, which has already lead to improvements in our ability to diagnose, treat, and prevent cancer, will remain publicly available for anyone in the research community to use.

This "data deluge" has been accompanied by a shift from hypothesis-driven research to data-driven research. The classical statistical approach to data analysis starts by defining a scientific hypothesis, collecting data, and performing inference. In contrast, the current practice in genomics starts by data collection and tends to become technology-or datadriven. Hypotheses are then defined based on this data, and inference is performed on these hypotheses 2 . An important consequence of this change of paradigm is the need for dedicated statistical methods. This point is illustrated in the next sections for the particular case of genomic data.

Instances of biomedical questions

In this section we introduce three biomedical contexts and questions, together with associated genomic data. We describe some of the statistical questions raised by the analysis of these data. The leading applications considered in this document will be related to these questions.

Differential expression studies

Differential gene expression studies in cancerology aim at identifying genes whose mean activity differs significantly between two (or more) cancer populations. The activity of a gene is obtained by gene expression measurements from individuals of these populations. These measurements are typically assessed by microarray or sequencing experiments, in which millions of features are collected and summarized in thousands of gene expression measurements. In this document, we specifically consider a Leukemia data set studied in [START_REF] Bourgon | Independent filtering increases detection power for high-throughput experiments[END_REF]. It consists of expression measurements for 12, 625 genes for biological samples from 79 individuals with B-cell acute lymphoblastic leukemia (ALL) [START_REF] Chiaretti | Gene expression profiles of B-lineage adult acute lymphocytic leukemia reveal genetic patterns that identify lineage derivation and distinct mechanisms of transformation[END_REF]. 37 of these individuals harbor a specific mutation called BCR/ABL. The goal of this study is to understand the molecular differences at the gene expression level between the mutated and non-mutated individuals in the population 3 . Perhaps the most basic question to ask is: for which genes is there a difference in the mean expression level of the mutated and non-mutated population? This question can be addressed, after relevant data preprocessing, by performing a statistical test of equality in means for each gene, and to derive a list of "differentially expressed" genes (DEG) as those passing some significance threshold. This is a typical example of multiple testing situation, which requires the definition of dedicated risk measures and associated methods to control these risks.

Other challenges raised by the analysis of gene expression data in cancer samples include the discovery of new cancer subtypes, or the prediction of clinical phenotypes such as survival or the severity of the disease. In statistical terms, these biomedical questions can be cast as problems of unsupervised or supervised classification, or prediction.

DNA copy number studies

Each normal human cell has 23 pairs of chromosomes. For each of them, one chromosome has been inherited from each biological parent. Tumor cells harbor numerous structural alterations of their DNA including point mutations, translocations, small insertion or deletion events, larger scale copy number changes, or amplifications. Figure 1.1 illustrates the effect of some of these events in the lung tumor cell line NCI-H1395-4W. The left panel shows copy number changes at the scale of the entire genome for one particular cell, while the right panel shows copy number changes at the level of one single chromosome for a large number of cells 4 . In these cells, chromosome 5 has four copies in the first 50 Mb and three copies in the rest of the chromosome.

Such copy number alterations can affect genes and regulatory transcripts, which may result in major cellular modifications and are associated with diagnostic and prognostic factors [START_REF] Hanahan | The hallmarks of cancer[END_REF][START_REF] Hanahan | Hallmarks of cancer: the next generation[END_REF]. An important issue in cancer research is to therefore to estimate the underlying copy number state (to be defined more formally in Chapter 9) at each position along the genome of a tumor sample. Microarray and sequencing-based technologies have been used in the last two decades to quantify copy numbers at a large number of genomic loci [START_REF] Chin | Translating insights from the cancer genome into clinical practice[END_REF][START_REF] Pinkel | High resolution analysis of DNA copy number variation using comparative genomic hybridization to microarrays[END_REF]. For instance, the copy number profile in Figure 1.1 was obtained from Affymetrix Genome-Wide Human SNP 6.0 Arrays, which contain more than 1.8 million genomic markers. In contrast to these "bulk" technologies which provide genomic information averaged over all cells in a sample, single-cell sequencing technologies have been developed in the last few years to provide genomic information at the level of individual cells [START_REF] Gawad | Single-cell genome sequencing: current state of the science[END_REF].

The analysis of DNA copy number data involves several steps. First, these data have to be preprocessed/normalized in order to make data points comparable across loci and possibly across samples. Then, they are segmented into regions of constant DNA copy number state, and these copy number states are labelled ("called") in a biologically-meaningful manner (e.g. normal, loss of one copy, etc). As for gene expression data, DNA copy number data also raise unsupervised and supervised classification questions, before or after segmentation.

Genome-Wide Association Studies

Genome-Wide Association Studies (GWAS) aim at identifying genomic markers associated with a phenotype of interest. This phenotype may be binary in the case of case-control studies, or continuous (e.g. the date of disease onset). In Chapter 10 we consider a GWA study with 615 patients infected by Human Immunodeficiency Viruses (HIV). One of the goal of this study was to detect genetic factors that influence the plasma viral load, that is, the level of HIV RNA in the patients blood. For each patient, more than 300, 000 genetic markers were assessed by an Illumina genotyping microarray experiment.

Like for differential expression analyses, state-of-the-art approaches start by testing the univariate association between each marker's genotype and the phenotype, and retaining the most significant markers using a multiple testing correction. This list is then interpreted at the genome scale using a Manhattan plot (as in Fig. 1.2, left panel, which is reproduced from [103, Figure 1]), where the (-log-transformed) p-values of the markers are plotted against their genomic position in order to identify which genome regions show an enrichment in significantly associated markers. The right panel in Fig. 1.2 (which is reproduced from [103, Figure 1]) shows a measure of statistical dependence (called linkage disequilibrium, LD) between pairs of genomic markers, together with relevant genomic annotation.

Other types of genomic information that can be assessed via high-throughput genomic experiments include:

• DNA methylation, a chemical transformation of one of the four letters of the DNA alphabet that can alter the expression of neighbor genes [START_REF] Chatterjee | Tools and strategies for analysis of genome-wide and gene-specific DNA methylation patterns[END_REF];

• the interaction between a protein of interest and genomic regions, which can be quantified by Chromatin Immuno-Precipitation (ChIP) [START_REF] Nakato | Recent advances in ChIP-seq analysis: from quality management to whole-genome annotation[END_REF];

• the intensity of physical interaction between two genomic regions, which can be quantified by Chromosome conformation capture techniques such as 3C, 4C, or Hi-C [START_REF] Dixon | Topological domains in mammalian genomes identified by analysis of chromatin interactions[END_REF]

Statistical challenges and opportunities

At first glance, the above-described biomedical questions can be cast as classical statistical problems of testing, prediction, unsupervised and supervised classification, or segmentation. However, genomic data have specificities that imply that classical statistical tools can generally not be readily applied to them, thus requiring new statistical developments. We refer to [42, Chapter 1] for a detailed typology of genomic data. In this section, we focus on specific features of these data -namely, their high-dimensionality, sparsity, heterogeneity and structuration -and try to analyze some of the challenges and opportunities raised by each of these features for statisticians.

High-dimensionality and sparsity

The number p of variables (genes, loci) is typically of the order of 10 3 to 10 6 , that is, several orders of magnitude larger than the number n of observations (biological samples, patients), which generally ranges from 1 to 10 3 . It is also generally the case that these data are sparse, in the sense that only a small (unknown) number of (unknown) variables or of combinations of variables actually contain signal. This high-dimensional setting is in sharp contrast with the classical statistical setting where we generally have n p and always n > p. As a consequence, even the most basic statistical tools like linear regression, or statistical tests have to be revisited to cope with this context. High-dimensionality is not specific to the field of genomic data, as it is the consequence of the general sophistication of data acquisition technologies not only in biomedical sciences (genomic and imaging data), but also in other fields such as physics, environmental sciences, economy, finance, to name but a few. A number of mathematical, statistical and computational tools dedicated to highdimensional data have been developed in the past 20-25 years. Two emblematic examples of these developments are the Least Absolute Shrinkage and Selection Operator (lasso) for penalized regression and variable selection [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] and the False Discovery Rate (FDR) for multiple testing [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF] 5 .

A possibly distinctive feature of genomic data compared to other data types is that each new technological advance results in a increase of one to several orders of magnitude for p, but with n typically remaining constant or even getting temporary smaller because of the prohibitive cost of new technologies. Moreover, the increase in the number of measured variables is unfortunately accompanied by an increasingly lower signal to noise ratio for each measurement. This can be illustrated with the example of sequencing technologies (from bulk to single cell sequencing) getting at the same time wider (larger p), sparser and noisier, see e.g. [START_REF] Hansen | Sequencing technology does not eliminate biological variability[END_REF] whose provocative title is: "Sequencing technology does not eliminate biological variability".

Curse or blessing?

The high-dimensional nature of genomic data is considered as a blessing by biomedical scientists (following the possibly misleading argument that "more data" implies "more knowledge"), and as a curse by mathematicians as nicely explained and illustrated in [46, Chapter 1] from several standpoints. In particular, n points in p-dimensional space with p n are typically isolated. This makes the notion of "neighborhood", which is fundamental in classical statistical inference, essentially useless in this setting. Moreover, the number of candidate models to explain a given phenomenon increases exponentially with the number of parameters of this model, which is typically indexed by p.

Fortunately, the sparsity of genomic data may also be seen as a blessing, because it implies that the biological signal is concentrated in a much lower-dimensional space than the p-dimensional observation space. For example, a typical assumption in prediction or variable selection tasks is that only a small number of variables, or of groups of variables, are actually relevant. Such assumptions are the basis of the methods recently developed for high-dimensional (genomic) data analysis. They are also assumptions under which statistical guarantees for such methods and efficient algorithms can be obtained, see e.g. [START_REF] Hastie | The Elements of Statistical Learning[END_REF][START_REF] Giraud | Introduction to high-dimensional statistics[END_REF].

Heterogeneity and structure

In the classical statistical framework, inference is generally performed on a single n×p table. By heterogeneous, we mean the situation where several levels of biological information are available for the same set of observations. For example, a tumor biopsy can be analyzed to study mutations, DNA copy numbers, DNA methylation, gene or protein expression, or gene regulation. Moreover, these genomic data can also be complemented by other data types such as imaging data (e.g. phenotyping experiments, or histology) or clinical data coming from electronic health records (EHR). This setting is called multi-view data in the machine learning literature, see e.g. [START_REF] Xu | A survey on multi-view learning[END_REF] for a survey. Moreover, the measured variables are often structured, that is, linked by networks, by similarity relationships, or simply by a natural ordering along a chromosome. This may be illustrated by the above examples: in differential expression studies, sets of genes can be codifferentially expressed because they belong to the same gene network or pathway; in DNA copy number studies, neighboring loci on a chromosome are expected to have identical copy number state (as in Figure 1.1, right); in GWAS, a SNPs can have its genotype associated to a phenotype because its belongs to the same block of LD (see Figure 1.2) as a causal SNP.

Curse or blessing? The above-described characteristics of genomic data represent a challenge for the statistician, because standard statistical tools are likely not to be suited to address a new biomedical question, and because devising a tailored method requires some basic understanding of the question at hand. I believe that a relevant way to tackle this complexity is precisely to take advantage of the heterogeneity and structuration of genomic data by considering them as constraints that alleviate the curse of dimensionality, and guide statistical methods toward solutions that are more plausible from a biological standpoint.

Overview of contributions

In the above-described context, my research objective has been to contribute to the development, mathematical study, and practical application in interdisciplinary research projects of statistical tools that take advantage of the complexity and structure of genomic data. My contributions generally result of the quest for an acceptable trade-off between mathematical rigor, computational efficiency, and biological interpretability. It is often the case that the intersection between what is mathematically justified, algorithmically feasible, and biologically interpretable is empty for a given problem. In such situations, I strive to make one or more of these three sets bigger, or the gap between these objectives smaller. Depending on the scientific question and the state of the art for a given problem, my contributions can be of different types:

• proposing new statistical methods, associated algorithms and their implementation

• establishing statistical properties for new or existing methods

• evaluating the statistical and computational performance of new or existing methods Organization of the manuscript. The remainder of this manuscript is organized in two main parts.

Part I summarizes my contributions to the field of multiple testing, and post-selection or post hoc inference. After an introduction to multiple testing (Chapter 2), I describe my contributions to the asymptotic properties of FDR controlling procedures (Chapter 3), to FDR thresholding for classification under sparsity assumptions (Chapter 4) and to post hoc inference (Chapter 5)

Part II gathers my other contributions to statistical inference for genomic data. As explained in a short introductory chapter (Chapter 6), which can be categorized in two broad themes:

• Inference from heterogeneous genomic data, where inference is performed by combining several data types, either corresponding to the same observations (Chapter 7) or to the same variables (Chapter 8);

• Inference from ordered genomic data, where one of the main statistical challenges is to segment a genome (or more precisely each of its chromosomes) into successive homogeneous regions. This is tackled via segmentation methods for DNA copy numbers (Chapter 9) and by constrained clustering methods for GWAS and Hi-C studies (Chapter 10).

Finally, some directions for future research are discussed in Chapter 11.

Part I

From multiple testing to post hoc inference Chapter 2

Introduction to multiple testing

Significance testing was introduced in the early 1900s with the works of Fisher, Neyman and Pearson. The examples of differential expression analyses and GWAS taken in Chapter 1 require thousands to millions of statistical tests to be performed simultaneously. In such situations, it is desirable to control a global risk measure associated to the entire set of tested hypotheses. Large-scale multiple testing is concerned with the definition of such risk measures, the formalization of mathematical assumptions under which these risks are effectively controlled, and the construction of dedicated algorithms (called multiple testing procedures). This chapter provides an overview of large-scale multiple testing theory, with emphasis on concepts that will be useful to the reader for the next chapters of this part. This chapter draws from several reference books or surveys on multiple testing [START_REF] Roquain | Contributions to multiple testing theory for high-dimensional data[END_REF][START_REF] Dickhaus | Simultaneous statistical inference[END_REF][START_REF] Goeman | Multiple hypothesis testing in genomics[END_REF][START_REF] Roquain | Type I error rate control in multiple testing: a survey with proofs[END_REF][START_REF] Lehmann | Testing statistical hypotheses[END_REF].

Remark: Throughout this part, the number of hypotheses tested is denoted by m whereas in Part II the number of variables will be denoted by p.

Statistical setting

We consider an observation x ∈ X , that is, a realization of a random variable X valued in a measurable space (X , X). Our statistical model is a family P of candidate probability distributions for the true distribution P of X on (X , X). Given a subset P 0 of P, we consider the null hypothesis H 0 : P ∈ P 0 . The corresponding alternative hypothesis is H 1 = P ∈ P \ P 0 unless otherwise noted. Given an observation x ∈ X , we aim at deciding whether P is compatible with the null hypothesis, that is, whether P ∈ P 0 . In order to do so, following e.g. [START_REF] Casella | Statistical inference[END_REF], we define a p-value as a random variable p(X) on [0, 1] satisfying1 :

∀P ∈ P 0 , ∀u ∈ [0, 1], P(p(X) ≤ u) ≤ u . (2.1)
By construction, the p-value p(x) is a quantitative measure of how unlikely observation x is if H 0 is true. This measure can be used to take a decision as to whether H 0 is true or not. Given α ∈ [0, 1], a test of level α "rejects" H 0 if and only if p(x) ≤ α. In order to compare the decision of rejection/acceptance of the null hypothesis to the truth, we introduce the classical vocabulary of true/false positive/negative in

Multiple testing setting

For m ∈ N, we consider a collection of subsets P 0,i of P and the associated null hypotheses: H 0,i : P ∈ P 0,i indexed by N m = {1, . . . , m}. The corresponding alternative hypotheses are H 1,i = P ∈ P \ H 0,i . each i ∈ N m . From the collection H = (H 0,i ) i∈Nm , we define a m-dimensional parameter of interest: (θ i (P)) i∈Nm , where

θ i (P) = 1 {P ∈ H 0,i } for i ∈ N m .
The sets of true null hypotheses and true alternative hypotheses are denoted by H 0 (P) and H 1 (P), respectively. When non ambiguous, we omit P in the notation and write θ i , H 0 and H 1 for simplicity. We let m 0 (m) = i∈Nm (1

-θ i ) = |H 0 | and m 1 (m) = i∈Nm θ i = |H 1 |
be the number of true null hypotheses and of true alternative hypotheses, respectively. Finally, π 0 (m) = m 0 (m)/m is the proportion of true null hypotheses.

Location model. We introduce a standard location model which will be used several times for illustration:

X i = µ i + ε i , i ∈ N m , (2.2)
where the ε i are identically distributed, with a common marginal distribution that is assumed to be continuous, and the location parameter µ i is null if and only if θ i = 0. A location model is characterized by the distribution of the test statistics conditionally on θ i = 0. The most common instances of location models are the Gaussian and Laplace (double exponential) models. Both are particular instances of the Subbotin location model, where the density of the test statistics conditionally on θ i = 0 is given by the ζ-Subbotin density:

d(x) = (L ζ ) -1 e -|x| ζ /ζ , with L ζ = +∞ -∞ e -|x| ζ /ζ dx = 2Γ(1/ζ)ζ 1/ζ-1 , (2.3)
for ζ ≥ 1. The Gaussian and Laplace models correspond to ζ = 2 and ζ = 1, respectively.

In Section 4 we will also consider scaling models based on the Subbotin density.

Multiple testing procedures. We define a multiple testing procedure as a mapping from an observation X to a subset of N m , corresponding to a set R m of rejected hypotheses.

For notation simplicity we will identify a multiple testing procedure with the associated rejection set R m . We denote by |R m | the corresponding number of rejections, and often omit the subscript m when m is fixed. We assume to be given a set of p-values associated to H 0,i for each i ∈ N m , that is, that there exists a random variable p i (X) on [0, 1] satisfying:

∀P ∈ P 0,i , ∀u ∈ [0, 1], P(p i (X) ≤ u) ≤ u . (2.4)
We consider multiple testing procedures that reject the hypotheses whose p-value is less than a (possibly random and data-driven) threshold. Such p-value thresholding-based multiple testing procedures (or thresholding procedures for short) may be written as

R m = {i ∈ N m , p i ≤ tm } ,
where tm is called the threshold of the multiple testing procedure. More precisely, we consider three classes of procedures. In "single-step" procedures, the rejection of a given hypothesis does not depend on the rejection of other hypotheses. In contrast, "step-up" and "step-down" procedures are part of the broader class of stepwise multiple testing procedures, where the rejection of a given hypothesis may depend on the rejection of other hypotheses. • The step-up procedure with critical values c is the procedure with threshold t = c ι↑ (c) , where

ι↑ (c) = max i ∈ N m , p (i) ≤ c i . (2.5)
• The step-down procedure with critical values c is the procedure with threshold t = c ι↓ (c) , where ι↓ (c) = max i ∈ N m , ∀j ≤ i, p (j) ≤ c j .

(2.6)

The threshold of a step-down procedure is the first crossing point between the ordered p-values and c, while the threshold of a step-up procedure is the last such crossing point. By definition, for a fixed family of critical values c, the associated step-up procedure rejects at least as many hypotheses than the associated step-down procedure.

Multiple testing risks

The (unobserved) number of false positives of a multiple testing procedure is |R m ∩ H 0 |, which implicitly depends on the distribution P of X through H 0 . Historically, the first risk measure considered in a multiple testing context is the Family-Wise Error Rate (FWER). It is defined as the probability of (at least) one false rejection:

FWER P (R m ) = P(|R m ∩ H 0 | > 0) . A natural generalization is k-FWER P (R m ) = P(|R m ∩ H 0 | ≥ k) which allows at most k false rejections. A multiple testing procedure R m is said to control k-FWER (strongly) at level α ∈ [0, 1] if k-FWER P (R m )
≤ α for all P ∈ P. A much less demanding criterion called weak k-FWER control consists in ensuring that k-FWER P (R m ) ≤ α for all P ∈ ∩ i∈Nm P 0,i . In other words, a weak k-FWER controlling procedure simply corresponds to a test of level α of the global null hypothesis ∩ i∈Nm H 0,i . Here, we focus on strong control. Another quantity of interest is the False Discovery Proportion (FDP), which is defined as the fraction of true null hypotheses among those rejected:

FDP P (R m ) = |R m ∩ H 0 | |R m | ∨ 1 .
As the FDP is a random quantity, one possibility to define an associated risk is to focus on its expectation under P, which is known as the False Discovery Rate (FDR), or its quantiles under P (see the FDX below). The FDR is defined by

FDR P (R m ) = E |R m ∩ H 0 | |R m | ∨ 1 ,
and a multiple testing procedure R m is said to control FDR (strongly) at level α ∈ [0, 1] if FDR P (R m ) ≤ α for all P ∈ P. A related quantity is the positive FDR [START_REF] Storey | The Positive False Discovery Rate: A Bayesian Interpretation and the q-value[END_REF], defined by

pFDR P (R m ) = E |R m ∩ H 0 | |R m | |R m | > 0 .
Finally, we also define the False Discovery Exceedance [START_REF] Genovese | Exceedance Control of the False Discovery Proportion[END_REF] as the tail probability of the FDP for a given q ∈ [0, 1]:

FDX(R m ) = P |R m ∩ H 0 | |R m | ∨ 1 ≥ q .
Naturally, (strong) pFDR and FDX control can be defined similarly as for FDR.

Assumptions on the joint p-value distribution

We consider several types of assumptions on the joint distribution of the p-value family (or, equivalently, on the set P of possible distributions for P). We introduce the notation p A = (p i ) i∈A for A ⊂ N m . Two extreme distributional assumptions are general dependence, where the p-value distribution is let arbitrary, and independence: p H(P) is a family of mutually independent variables (indep)

Remark 2.2. As the statistical risks defined above only focus on type I errors (false positives), any result stated below under (indep) is also valid under the weaker assumption that p H0(P) is a family of mutually independent variables, and is also independent from p H1(P) . In practice, H 0 is unknown so we chose to use Assumption (indep) for simplicity.

Assumption (indep) is useful in theory because it simplifies the statistical study of multiple testing risks, but it is unrealistic in applications. A weaker assumption is Positive regression dependency on a subset of hypotheses (PRDS) 2 . The set S ⊂ [0, 1] m is nondecreasing if for all (q, q ) ∈ ([0, 1] m ) 2 such that ∀i ∈ N m , q i ≤ q i , we have: q ∈ S implies q ∈ S. We assume that the p-value family is PRDS on the subset H 0 (P) of true null hypotheses, that is:

for any i 0 ∈ H 0 (P) and any non-decreasing set S ⊂ [0, 1] m , the function u → P((p i ) i∈Nm ∈ S|p i0 ≤ u) is non-decreasing (PRDS(H 0 ))
Assumption (PRDS(H 0 )) is weaker than independence, and it is considered as realistic in genomics [START_REF] Goeman | Multiple hypothesis testing in genomics[END_REF], although it is generally not possible to check whether it holds for a particular application. A classical example of PRDS distribution of the p-value family is the following equi-correlated model:

Example 2.7. The test statistics are multivariate Gaussian, with a covariance matrix whose non-diagonal entries are all equal to ρ ∈ [0, 1].

More generally, Assumption (PRDS(H 0 )) holds whenever the test statistics are multivariate Gaussian, with a covariance matrix Σ whose entries are all non-negative. Example 2.7 will be used repeatedly for illustration in this part.

The main reason for the popularity of the PRDS assumption in multiple testing is that the Simes inequality [START_REF] Simes | An improved Bonferroni procedure for multiple tests of significance[END_REF] is valid under this assumption. We recall this inequality, together with Hommel's inequality [START_REF] Hommel | Tests of the overall hypothesis for arbitrary dependence structures[END_REF], a more conservative inequality valid under general dependence. We denote by p (1:I) ≤ • • • ≤ p (|I|:I) the ordered p-values of I.

Proposition 2.3 (Simes and Hommel's inequalities). We have:

P ∃i ∈ H 0 : p (i:H0) ≤ αi m ≤ α|H 0 | m c(m) , (2.8)
in the two following cases:

• c(m) = 1 and the p-value family satisfies (PRDS(H 0 )) (Simes' inequality);

• c(m) = C(m) := m i=1 1/i (Hommel's inequality). Moreover (2.8
) is an equality (with c(m) = 1) when the p i , i ∈ H 0 (P ), are i.i.d. U (0, 1).

In particular, the right-hand side of (2.8) is upper bounded by αc(m). The family of thresholds (αi/m) 1≤i≤m is called the Simes threshold family. The factor C(m) is of the order of log(m) and quantifies "the price to pay" for allowing for general dependence. The last statement in the above Proposition illustrates the sharpness of the Simes inequality, in the sense that they cannot be uniformly improved on the class of PRDS distribution. The same holds for Hommel's and for arbitrary distributions, respectively 3 . The Simes inequality will is also the fundamental ingredient of the original post hoc procedures proposed by [START_REF] Goeman | Multiple testing for exploratory research[END_REF], as will be discussed in Section 5.

Family-Wise Error Rate control

FWER control under general dependence can be obtained using a simple union bound argument, which is sometimes referred to as the Bonferroni inequality. For any I ⊂ N m and any t > 0,

P (∃i ∈ I, p i ≤ t) ≤ i∈I P (p i ≤ t) ≤ t|I| (2.9)
where the last inequality follows from the p-value property. Applying this inequality with I = H 0 and t = α/m ensures that

P (∃i ∈ H 0 , p i ≤ α/m) ≤ m 0 m α ≤ α (2.10)
The left-hand side is the FWER of the Bonferroni procedure with threshold t Bonf = α/m. Equation 2.10 demonstrates that the Bonferroni procedure controls FWER under general dependence, but also that it is conservative, in the sense that it controls FWER at level π 0 α for a target level α. Therefore, it is possible to obtain more powerful FWER controlling procedures by bridging the gap between π 0 α and α. The (Bonferroni-)Holm procedure [START_REF] Holm | A simple sequentially rejective multiple test procedure[END_REF] is a stepwise modification of the Bonferroni procedure which also controls FWER under general dependence, but possibly with an increased number of rejections. The Holm procedure is defined as the step-down procedure with critical values:

c i = α m -(i -1) , i ∈ N m . (2.11)
Another way to construct FWER controlling procedures that are less conservative than the Bonferroni procedure is to make dependency assumptions on the family of p-values. For completeness we mention that under independence, FWER is controlled by the Šidák procedure [START_REF] Šidák | Rectangular confidence regions for the means of multivariate normal distributions[END_REF], whose threshold is defined as

t Sidak = 1 -(1 -α) 1/m . Because 1 -(1 -α) 1/m ≥ α/m
, the Šidák procedure is less conservative than the Bonferroni procedure. However, when m is large, both thresholds are equivalent when α/m is small, so the Šidák procedure may only improve the number of rejections marginally, at the price of a much narrower applicability. Under (PRDS(H 0 )), the FWER is controlled by the step-up procedure with the same critical values as the Holm procedure (2.11), which is called the Hochberg procedure. The decreased conservativeness of the Hochberg procedure compared to the Holm procedure comes at the price of a narrower applicability.

Connection to closed testing

The Holm and Hochberg procedures can be derived from a general method for obtaining FWER control, known as closed testing [START_REF] Marcus | On closed testing procedures with special reference to ordered analysis of variance[END_REF]. This method can also be viewed as the basis of the post hoc procedures proposed by [START_REF] Goeman | Multiple testing for exploratory research[END_REF], that inspired our work [J2] described in Section 5. Closed testing focuses on the larger multiple testing problem of testing all possible intersections of null hypotheses. Formally, let us define the closure of the family of hypotheses H = (H 0,i ) i∈Nm as the set H = (H I ) I⊂Nm,I =∅ , where H I denotes the intersection hypothesis associated with I:

H I = i∈I H 0,i .
We denote for short by Nm = P(N m ) \ {∅} the set of all non empty subsets of N m . The set of true null intersection hypotheses is then H0 = I ∈ Nm , H I is true . Since H {i} = H 0,i for all i, we can write

P R ∩ H0 > 0 ≤ α
As H ⊂ H, R ∩ H can be seen as a multiple testing procedure for H. By Proposition 2.4, this procedure controls FWER (in H) at level α. To see why Proposition 2.4 holds, simply note that as H 0 ∈ Nm , φ H0 is an α-level test of H H0 . Therefore, there exists an event of probability larger than 1-α under which H H0 is not rejected by the closed testing procedure. Under this event, any true null intersection hypothesis, which is by definition a subset of H 0 , is also not rejected by the closed testing procedure.

The closed testing method thus provides an elegant and generic construction of FWERcontrolling procedures, where the properties of the procedure are inherited from the properties of the local tests. Closed testing can itself be seen as a consequence of the sequential rejection principle [START_REF] Goeman | The sequential rejection principle of familywise error control[END_REF], which provides generic sufficient conditions to build stepwise FWER controlling procedures. A caveat to the practical application of closed testing is that it implies testing all 2 m -1 possible non-empty intersections between m hypotheses. However, depending on the form of the local test, it may be possible to avoid testing all 2 m -1 hypotheses explicitly. This is the case in the procedures listed in the next paragraph.

The closed testing principle may be used to recover some of the procedures described above. First, the Holm procedure coincides with the closed testing procedure associated with local Bonferroni tests, which are defined by φ I = 1 {∀i ∈ I, p i ≤ α/|I|} for the null hypothesis H I . These local tests are of level α under general dependence by (2.9), so the closed testing principle ensures that the Holm procedure controls FWER under general dependence. Similarly, the Hochberg procedure can be seen as a closed testing procedure associated with local Simes tests, which are defined by φ I = 1 ∀i ∈ I, p (i) ≤ αi/|I| . These tests are of level α under (PRDS(H 0 )) by (2.8), so the closed testing principle ensures that the Hochberg procedure controls FWER under (PRDS(H 0 )) 5 . The procedures described in this section are summarized in Table 2.2. Although the Holm procedure is always at least as powerful as the Bonferroni procedure, the latter is much more widely used, for example in bio-medical applications.

Adaptivity to dependence via randomization

The above-defined FWER controlling procedures rely on one of the assumptions on the joint distribution of the p-values formulated in Section 2.1.3, namely general dependence, (PRDS(H 0 )) or (indep). By construction, these procedures are not adaptive to the dependency structure at hand in a specific context. To illustrate this point let us reformulate the FWER of a thresholding procedure R = {i ∈ N m , p i ≤ t}:

FWER P (R) = P ∃i ∈ H 0 (P), p i ≤ t = P inf i∈H0(P) p i ≤ t .
Denoting by q α (A) the (1 -α)-quantile of the distribution of inf{p i : i ∈ A}, an optimal choice is t = q α (H 0 (P)). Using this formulation, the above-described FWER-controlling procedures provide lower bounds for this quantile under specific dependency assumptions.

Rather than making such assumptions, Westfall and Young [START_REF] Westfall | Resampling-Based Multiple Testing[END_REF] proposed to use procedures based on permutation, called minP and maxT, to build adaptive bounds for q α (N m ), where q α (N m ) ≤ q α (H 0 (P)). They also introduced step-down versions of these procedures. Romano and Wolf [START_REF] Romano | Exact and approximate stepdown methods for multiple hypothesis testing[END_REF] introduced a randomization assumption under which these procedures are proved to control FWER, bypassing the need for a technical condition called subset pivotality as in the original results of Westfall and Young [START_REF] Westfall | Resampling-Based Multiple Testing[END_REF]. Here, we use a recent formulation of this randomization assumption due to Hemerik and Goeman [START_REF] Hemerik | Exact testing with random permutations[END_REF], which is slightly weaker than the assumption of Romano and Wolf [START_REF] Romano | Exact and approximate stepdown methods for multiple hypothesis testing[END_REF]. Specifically, we assume that there exists a finite group of transformations G acting onto the observation space X , is such a way that the joint distribution of the transformed null p-values is invariant under the action of any g ∈ G. Formally,

∀P ∈ P, ∀g ∈ G, (p H0 (g .X)) g ∈G ∼ (p H0 (g .g.X)) g ∈G , (Rand)
where g.X denotes X that has been transformed by g. We refer to [START_REF] Hemerik | Exact testing with random permutations[END_REF] for examples of situations where (Rand) holds. Examples based on sign-flipping in the location model (2.2) and permutation testing for two-sample tests are given in Section 5.

Extension to k-FWER control

We refer to [START_REF] Roquain | Type I error rate control in multiple testing: a survey with proofs[END_REF][START_REF] Lehmann | Generalizations of the familywise error rate[END_REF] for a survey of k-FWER controlling procedures. In particular, [START_REF] Lehmann | Generalizations of the familywise error rate[END_REF] provides single step and step-down generalizations of the adaptive FWER controlling procedures under (Rand) to k-FWER control. Here, we simply note that a natural generalization of the Bonferroni procedure to k-FWER control may be obtained by elementary arguments. Indeed, the k-FWER of the procedure with threshold t may be written as

P i∈H0 1 {p i ≤ t} ≥ k ≤ k -1 E i∈H0 1 {p i ≤ t} = k -1 i∈H0 P (p i ≤ t) ≤ |H 0 | k t ,
where the first inequality holds by Markov's inequality, and the last inequality holds by the definition of the p-values. Therefore, the "Generalized Bonferroni" procedure with threshold t k-Bonf = αk/m controls k-FWER under general dependence. Similar to the improvement of the Bonferroni procedure for FWER control by the Holm procedure, it is possible to obtain a more powerful k-FWER controlling procedure by building a step-down procedure with critical values:

αk m -(i -k)1 {i > k} .
Just as for FWER control, the threshold of this procedure depends on the p-values, contrary to the deterministic threshold of the generalized Bonferroni procedure.

False

(α) = max i ∈ N m , p (i) ≤ αi m . (2.13)
The BH procedure controls FDR at level α under independence [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF]. Benjamini and Yekutieli [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF] later proved that the BH(α) procedure controls FDR at level π 0 α under (PRDS(H 0 )), which is a much weaker assumption than (indep). Moreover, the BH(α/C(m))

procedure controls FDR at level π 0 α under general dependence [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF], where C(m) is Hommel's correction factor for dependence defined in Proposition 2.3. This procedure is called the BY procedure at level α. Although the BH procedure is based on Simes critical values, the results obtained in [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF] are not a consequence of the Simes inequality. Sarkar et al. [START_REF] Sarkar | Beyond parametrics in interdisciplinary research: Festschrift in honor of Professor Pranab K. Sen[END_REF] noted that the validity of the Simes inequality is in fact a necessary condition for the BH procedure to control FDR.

Estimation of the proportion of true null hypotheses

The FWER control by the Bonferroni procedure and the control offered by the Simes or Hommel inequalities are based on the behavior on the true null hypotheses, and therefore imply the unknown factor π 0 . Similarly, the FDR control provided by the BH procedure at the target level α is conservative by a factor π 0 ≤ 1. Several procedures have been proposed to fill the "conservativeness gap" between FDR control at π 0 α and α. Such procedures incorporate an explicit or implicit estimation of the proportion π 0 of true null hypotheses [START_REF] Benjamini | Adaptive linear step-up procedures that control the false discovery rate[END_REF][START_REF] Gavrilov | An adaptive step-down procedure with proven FDR control under independence[END_REF][START_REF] Genovese | A Stochastic Process Approach to False Discovery Control[END_REF][START_REF] Storey | Strong control, conservative point estimation and simultaneous conservative consistency of false discovery rates: a unified approach[END_REF][START_REF] Storey | A direct approach to false discovery rates[END_REF]. One important class of such procedures is plug-in procedures, which apply the BH procedure at level α/ π0 , where π0 is an estimator of π 0 . We give the definition of one of the most widely used estimator of π 0 in applications is the Storey-λ estimator [START_REF] Storey | Strong control, conservative point estimation and simultaneous conservative consistency of false discovery rates: a unified approach[END_REF][START_REF] Storey | A direct approach to false discovery rates[END_REF]:

π0 (λ) = 1 + 1/m -G m (λ) 1 -λ . (2.14)
The FDR control of plug-in procedures such as the Storey-λ typically only holds under (indep). Some asymptotic properties of such plug-in procedures under (indep) are studied in Section 3.

FDR control vs FWER control

As the FDP is between 0 and 1, we have FDR = E(FDP) ≤ P(FDP > 0) = FWER. Therefore, any procedure controlling FWER at a given level also controls FDR at the same level. In the multiple testing literature, FDR controlling procedures are often said to be "more powerful" than FWER controlling procedures, because FDR controlling procedures generally reject more hypotheses than FWER controlling procedures. As an illustration, for the Leukemia study described in Chapter 1, 20 genes are called differentially expressed at a FWER level of 0.05 (Bonferroni, Holm, Hochberg and Hommel procedures produced the same results for this data set); 163 genes (including these 20) are called differentially expressed at a FDR level of 0.05 by the BH procedure, 30 of them being also rejected by the BY procedure at the same level. We believe that the above interpretation in terms of "power" is quite misleading, because it does not make sense to compare the number of rejections (or power) of two statistical procedures that aim at controlling different type I error risks. By design, FWER control is adapted to confirmatory analyses, where a statement on the rejected hypotheses needs to be made with high probability, whereas FDR control has been developed for exploratory analyses, where a statement in expectation is acceptable.

Perhaps surprisingly, it may be the case that a FWER controlling procedure rejects more hypotheses than a FDR controlling procedure. As noted by [START_REF] Goeman | Multiple hypothesis testing in genomics[END_REF], in the very sparse situation where the expected number of false null hypotheses is o(log(m)), the Bonferroni (or Holm) procedure may reject more hypotheses than the BY procedure, because the log(m) first critical values of the BY procedure are smaller than the Bonferroni threshold α/m. In this extreme example, if one is not willing to make any assumption on the dependency within the p-value family, it may make sense to use the Bonferroni (or Holm) procedure rather than the BY procedure in order to control FDR. This example is mainly an illustration of the conservativeness of the BY procedure, which is seldom used in practice.

FDR control vs FDX control

The great popularity of FDR control and of the BH procedure in particular can probably be explained in part by the fact that the FDP is an appealing and intuitive quantity. An obvious caveat from the statistical perspective is that controlling FDR = E(FDP) does not tell much about the distribution of the underlying FDP. Here we illustrate this point numerically in the model of Example 2.7. One simulation run consists in m = 1, 000 null hypotheses distributed as Gaussian equi-correlated with a given correlation ρ. Among those, m 1 = 200 (corresponding to π 0 = 0.8) are true alternatives, with marginal distribution N (2, 1) while the other m 0 = 800 are true nulls with marginal distribution N (0, 1). For each ρ ∈ {0, 0.1, 0.2, 0.3, 0.4}, we have performed B = 1, 000 simulation runs, and apply the BH procedure at level α = 0.25. The results are summarized in Figure 2.1, where the empirical distribution of the FDP actually achieved on the B replications is represented as a "violin plot", ie a mirrored and rotated kernel density plot. The empirical FDR achieved is represented by a triangle, and the median FDP by a diamond. The oracle BH procedure at level α under independence (ρ = 0) achieves a FDR of π 0 α = 0.2, which is represented by a dashed horizontal line. As expected in this PRDS setting, FDR is controlled in all situations by the BH procedure. When ρ = 0 (independence), the distribution of the FDP of the BH procedure has a Gaussian-like distribution centered at π 0 α and well-concentrated 6 . In this regime, forgetting the random nature of the FDP of the BH procedure is not too problematic, as most of the mass is below the target level α. However, for larger values of ρ the FDP of the BH procedure has a strikingly different empirical distribution. First, FDR control is conservative, in the sense that the FDR achieved is substantially below π 0 α = 0.2, and the distribution of the FDP is heavily shifted toward smaller values. More problematic is the influence of ρ on the dispersion of the FDP achieved by the BH(0.25) procedure: for ρ > 0 the range of the FDP achieved by the BH(0.25) procedure is pretty much the entire [0, 1] interval, with more than 20% of replications having a FDP larger than α = 0.25. In such a situation of positively dependent tests, which is expected to be common in applications, we believe that FDR control is not very helpful, and can be misleading. Indeed, most practitioners assume that FDR ≤ α guarantees that the proportion of false positives in their particular experiment is no more than α. This experiment shows that this may be far from true, even under moderate positive dependence. Although illustrated with the BH procedure for simplicity, this limitation is intrinsic to FDR control, and not specific of a particular FDR controlling procedure.

From this perspective, it seems that a more sensible objective than FDR control would be to control the FDX, that is, quantiles of the FDP instead of only controlling its expectation. A straightforward way to formalize the connection between FDR and FDP control has been proposed by [START_REF] Lehmann | Generalizations of the familywise error rate[END_REF]: by Markov's inequality, for any q ∈ (0, 1), FDX ≤ FDR(R)/q (2.15)

As a consequence, any FDR controlling procedure can be turned in a FDX-controlling procedure. For example, under Assumption PRDS(H 0 ), the BH procedure at level qα yields FDR ≤ qπ 0 α ≤ qα, so that FDX R BH(qα) ) ≤ α .

This simple remark is interesting in that it explicitly provides a quantification of the price to pay for moving from the control of FDR, the expected FDP, to the control of the corresponding tail probability: specifically, for the FDP to be less than some q with probability greater than 1 -α, one needs to apply the FDR controlling procedure at the smaller level qα. We believe this type of statement to be much more informative than FDR control. However, it is also more demanding, in the sense that qα will typically be very small if one desires a guarantee with small FDP (small q) and high-confidence (small α). Continuing the example of DGE studies, taking α = 0.1 and q = 0.5, we can guarantee that FDX(R BH(0.05) ≥ 0.5) ≤ 0.1, that is, with 90% confidence, the FDP among the 163 genes selected by the BH(0.05) procedure is less than 1/2. Or, taking q = 0.1 and α = 0.5, we can guarantee that FDX(R BH(0.05) ) ≥ 0.1) ≤ 0.5 .

That is, the median FDP of the 163 genes selected by the BH(0.05) procedure is less than 1/10. Another approach to FDX control is to simply "augment" the set R of rejections of a FWER controlling procedure by any set A of hypotheses whose cardinality is such that |A|/(|A| + |R|) ≤ q [START_REF] Van Der Laan | Augmentation procedures for control of the generalized family-wise error rate and tail probabilities for the proportion of false positives[END_REF]. Related works on the construction of confidence envelopes for the FDP process [START_REF] Genovese | A Stochastic Process Approach to False Discovery Control[END_REF][START_REF] Genovese | Exceedance Control of the False Discovery Proportion[END_REF][START_REF] Katsevich | Simultaneous high-probability bounds on the false discovery proportion in structured, regression, and online settings[END_REF] are discussed in the next section.

From selective inference to simultaneous inference

In this section, we take one step back and look at existing methods in terms of their ability to address the problem of data snooping, in a multiple testing context (Section 2.4.1) and for linear models (Section 2.4.2). In both cases, we distinguish two broad types of approaches pertaining to selective inference, which is described by [START_REF] Taylor | Statistical learning and selective inference[END_REF] as "the assessment of significance and effect sizes from a dataset after mining the same data to find these associations" and to simultaneous inference in the sense of [START_REF] Scheffé | A method for judging all contrasts in the analysis of variance[END_REF], that is, inference valid for any possible look at the data.

Post hoc inference via multiple testing

As noted by Goeman and Solari [START_REF] Goeman | Multiple testing for exploratory research[END_REF] and more recently by Katsevich and Ramdas [START_REF] Katsevich | Simultaneous high-probability bounds on the false discovery proportion in structured, regression, and online settings[END_REF], controlling multiple testing risks such as FDR -and even FDX -provides statistical guarantees on a specific set of hypotheses selected by the procedure. Therefore, there is a substantial gap between the statistical guarantees provided by state-of-the-art multiple testing procedures and the actual needs of practitioners. To illustrate this important point, let us go back to our example of differential expression analyses. The state-of-the-art approach to this problem consists in performing one statistical test of no difference between means for each gene, and to derive a list R of "significant" genes according to a multiple testing criterion, usually the FDR. This list is then typically refined and/or interpreted using prior knowledge on the problem at hand. For example, as smaller effects are generally of a less relevant from a biological perspective, a typical practice is then to only retain those genes whose "fold change" (that is, the ratio of mean expression levels between the two groups) exceeds a prescribed level [START_REF] Cui | Statistical tests for differential expression in cDNA microarray experiments[END_REF]. Another example is to retain only those genes that belong to a specific biological pathway of interest.

In the above examples, multiple testing procedures provide no statistical guarantee as to the number or proportion of false positives in these user-refined gene lists. This illustrates a gap between theory and practice: multiple testing procedures have been built as tools for statistical inference, but they are commonly used as tools for exploratory data analysis (EDA) tools. From a statistician's perspective, these examples are instances of data snooping, which is inherent to current data-driven research. The statistical community should not only warn investigators about the caveats of this type of practice, but also develop inference methods dedicated to this new paradigm. Two recent research areas have been developed within the multiple testing community: False Coverage Rate (FCR) control in the framework of selective inference, and post hoc inference in the framework of simultaneous inference.

The False Coverage Rate (FCR) is a multiple testing criterion inspired by FDR, but specifically dedicated to inference on confidence intervals for selected parameters [START_REF] Benjamini | False discovery rate-adjusted multiple confidence intervals for selected parameters[END_REF]. The FCR is defined as the average proportion of covering intervals among those selected. Benjamini and Yekutieli [START_REF] Benjamini | False discovery rate-adjusted multiple confidence intervals for selected parameters[END_REF] introduced a generic procedure controlling FCR under (indep) or (PRDS(H 0 )). This procedure consists in adjusting the level of each selected marginal confidence interval by a scaling factor depending on the selection. Importantly, it can be applied to any selection rule. This work has then been generalized to the case where families of hypotheses are selected, instead of individual ones [START_REF] Benjamini | Selective inference on multiple families of hypotheses[END_REF].

An important contribution to the field of simultaneous inference is the work of Goeman and Solari [START_REF] Goeman | Multiple testing for exploratory research[END_REF] to construct procedures that provide confidence statements for the number (or proportion) of false positives in such arbitrary, possibly data-driven subsets of hypotheses. Formally, the aim is to find a functional V α satisfying

P ∀S ⊂ N m , |S ∩ H 0 (P )| ≤ V α (S) ≥ 1 -α , (PH α )
The bound V α (S) is an (1 -α)-upper confidence bound on the number of false positives in S. The aim formalized in (PH α ) can be equivalently formulated in terms of a 1 -α lower confidence bound on the number of true positives in S; we focus on V α for conciseness. A bound V α satisfying (PH α ) is called a post hoc bound Goeman and Solari [START_REF] Goeman | Multiple testing for exploratory research[END_REF], as the set of selected hypotheses may be defined by an investigator after "seeing the data". An earlier contribution in this direction is the work of Genovese and Wasserman [START_REF] Genovese | Exceedance Control of the False Discovery Proportion[END_REF], where it was noted that a bound satisfying (PH α ) could be derived from a multiple testing procedure controlling k-FWER. The derivation of confidence envelopes for the process (FDP(t), t ∈ [0, 1]) under (indep) [START_REF] Genovese | A Stochastic Process Approach to False Discovery Control[END_REF][START_REF] Meinshausen | Lower bounds for the number of false null hypotheses for multiple testing of associations[END_REF][START_REF] Katsevich | Simultaneous high-probability bounds on the false discovery proportion in structured, regression, and online settings[END_REF] or arbitrary dependence [START_REF] Meinshausen | False discovery control for multiple tests of association under general dependence[END_REF] is also related to post hoc inference, because these envelopes are uniform in t. A more precise comparison between these approaches is given in Chapter 5. General post hoc bounds have been obtained by Goeman and Solari [START_REF] Goeman | Multiple testing for exploratory research[END_REF] using a construction based on closed testing. This construction takes advantage of so-called non-consonant rejections in the closed testing procedure. It inherits the elegance and genericity of closed testing, but also its heavy computational complexity (see Section 2.2.1). To address this issue, Goeman and Solari [START_REF] Goeman | Multiple testing for exploratory research[END_REF] have introduced computational shortcuts to obtain a computationally efficient post-hoc procedure in the particular case of Simes local tests, thus under Assumption (PRDS(H 0 )).

Inference after model selection

In this section we follow the notation of [START_REF] Berk | Valid post-selection inference[END_REF] and focus on the case of a linear model of the form

y = µ + , (2.16)
where y is a n-dimensional observation, µ a fixed vector in R n and ∼ N (0, σI n ). We consider a n × p fixed design matrix X, whose columns correspond to explanatory variables for µ. It is not necessarily assumed that there exists β ∈ R p such that Xβ = µ. A model corresponds to a subset of selected variables in {1, . . . , p}. For a fixed model M ⊂ {1, . . . , p}, assuming that X M X M is invertible, we define β M = X † M µ as the target of inference, where X † M = (X M X M ) -1 X M . The ordinary least-squares (OLS) estimate of β M is βM = X † M y. By definition, βM is an unbiased estimator of β M . We use the "full model indexing" convention defined in [START_REF] Berk | Valid post-selection inference[END_REF]: for j ∈ M , we denote by j • M ∈ {1, . . . , |M |}) the rank of j in M , and β j•M = E( βj•M ) is the coordinate of the vector β M corresponding to the j-th predictor in X. A candidate confidence interval on β j•M , centered at βM , is denoted by:

CI j•M (K) = βj•M ± K σ v j•M , (2.17)
where

v j•M is the (j • M )-th row of X † M , which is such that βj•M = v j•M y.
Here, σ is assumed to be an estimator of σ that is independent of the estimates βj•M for all M , and such that σ2 ∼ σ 2 χ 2 r /r some r > 0, as discussed e.g. in [START_REF] Berk | Valid post-selection inference[END_REF]. Let K = t r,1-α/2 be the 1 -α/2 quantile of Student's t-distribution with r degrees of freedom. Then CI j•M (K) is a valid (1 -α)-confidence interval on β j•M , in the sense that

P (β j•M ∈ CI j•M (K)) ≥ 1 -α .
In practice however, the model M is not fixed, as it is typically the result M of a data-dependent selection step. This extra layer of randomness implies that CI j• M (K) is not a valid (1 -α)-confidence interval anymore. Constructing valid inference after model selection is recognized as a difficult problem, see e.g. [START_REF] Leeb | On various confidence intervals postmodel-selection[END_REF][START_REF] Leeb | Can one estimate the unconditional distribution of post-model-selection estimators?[END_REF][START_REF] Kabaila | On the large-sample minimal coverage probability of confidence intervals after model selection[END_REF][START_REF] Leeb | Can one estimate the conditional distribution of post-model-selection estimators?[END_REF][START_REF] Leeb | Model selection and inference: Facts and fiction[END_REF]. The numerous contributions to this problem can be categorized into sample splitting, inference conditional on model selection, and inference uniformly over all model selections.

Sample splitting : single sample splitting consists in using part of the data for model selection, and the rest for inference. This idea is probably the oldest method for inference after model selection [START_REF] Cox | A note on data-splitting for the evaluation of significance levels[END_REF][START_REF] Picard | Data splitting[END_REF]. More recently, this idea has been extended to high-dimensional contexts (see e.g. [START_REF] Bühlmann | High-dimensional variable screening and bias in subsequent inference, with an empirical comparison[END_REF][START_REF] Dezeure | High-dimensional inference: Confidence intervals, p-values and R-software hdi[END_REF]5]). In order to address the fact that the results depend on the split, multi-sample splitting methods [START_REF] Meinshausen | P-values for high-dimensional regression[END_REF] have been proposed, that consist in successively applying single sample splitting to B different splits, and then aggregating the resulting p-values. An intrinsic limitation of this kind of methods is the loss of power inherent to splitting.

Inference conditional on model selection : in the context of penalized linear regression, an asymptotic "covariance test" has been introduced by Lockhart, Taylor, Tibshirani, and Tibshirani [START_REF] Lockhart | A significance test for the lasso[END_REF]. For each k along the regularization path of the lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], it tests for signal evidence in the k th predictor in the model, conditionally on the event that k -1 predictors are already present in the model.

This theory has been extended an consolidated, notably with a test for more generic sequential selection rules [START_REF] Tibshirani | Exact Post-Selection Inference for Sequential Regression Procedures[END_REF], a non-asymptotic "spacing test" [START_REF] Taylor | Inference in adaptive regression via the Kac-Rice formula[END_REF], and a test allows inference to be made at any fixed value of the regularization parameter [START_REF] Lee | Exact post-selection inference, with application to the lasso[END_REF] among other contributions. The power of (a studentized version of) the spacing test is also studied in [START_REF] Azaïs | Power of the spacing test for leastangle regression[END_REF]. By construction, the resulting confidence intervals are conditional on the outcome of a model selection step of the form CI j• M (K), and satisfy guarantees of the form:

P β j•M ∈ CI j• M (K)|M = M ≥ 1 -α .
A recent simulation-based study suggests that the expected length of the associated confidence intervals is typically large and can even be infinite [START_REF] Kivaranovic | Expected length of post-model-selection confidence intervals conditional on polyhedral constraints[END_REF].

Inference uniformly over all model selections : the post hoc adjustment method of Scheffé [START_REF] Scheffé | A method for judging all contrasts in the analysis of variance[END_REF] yields confidence intervals for predictors that are valid for all possible model selections in low-dimensional linear models. This method was introduced in the specific context of analysis of variance (ANOVA) for group mean comparison [START_REF] Scheffé | A method for judging all contrasts in the analysis of variance[END_REF], which can be cast as a specific linear model (2.16). Scheffé's method consists in finding a constant K such that the associated confidence intervals (2.17) are simultaneously valid for all possible contrasts involving group means7 . [START_REF] Scheffé | A method for judging all contrasts in the analysis of variance[END_REF]'s method can be seen as a method to control the FWER, where individual tests would correspond to (infinitely many) contrasts. This method is known as "post hoc analysis" and the type of guarantee it provides is very similar in spirit to the post hoc goals of Section 2.4.1, hence the name "post hoc inference" coined by [START_REF] Goeman | Multiple testing for exploratory research[END_REF] for the latter approach.

Elaborating on this idea of casting the problem of inference after model selection as a problem of simultaneous inference, Berk, Brown, Buja, Zhang, Zhao, et al. [START_REF] Berk | Valid post-selection inference[END_REF] have proposed so-called Post-Selection Inference (PoSI) confidence intervals of the form (2.17), which are simultaneously valid for all possible model selections, that is,

P (∀M ∈ M, ∀j ∈ M, β j•M ∈ CI j•M (K)) ≥ 1 -α ,
where M denotes the set of all possible models, corresponding to all 2 p -1 non-empty subsets of {1, . . . , p}. The corresponding constant K is called the PoSI constant. It depends on the design matrix, which makes the associated confidence possibly sharper than the ones obtained by the Scheffé method, but also much more computationally demanding. Extensions to this work have removed the Gaussian homoscedasticity assumption [START_REF] Bachoc | Uniformly valid confidence intervals post-model-selection[END_REF] and generalized PoSI inference to the problem of prediction after model selection [START_REF] Bachoc | Valid confidence intervals for post-modelselection predictors[END_REF]. Using a reduction of the PoSI problem to a simultaneous inference problem, Kuchibhotla, Brown, Buja, Cai, George, and Zhao [4] have recently introduced confidence intervals of the following form:

CI * M = θ ∈ R |M | , Σ M ( βM -θ) ∞ ≤ C xy (α) + C xx (α) θ 1 ,
where Σ M = n -1 X M X M is the empirical covariance matrix associated with the submodel indexed by M and C xy (α) and

C xx (α) denote (1 -α) joint quantiles of n -1 ( n i=1 X i Y i -E(X i Y i )) ∞ and n -1 ( n i=1 X i X i -E(X i X i )) ∞ .
These confidence intervals enjoy the desired uniform coverage property:

P (∀M ∈ M, β M ∈ CI * M ) ≥ 1 -α .
When compared to the original PoSI intervals, these confidence intervals enjoy two remarkable properties: they have a much reduced computational cost of O(p 2 ), and their volume scales with the size of the selected model and not with the largest possible model size considered.

Contributions

In the next chapters of this part, we describe several contributions to the multiple testing and post hoc inference literature:

• Chapter 3 is motivated by the frequent misinterpretation of FDR controlling procedures in (genomic) applications as "controlling FDP", which is not a well-defined property as the FDP is a random variable. This chapter is a synthesis of my contributions [J24] and [J13], which were done during my PhD thesis and my post-doc.

• Chapter 4 is a synthesis of a joint work Etienne Roquain [J17]] in the framework of classification in sparse high-dimensional models, toward the characterization of the asymptotic properties of FDR controlling procedures viewed as binary classification procedures.

• Chapter 5 describes our contributions to the field of post hoc inference. The contributions described in this chapter [J2, J3, P1] were made in collaboration with Gilles Blanchard and Etienne Roquain in the context of the JCJC ANR project SansSouci (2016-2020). The works described in the last section were performed by Guillermo Durand during the last year of his PhD thesis [START_REF] Durand | Multiple testing and post hoc bounds for heterogeneous data[END_REF]. All of the procedures described in this chapter are implemented in the R package sansSouci [S1], which is available from https://github.com/pneuvial/sanssouci and described at https://pneuvial.github.io/sanssouci.

Before describing these contributions, we briefly mention a contribution to the PoSI literature introduced in Section 2.4.2, which is not described in the rest of this document. Let us denote by K(X, M) the PoSI constant associated with n × p design matrix X, where M denotes the set of all 2 p -1 non-empty submodels of {1, . . . , p}. As mentioned in Section 2.4.2, the computation of K(X, M) for a generic X is a priori of exponential complexity as it involves a maximization over all possible submodels M ∈ M. This motivates attempts to find upper bounds on the PoSI constant. If X is orthogonal, then K(X, M) = Ω( log(p)) [START_REF] Berk | Valid post-selection inference[END_REF]. Moreover, it was shown in an intermediary version of [START_REF] Zhang | Confidence intervals for low dimensional parameters in high dimensional linear models[END_REF] using a cardinality-based argument that the PoSI constant restricted to the set M s of ssparse models satisfies K(X, M s ) = O( s log(p/s)), with no assumption on the design X. In Bachoc, Blanchard, and Neuvial [J6], we have proposed a new upper bound on K(X, M s ) in the case of design matrices satisfying a Restricted Isometry Property (RIP) condition. This upper bound is an explicit function of the RIP constant δ of the design matrix, and it can be seen as an interpolation between the orthogonal setting and the generic sparse setting. In particular, when X is such that δ → 0, our results imply that K(X, M s ) = O( log(p) + δ s log(p/s)). This corresponds to the intuition that for such design matrices, any subset of s columns of X is "approximately orthogonal". Moreover, we have shown that this upper bound is asymptotically optimal in many settings by constructing a matching lower bound.

Chapter 3

Asymptotics of FDR controlling procedures

We introduce a mathematical formalism to study the concentration of the FDP around the FDR asymptotically as the number of hypotheses tested tends to infinity. This formalism is applied to derive central limit theorems for the FDP of a wide class of FDR controlling procedures including the BH procedures and π 0 -adaptive procedures. This chapter is a synthesis of my contributions [J24] and [J13], which has required the adaptation of the original results of [J24] to the "random effects setting" used in [J13]. Moreover, I have chosen to focus here on the results of [J24] obtained for plug-in procedures, because the procedures studied in [J13] also fall into this category. I have also added a section dedicated to FDR control under dependency. 3.1 FDP as a stochastic process of a random threshold

Random effects setting

In the setting defined in Chapter 2, the collection H of hypotheses is deterministic and the p-values (p i ) i∈Nm are such that p i ∼ U[0, 1] when θ i = 0. During my PhD thesis, I studied the asymptotic properties of FDR controlling procedures in this setting [T1, J24].

In this chapter however, we consider the "random effects" setting introduced by [START_REF] Efron | Empirical Bayes Analysis of a Microarray Experiment[END_REF]. This setting, also known as the two-group mixture model, has been widely used in the multiple testing literature, see, e.g., [START_REF] Chi | On the performance of FDR control: constraints and a partial solution[END_REF][START_REF] Genovese | A Stochastic Process Approach to False Discovery Control[END_REF][START_REF] Storey | The Positive False Discovery Rate: A Bayesian Interpretation and the q-value[END_REF], which makes it easier to replace my contributions in the context of this literature. In the present chapter, H = (H 0,i , i ∈ N m ) is a sequence of random indicators, independently and identically distributed as B(1 -π 0 ), where π 0 ∈ (0, 1). Conditional on H, the p-values satisfy p i |θ i ∼ G θi , where G θi denotes the cumulative distribution function of the p-values under θ i . Thus the p-values are independently, identically distributed as

G = π 0 G 0 + (1 -π 0 )G 1
, and m 0 (m) follows the binomial distribution Bin(m, π 0 ). We make the following assumptions on the conditional distributions of the p-values: G 0 is uniform (that is, G 0 (u) = u for any u ∈ [0, 1]), and G 1 is concave. This last assumption is quite natural because the concavity of G 1 is equivalent to the fact that the likelihood ratios of the test statistics under the null versus under the alternative is non-decreasing. We refer to [J13, section 2.1] for a discussion on this assumption.

We consider an asymptotic setting where m → +∞. Therefore, in this chapter the phrase "multiple testing procedure" refers to a collection R = (R m ) m∈N of rejection sets R m rather than a single one as in chapter 2. Consequently, a thresholding procedure is here defined by a collection t = ( tm ) m≥1 of thresholds such that for all m ∈ N, R m = {i ∈ N m , p i ≤ tm }. Following [START_REF] Genovese | A Stochastic Process Approach to False Discovery Control[END_REF], we note that the False Discovery Proportion can be viewed as a stochastic process. Let G 0,m and G 1,m denote the (unobserved) empirical cumulative distribution function of the p-values under the null and alternative hypotheses:

G 0,m (t) = m 0 (m) -1 m i=1 (1 -θ i )1 {p i ≤ t} G 1,m (t) = (m -m 0 (m)) -1 m i=1 θ i 1 {p i ≤ t} .
With this notation, and letting

π 0,m = m 0 (m)/m, G m = π 0,m G 0,m + (1 -π 0,m ) G 1,m is the (observable) empirical cumulative distribution function of the p-values. Moreover, for any t ∈ [0, 1], |R m (t)| = m i=1 1 {p i ≤ t} = G m (t) and |R m (t) ∩ H 0 | = m i=1 (1 -θ i )1 {p i ≤ t} = π 0,m G 0,m (t), so that FDP m (t) = π 0,m G 0,m (t) G m (t) ∨ 1 m (3.1)
is the False Discovery Proportion achieved at the deterministic threshold t. The asymptotic properties of the stochastic process (FDP m (t)) 0≤t≤1 were analyzed by Genovese and Wasserman [START_REF] Genovese | A Stochastic Process Approach to False Discovery Control[END_REF]. The FDR achieved at threshold t, FDR m (t) = E (FDP m (t)), may be written as

FDR m (t) = p(t) (1 -(1 -G(t)) m ), where p(t) = π 0,m t/G(t)
is the positive False Discovery Rate at t defined above in (2.1.2). Using the functional Delta method [START_REF] Van Der | Cambridge Series in Statistical and Probabilistic Mathematics[END_REF], they proved that the FDP m process converges to pFDR at a rate m -1/2 , and built confidence envelopes for the FDP process using this result. However, this result is not sufficient to characterize the behavior of the FDP actually achieved by a given multiple testing procedure, that is, by the random variable FDP m ( tm ). We are interested in the asymptotic behavior of this variable and, in particular, its fluctuations around the asymptotic FDR achieved by the procedure with threshold tm . Following [START_REF] Genovese | A Stochastic Process Approach to False Discovery Control[END_REF], we note that the threshold tBH (α) of the BH procedure defined in Section 2.3 may also be written as Sorted p-values q q q q q q q q q q q q q q q y = x y= αx τ ατ q q False positive False negative q q False positive False negative α=0. Empirical cdf of the p-values q q q q q q q q q q q q q q q y = x y= x α ατ τ ^q q False positive False negative q q False positive False negative α=0.2 ; FDP=3/17 Therefore, defining

tBH (α) = sup u ∈ [0, 1], G m (u) ≥ u/α . ( 3
U(F, α) = sup{u ∈ [0, 1], F (u) ≥ u/α} , (3.3)
we have tBH (α) = U( G m , α). Figure 3.1 illustrates the duality between the classical definition of the BH procedure, and its interpretation using threshold functions. Following [START_REF] Finner | On the false discovery rate and an asymptotically optimal rejection curve[END_REF], u → u/α is called the rejection curve of the BH procedure (also known as the Simes line [START_REF] Simes | An improved Bonferroni procedure for multiple tests of significance[END_REF]).

Two types of plug-in procedures

The plug-in procedures introduced in Chapter 2 are a class of procedures that aim at filling the "conservativeness gap" between the level π 0 α, which is the FDR control actually provided by the BH procedure, and the target FDR level α. Plug-in procedures mimic the Oracle BH(α/π 0 ) procedure by applying the BH procedure at level α/ π0 , where π0 is an estimator of π 0 . Therefore, the threshold of such procedures may be written as

tm (α) = U( G m , α/π 0 ) . (3.4)
We consider two types of such estimators of π 0 : estimators that are based on the empirical cumulative distribution function G m of the p-values, and estimators that are based on the estimation of the density g(1) of the p-values at 1. We begin with the Storey-λ procedure introduced in Section 2.3. The threshold of this procedure may be written as

t0,λ m (α) = U G m , α/Π λ ( G m ) , (3.5)
where Π λ : F → (1 -F (λ))/(1 -λ) only depends on λ1 . More generally, we define plug-in procedures of order 0 (PI-0) as those whose threshold can be written as

t0 m (α) = U G m , α/Π( G m ) . (3.6)
for a given Π :

D[0, 1] → [0, 1],
where D[0, 1] is the set of càdlàg functions on [0, 1], that is, the set of all real-valued functions on [0, 1] that are right-continuous at each point of [0, 1) and have left limits in each point of (0, 1]. In particular, the BH procedure can be viewed as a PI-0 procedure with a constant Π = 1. The BKY procedure [START_REF] Benjamini | Adaptive linear step-up procedures that control the false discovery rate[END_REF] is another example of PI-0 procedure studied in [J24]. One limitation of PI-0 estimators in general and the Storey-λ estimator in particular, is that these estimators are generally not consistent estimators of π 0 . For λ ∈ (0, 1), the expectation of the Storey-λ estimator is

Π λ (G) = π 0 + (1 -π 0 ) 1 -G 1 (λ) 1 -λ , (3.7)
As G 1 is concave, Π λ (G) is a non-increasing function of λ, which is therefore lower-bounded by

π 0 = π 0 + (1 -π 0 )g 1 (1) = g(1)
. Therefore, estimators of g(1) are relevant candidates to estimate π 0 . The second class of estimators considered here is a general class of kernel estimators of g(1) = π 0 , which we have introduced in [J13].

Definition 3.1.

1. A kernel of order ∈ N is a function K : R → R such that the functions u → u j K(u) are integrable for any j = 0 . . . , and satisfy R K = 1, and R u j K(u)du = 0 for j = 1 . . . .

2.

The kernel estimator of a density g at x 0 based on m independent, identically distributed observations x 1 , . . . x m from g is defined by

ĝm (x 0 ) = 1 mh m i=1 K x i -x 0 h ,
where h > 0 is called the bandwidth of the estimator, and K is a kernel.

The threshold of the corresponding FDR controlling procedures may be written as

t1 m (α) = U G m , α/ĝ m (1) . (3.8)
We call procedures whose threshold can be written as (3.8) plug-in procedures of order 1 (PI-1), because their threshold depends on the observations through both G m and ĝm (1), where g is the derivative (of order 1) of G. Note that Storey's estimator may technically be viewed as a kernel estimator of order 0, with kernel function

K(t) = 1 {[-1, 0]} (t).

Criticality

The notion of criticality has been introduced by [START_REF] Chi | On the performance of FDR control: constraints and a partial solution[END_REF] in the context of the BH procedure.

It is necessary to define this notion here because it is tightly connected to the asymptotic properties of FDR controlling procedures. Letting α = lim u→0 u/G(u), if α < α , the number of discoveries made by the BH procedure is stochastically bounded as the number of tested hypotheses increases. Conversely, if α > α , the proportion of discoveries converges in probability to a positive value. This property is quite intuitive when we recall that according to . This restriction is intrinsic to the multiple testing problem, in the sense that it holds regardless of the considered multiple testing procedure. The bound α can be interpreted as the critical value of the abovementioned Oracle BH procedure. Obviously, when α = 0, there is no criticality, e.g. the BH procedure has positive asymptotic power for any α > 0. For instance, this is the case in the Gaussian location model, but not in the Laplace location model, as illustrated in Figure 3.2. More examples of critical and non-critical distributional settings can be found in [J13, Section 5].

The critical value of a generic multiple testing procedure may be defined as the largest target FDR level for which the corresponding asymptotic power is null [J13, Definition 11].

The following Lemma gives a lower bound on the critical value of plug-in procedures.

Lemma 3.2 (Criticality of plug-in procedures [J13]). Let α m be a sequence of (possibly datadependent) levels that converges in probability to α ∞ ∈ (0, 1) as m → +∞. If α ∞ < α , then the threshold tm (α m ) of the BH(α m ) procedure converges in probability to 0 as m → +∞. If the convergence of α m to α ∞ holds almost surely, then the convergence of tm (α m ) to 0 holds almost surely as well.

This property holds in particular for both PI-0 and PI-1 procedures. In the next sections, we give results specific to each type of procedure.

Results for PI-0 procedures

The asymptotic distribution of the FDP of a wide class of procedures including PI-0 procedures has been obtained in [J24] in the original setting of Benjamini and Hochberg [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF]. A fundamental result to derive such asymptotic distributions is given in the next Theorem. Theorem 3.3 (Adapted from [J24]). Consider a multiple testing procedure with threshold function T . Assume that:

(i) T (G) > 0, and is Hadamard-differentiable2 at G;

(ii) (π 0,m G 0,m , (1 -π 0,m ) G 1,m ) converges in distribution to (π 0 G 0 , (1 -π 0 )G 1 ) at rate r m as m → +∞, where r m → 0;
Then the FDP of the procedure, FDP m (T ( G m )), converges in distribution to π 0 τ /G(τ ) at rate r m as m → +∞.

The limit distribution is an explicit function of the limit distribution Z of (π 0,m G 0,m , (1π 0,m ) G 1,m ) and the Hadamard derivative at G taken at Z [J24, Theorem 3.2]. The proof relies on the functional Delta method [START_REF] Van Der | Cambridge Series in Statistical and Probabilistic Mathematics[END_REF]. The interest of this result is that it breaks down the asymptotic analysis of the FDP of a multiple testing procedure with threshold function T into (i) the Hadamard differentiability of the mapping T , which only depends on the procedure, and (ii) the asymptotic behavior of the empirical distribution functions of the p-values, which only depends on the model assumptions. We will go back to this modularity in Section 3.4.

In particular, as for (ii), (π 0,m G 0,m , (1 -π 0,m ) G 1,m ) classically satisfies a Donsker-type theorem (with r m = m -1/2 ) in the random effects setting, as noted by [START_REF] Genovese | A Stochastic Process Approach to False Discovery Control[END_REF]Theorem 4.1]. Moreover, in the case of a PI-0 procedure associated to Π, (i) holds as soon as Π is itself Hadamard-differentiable at G, tangentially to C[0, 1], and satisfies α > Π(G)α [J24, Proposition 7.11 and Corollary 7.12]. Under these assumptions, Theorem 3.3 yields the convergence in distribution as m → +∞ of

√ m FDP m ( t0 m (α)) -π 0 α/Π(G) , where t0 m (α) = U( G m , α/Π( G m ))
is the threshold of the PI-0 procedure associated to Π. The Hadamarddifferentiability of Π is a technical condition which is typically verified for classical PI-0 procedures such as the Storey-λ procedure or the BKY procedure. The important assumption is the condition: α > Π(G)α . Recalling that by the definition of Π, the PI-0 procedure consists in applying the BH procedure at level α m = α/Π( G m ), the above result combined with Lemma 3.2 implies that Π(G)α is the critical value of the PI-0 procedure associated to Π.

As both the BH and the Storey-λ procedures can be viewed as PI-0 procedures, the asymptotic distribution of the FDP of these procedures are obtained below as consequences of [J24, Theorem 4.12]; these results were also stated in the Appendix of [J13].

Corollary 3.4 (BH procedure -adapted from [J24]). Let tm (α) = U( G m , α) denote the threshold of the BH(α) procedure, and t ∞ (α) = U(G, α). For any α ≥ α , we have

√ m FDP m ( tm (α)) -π 0 α N 0, (π 0 α) 2 1 π 0 t ∞ (α) -1 (3.9)
This central limit theorem provides a theoretical explanation for the Gaussian shape of the empirical distribution of the FDP observed in the numerical experiments reported in Section 2.3.4 the particular case of independence (corresponding to ρ = 0 in the Gaussian equi-correlated model of Example 2.7).

Corollary 3.5 (Storey-λ procedure -adapted from [J24]). For any λ ∈ [0, 1), and α ∈

[0, 1], let Π λ : F → (1 -F (λ))/(1 -λ). Let t0,λ m (α) = U( G m , α/Π( G m
)) be the threshold of the Storey-λ procedure at level α, and

t λ ∞ (α) = U(G, α/Π(G)). Then, for any α > Π λ (G)α , we have: √ m FDP m ( t0,λ m (α)) -π 0 α/Π λ (G) N 0, σ 2 λ , (3.10)
where

σ 2 λ = π 0 α Π λ (G) 2 1 π 0 t λ ∞ (α) + 2 t λ ∞ (α) ∧ λ t λ ∞ (α)(1 -G(λ)) - 1 1 -G(λ)
Note that Corollary 3.5 with λ = 0 recovers Corollary 3.4. Moreover, the asymptotic properties of the BH Oracle procedure can be obtained by applying Corollary 3.4 at level α/π 0 .

Results for PI-1 procedures

Corollary 3.5 establishes that the FDP of the Storey-λ procedure converges in distribution at the parametric rate m -1/2 to π 0 α/Π λ (G). As noted above, we have Π λ (G) > g(1) ≥ π 0 . In this section, we show that PI-1 procedures, which are based on estimators of g(1), may provide an asymptotic FDP control closer to the target FDR level α. However, this tighter control comes at the price of slower convergence rates for the FDP. Indeed, PI-0 estimators involve estimators of π 0 that depend on the observations only through G m , and therefore achieve parametric convergence rates. Conversely, PI-1 estimators involve estimators of π 0 that rely on the estimation of g(1), and such estimators are known to yield non-parametric convergence rates that depend on the regularity of g at 1. For a kernel estimator of g(1), the convergence rate is typically of m -k/(2k+1) when g is k times differentiable at 1 [START_REF] Tsybakov | Introduction to Nonparametric Estimation[END_REF]. We state in Proposition 3.6 the asymptotic properties of a particular kernel estimator: the Storey-λ estimator with λ = 1 -h m tending to 1 as m → +∞. This estimator is denoted by π0,m (1 -h m ).

Proposition 3.6 (Asymptotic properties of π0,m (1 -h m ) -[J13]). Assume that g is k times differentiable at 1 for k ≥ 1, with g (l) (1) = 0 for 1 ≤ l < k.
1. If g (k) (1) = 0, then the asymptotically optimal bandwidth for π0,m (1 -h m ) in terms of Mean Squared Error (MSE) is of order m -1/(2k+1) , and the corresponding MSE is of order m -2k/(2k+1) .

2. Let η m be any sequence such that η m → 0 and m k/(2k+1) η m → +∞ as m → +∞.

Then, letting h m (k) = m -1/(2k+1) η 2 m , we have, as m → +∞:

m k/(2k+1) η m (π 0,m (1 -h m (k)) -g(1)) N (0, g (1)) (3.11)
This result cannot be derived from classical results on kernel estimators (e.g. [START_REF] Tsybakov | Introduction to Nonparametric Estimation[END_REF]) as such results typically require that the order of the kernel matches the regularity k of the density, whereas the kernel of Storey's estimator, K(t) = 1 {[-1, 0]} (t), is of order 0. The proof is based on a formulation of π0,m (1-h m ) as a sum of m independent random variables that satisfy the Lindeberg-Feller conditions for the Central Limit Theorem [START_REF] Pollard | Convergence of Stochastic Processes[END_REF].

In order to derive the asymptotic distribution of the FDP of PI-1 procedures, we begin by stating a slightly more general result which holds for any plug-in procedure of the form BH(α/π 0,m ), where π0,m is a generic estimator of π 0 whose convergence rate is slower than m -1/2 . Theorem 3.7 (Plug-in procedures with non-parametric rates - [J13]). Let π0,m be an estimator of π 0 whose asymptotic distribution is given by

mh m (π 0,m -π 0,∞ ) N (0, s 2 0 )
for some s 0 , with h m = o(1/ log log m) and mh m → +∞ as m → +∞. Denote by tm (α) the threshold of the associated plug-in procedure. Then, for any α > π 0,∞ α , The asymptotic distribution of the FDP achieved by the BH(α/π 0,m ) procedure is given by

mh m FDP m ( tm (α)) - π 0 α π 0,∞ N   0, π 0 αs 0 π 2 0,∞ 2   . Corollary 3.8 (Asymptotic distribution of the FDP of PI-1 procedures [J13]). Assume that g is k times differentiable at 1 for k ≥ 1. Define h m (k) = m -1/(2k+1) η 2 m
, where η m → 0 and m k/(2k+1) η m → +∞ as m → +∞. Denote by πk 0,m one of the following two estimators of π 0 :

• Storey's estimator πSto 0,m (1 -h m (k))
; in this case, it is further assumed that g (l) (1) = 0 for 1 ≤ l < k;

• A kernel estimator of g(1) associated with a k th order kernel with bandwidth h m (k).

Then α 0 = g(1)α BH is the critical value of the BH(α/π k 0,m ) procedure, and for any α > α 0 ,

m k/(2k+1) η m FDP m ( τ 0 m (α)) - π 0 α g(1) N 0, π 2 0 α 2 g(1) 3 .

Extensions to other dependency settings

The results obtained in [J24] give a characterization of the asymptotic distribution of the FDP actually achieved by a wide class of multiple testing procedures. In particular, these results make it possible to go beyond FDR control by constructing asymptotic confidence intervals for FDP m ( tm ). The obtained central limit theorems are stated here in the random effects setting of [START_REF] Efron | Empirical Bayes Analysis of a Microarray Experiment[END_REF], and in [J24] in the original setting of [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF]. In both cases, the null hypotheses are assumed to be independent, which is a strong assumption, as discussed in Chapter 2. By Theorem 3.3, the results stated above for the BH and PI(0) procedures may be extended to any multiple testing setting where the asymptotic distribution of (π 0,m G 0,m , (1π 0,m ) G 1,m ) is characterized. For example, Wu [114, Theorem 1] obtained a functional central limit theorem for (π 0,m G 0,m , (1-π 0,m ) G 1,m ) under the assumption that H is stationary and π 0,m itself satisfies a central limit theorem [START_REF] Wu | On false discovery control under dependence[END_REF]Condition 1]. These assumptions cover the case of Ising models in Z 2 .

Two other examples of uses of this formalism to obtain asymptotic results under dependence are Delattre and Roquain [START_REF] Delattre | On the false discovery proportion convergence under Gaussian equi-correlation[END_REF][START_REF] Delattre | On empirical distribution function of high-dimensional Gaussian vector components with an application to multiple testing[END_REF]. Delattre and Roquain [START_REF] Delattre | On the false discovery proportion convergence under Gaussian equi-correlation[END_REF] consider the Gaussian ρ m -equi-correlated model introduced in Example 2.7, additionally assuming that the equicorrelation parameter ρ m satisfies ρ m → 0 as m → +∞. In this model, the empirical process ( G 0,m , G 1,m ) can be shown to converge in distribution at rate r m = min(m, 1/ρ m ) -1/2 [80, Lemma 3.3]. This rate is different from the standard convergence rate m -1/2 holding under independence. A consequence of the results obtained in [J24, J13] is that the FDP of the BH procedure converges to the corresponding false discovery rate (FDR) at the same rate min(m, 1/ρ m ) -1/2 [80, Theorem 2.1]. Using our above-described formalism, the same convergence rates could also be obtained for all of the procedures studied in [J24], including PI(0) procedures.

Delattre and Roquain [START_REF] Delattre | On empirical distribution function of high-dimensional Gaussian vector components with an application to multiple testing[END_REF] have obtained a functional central limit theorem for the process (π 0,m G 0,m , (1 -π 0,m ) G 1,m ) under a more general dependency model, where the test statistics are Gaussian, with a covariance matrix Γ (m) asymptotically lying in a neighborhood of the identity matrix I m3 . The convergence rate of (π 0,m G 0,m , (1 -π 0,m ) G 1,m ) naturally depends on the rate at which the elements of Γ (m) vanish. Interestingly, the asymptotic distribution of the FDP of the BH procedure cannot be deduced directly from the results of [J24] because the joint convergence of (π 0,m G 0,m , (1 -π 0,m ) G 1,m ) is required. Delattre and Roquain [START_REF] Delattre | On empirical distribution function of high-dimensional Gaussian vector components with an application to multiple testing[END_REF] have cleverly noted that the Hadamard derivative of the functional Ψ such that FDP

m ( t) = Ψ(π 0,m G 0,m , (1 -π 0,m ) G 1,m ) at (π 0 G 0 , (1 -π 0 )G 1 )
only depends on the first of its two coordinates. This property makes it possible to extend the original result of [J24] by deriving the asymptotic distribution of FDP of the BH procedure from only the marginal convergence of G 0,m and G 1,m . To our knowledge however, the result obtained in Delattre and Roquain [START_REF] Delattre | On empirical distribution function of high-dimensional Gaussian vector components with an application to multiple testing[END_REF] cannot be extended to PI(0) procedures because the above property of the mapping Ψ is specific to the BH procedures. Indeed, our calculations show that for PI(0) procedures, the Hadamard derivative is given by

α H 0 T (G) -α ΠG (H 0 + H 1 ) Π(G) ,
where the second term (induced by the estimation of π 0 ) does depend on both H 0 and H 1 .

Chapter 4

FDR thresholding for classification under sparsity

Albeit motivated by pure testing considerations, the BH procedure has been shown to enjoy remarkable properties as an estimation procedure in Gaussian location [START_REF] Abramovich | Adapting to Unknown Sparsity by controlling the False Discovery Rate[END_REF] and Laplace scaling [START_REF] Donoho | Asymptotic minimaxity of false discovery rate thresholding for sparse exponential data[END_REF] problems. More specifically, it turns out to be adaptive to the amount of signal contained in the data, which has been referred to as "adaptation to unknown sparsity". More recently, Bogdan, Chakrabarti, Frommlet, and Ghosh [START_REF] Bogdan | Asymptotic Bayes optimality under sparsity of some multiple testing procedures[END_REF] have studied FDR thresholding with respect to the mis-classification risk, and proved that FDR thresholding is asymptotically optimal (as the number m of tests goes to infinity) with respect to that risk in the case of a Gaussian scaling model. A natural question is whether this asymptotic optimality is specific to the Gaussian scaling model or whether it also holds in more general settings. We have generalized and extended the work of [START_REF] Bogdan | Asymptotic Bayes optimality under sparsity of some multiple testing procedures[END_REF] by studying the asymptotic properties of the BH procedure with respect to the mis-classification risk in more general models that encompass Subbotin location and scale models.

θ i = 1 is given by G 1,m (t) = D(D -1 (t) -µ m ).
In the Subbotin scale model, we have

d 1,m (x) = d(x/σ m )/σ m ,
for some (unknown) scale parameter σ m > 1, and the p-values are obtained as

p i = 2D(|X i |), which yields G 1,m (t) = 2D(D -1 (t/2)/σ m )
. This setting covers the Gaussian location model studied in [START_REF] Abramovich | Adapting to Unknown Sparsity by controlling the False Discovery Rate[END_REF] and the Laplace scale model studied in [START_REF] Donoho | Asymptotic minimaxity of false discovery rate thresholding for sparse exponential data[END_REF] for estimation, and the Gaussian scale model studied in [START_REF] Bogdan | Asymptotic Bayes optimality under sparsity of some multiple testing procedures[END_REF] for classification.

Classification risk

We are interested in classification rules that predict the labels (θ i ) i∈Nm from the sample of test statistics (X i ) i∈Nm . We define the classification risk of such a rule as the expected proportion of misclassified items. In our setting where the distribution of the X i under the null distribution is known, a classification procedure can be identified with a threshold tm ∈ [0, 1], that is, a measurable function of the p-value family (p i , i ∈ {1, ..., m}). The corresponding procedure chooses label 1 whenever the p-value is smaller than tm . The corresponding classification risk may be written in terms of p-values:

R c m ( tm ) = m -1 m i=1 P(p i ≤ tm , θ i = 0) + m -1 m i=1 P(p i > tm , θ i = 1)
In particular, for a deterministic threshold

t m ∈ [0, 1], we have R c m (t m ) = π 0,m t m + π 1,m (1 - G 1,m (t m )). Classically, the rule minimizing R c m corresponds to the threshold t B m = f -1 m (τ m )
, where τ m = π 0,m /(1-π 0,m ) is called the mixture parameter of the model. This rule is called the Bayes rule. It depends on the unknown value of the model parameters. Our goal is to find a classification rule whose risk is "comparable" to R c m (t B m ), in the following sense. A classification rule is said to be asymptotically optimal at rate r m = o(1) if and only if there exists D > 0 such that for large m,

R c m ( tm ) -R c m (t B m ) ≤ D × R c m (t B m ) × r m . (4.1)

Power of the Bayes rule under sparsity

Following [START_REF] Bogdan | Asymptotic Bayes optimality under sparsity of some multiple testing procedures[END_REF][START_REF] Abramovich | Adapting to Unknown Sparsity by controlling the False Discovery Rate[END_REF][START_REF] Donoho | Asymptotic minimaxity of false discovery rate thresholding for sparse exponential data[END_REF], we consider a sparse situation, in the sense that τ m → +∞ as m → +∞. Typically, we consider scenarios where τ m = m β , with 0 < β ≤ 1. To compensate for this weakening of the signal, the distance between distribution of the test statistics under the alternative is allowed to grow as m → +∞, in such a way that µ m → +∞ and σ m → +∞ in the location and scale models, respectively. Formally, we make the general assumption that the power C m = G 1,m (t B m ) of the Bayes rule is bounded away from 0 and 1. This assumption makes the classification problem "just solvable" under the sparsity constraint.

Asymptotic optimality of FDR thresholding

Although we are not directly interested in the control of the FDR or the FDP in this chapter, the BH procedure can be viewed as a classification rule with threshold tm = tBH m (α). We have proved in [J17] that there exists a choice of the target level α m , such that the BH procedure at level α m is asymptotically optimal (with an explicit rate) in terms of the above classification risk. Our analysis of the risk of the BH procedure relies on two key elements: a finite sample Oracle inequality on the excess risk of the BH procedure [J17, Corollary 4.3], and an argument of concentration of the threshold of the BH procedure at level α m . We give a short explanation of the second argument. As already noted above, we have tBH

m (α) = sup u ∈ [0, 1], G m (u) ≥ u/α (3.
2). Therefore, the threshold of the BH procedure at level α m is close to the threshold of its deterministic counterpart where G m is replaced by G:

t m (α) = sup {u ∈ [0, 1], G(u) ≥ u/α} .
The classification procedure with threshold t m (α) controls the pFDR. This suggests that a relevant choice for the target (p)FDR level is the level α opt m satisfying

t m (α opt m ) = t B m (4.2)
The main result can be stated as follows for the location and scale models with Subbotin distribution: When applied to the Gaussian scale model, Theorem 4.1 recovers the asymptotic optimality results obtained by [START_REF] Bogdan | Asymptotic Bayes optimality under sparsity of some multiple testing procedures[END_REF], and extends these results by providing explicit convergence rates. Moreover, our result shows that these theoretical properties are not specific to the Gaussian scale model, but carry over to Subbotin location and scale models. These models include the Gaussian location model studied in [START_REF] Abramovich | Adapting to Unknown Sparsity by controlling the False Discovery Rate[END_REF], and the Laplace scale model studied in [START_REF] Donoho | Asymptotic minimaxity of false discovery rate thresholding for sparse exponential data[END_REF]. Theorem 4.1 suggests an asymptotic choice of α m ∝ 1/(log m) γ for FDR thresholding, but does not prescribe an explicit value for α m for a given m. However the level satisfying (4.2) may be explicitly written as a function of the model parameters:

α opt m = (1 + m -β C m /G -1 1,m (C m )) -1 . (4.4)
Therefore, one possible choice for α m in a given model is α m = α opt m (β 0 , C 0 ), where (β 0 , C 0 ) is an a priori value for the unknown model parameters (β, C m ).

Numerical experiments

We report the results of some of the numerical experiments performed in [J17] the Gaussian location model in order to illustrate the adaptation to unknown sparsity by the FDR classification rule, and to discuss the relevance of the choice of α m proposed above. We quantify the quality of a thresholding procedure using the relative excess risk Bayes0 performs well when the sparsity parameter β is correctly specified, and its performance is fairly robust to C m . However, it performs poorly when β is misspecified, and increasingly so as m increases. The results is markedly different for FDR: FDR thresholding is less adaptive to C m than Bayes0, but much more adaptive to the sparsity parameter β, as illustrated by the fact that the configurations with low relative excess risk span the whole range of β.

E m ( t) = (R c m ( t) -R c m (t B m ))/R c m (t B m ).

Extensions

We have also shown in the Supplementary Materials of [J17] that these results can be extended to a more general model where the density of the test statistics under the alternative is log-concave. Following our work [J17], several papers have studied the performance of the BH procedure for classification under sparsity in other settings. In Frommlet and Bogdan [START_REF] Frommlet | Some optimality properties of FDR controlling rules under sparsity[END_REF], the test statistics are independently and identically distributed as a mixture between N (0, 1) (corresponding to a null hypothesis) and a convolution between a measure ν and the N (0, σ 2 /n) distribution (corresponding to the alternative distribution), with a known σ 2 . Other contributions have studied monotone polynomial tail distributions including the Student's t, Pareto, and Inverse Gamma distributions [START_REF] Li | Bayesian multiple testing under sparsity for exponential distributions[END_REF][START_REF] Tang | Properties of multiple testing procedures for Student's t distributions[END_REF][START_REF] Tang | Bayesian Multiple Testing Under Sparsity for Polynomial-Tailed Distributions[END_REF] 1 . Another interesting contribution in the same area of research is [START_REF] Ghosh | Asymptotic Properties of Bayes Risk of a General Class of Shrinkage Priors in Multiple Hypothesis Testing Under Sparsity[END_REF], where the Gaussian scale models with so-called one-group shrinkage priors. In all these contributions, the authors establish the asymptotic optimality of the BH procedure (without studying the convergence rates). 1 Two of these manuscripts seem to have merged into the third one [START_REF] Tang | Bayesian Multiple Testing Under Sparsity for Polynomial-Tailed Distributions[END_REF].

Chapter 5

Post hoc inference via multiple testing

We introduce a generic framework to perform simultaneous inference, based on the control of a risk called the Joint Error Rate. This framework provides a unified view of post hoc inference methods. We propose two types of JER controlling procedures: procedures that are valid under arbitrary dependence under a randomization assumption, and procedures that are valid under independence but dedicated to the case of locally-structured signals.

Goeman and Solari [START_REF] Goeman | Multiple testing for exploratory research[END_REF] have introduced the appealing notion of post hoc inference as a means to close the gap between the exploratory nature of data-driven research, and to address the relative inflexibility and difficulty of interpretation of the control of classical multiple testing risks such as the FDR. We refer so Section 2.3.4 and Figure 2.1 for an illustration of the latter point. Their proposed procedure rely on the availability of (i) a valid local test for all possible intersection hypotheses, and (ii) a computational shortcut to bypass the exponential complexity of the calculation of the bound.

In this chapter, we propose another construction of post hoc inference procedures, where the post hoc guarantee is obtained, by interpolation, from the control of a new multiple testing risk called the Joint Error Rate (JER). After taking a closer look at the construction by Goeman and Solari [START_REF] Goeman | Multiple testing for exploratory research[END_REF] in Section 5.1, we introduce the notion of JER and its properties regarding post hoc inference (Section 5.2), and make a number of connections with previous approaches. In Section 5.4 we introduce a relatively generic way to obtain sharp JER control by adapting to the unknown dependency in the data. We have tried in this presentation to emphasize the ideas and practical implications of the results of [J2]. Finally, in Section 5.5 we describe an extension of this work to the specific case of structured hypotheses (e.g. in time or space)

State of the art and motivation

To motivate this chapter, we go back to the example of differential gene expression analysis of leukemia samples, introduced in Chapter 1. The state-of-the-art approach to this problem consists in performing one statistical test of no difference between means for each gene, and to derive a list of "significant" genes according to a multiple testing criterion, usually the FDR. As argued in Section 2.4.1, a common practice in the biomedical literature is then to only retain those genes whose "fold change" (that is, the ratio of mean expression levels between the two groups) exceeds a prescribed level [START_REF] Cui | Statistical tests for differential expression in cDNA microarray experiments[END_REF]. This is illustrated by Figure 5.1, where each gene is represented as a point in the (log(fold change), -log(p)) plan. This representation is called a "volcano plot" in the biomedical literature. In this example, 163 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -1. genes were selected by the BH procedure at level 0.05. 151 of these 163 genes have an absolute log fold change larger than 0.3, but because FDR is an expected proportion, it is not stable by selection, so that we have no statistical guarantee on these 151 genes. Following [START_REF] Goeman | Multiple testing for exploratory research[END_REF], our goal in this chapter is to provide statistical guarantees on such userdefined selections, in the form of (PH α ), which was already stated in Section 2.4.1.

P ∀S ⊂ N m , |S ∩ H 0 (P )| ≤ V α (S) ≥ 1 -α . (PH α )
We recall that for P ∈ P, H 0 (P ) denotes the set of (indices of) true null hypotheses satisfied by P , that is, H 0 (P ) = {i ∈ N m : P ∈ H 0,i }, and H 1 (P ) = N m \ H 0 (P ). In words, V α (S) is an (1 -α)-upper confidence bound on the number of false positives in S. Using the closed testing machinery, Goeman and Solari [START_REF] Goeman | Multiple testing for exploratory research[END_REF] have proposed a generic bound:

V GS α (S) = max J ⊂R ct |S ∩ J| = max |I| : I ⊂ S, I / ∈ R ct (5.1)
where the second equality holds because R ct (the set of intersection hypotheses rejected by the closed testing procedure, defined in (2.12)) is stable by the superset operation. This bound exploits non-consonant rejections of the closed testing procedure. A non-consonant rejection is a subset I ⊂ N m such that H I is rejected by closed testing, but no elementary hypothesis H i = H 0,i for i ∈ I is rejected by closed testing. Non-consonant rejections are a waste from the FWER control perspective, because FWER control only retains elementary hypotheses, that is, elements of R ct ∩ H. However, they can be useful for post hoc inference, as we now illustrate by reproducing Figure 1 Goeman and Solari [START_REF] Goeman | Multiple testing for exploratory research[END_REF] have introduced shortcuts to obtain a computationally efficient post-hoc procedure in the particular case of Simes local tests, thus valid under Assumption (PRDS(H 0 )). An improved shortcut in linearithmic time (that is, O(m log(m))) complexity has recently been proposed [START_REF] Meijer | Hommel's procedure in linear time[END_REF]. This shortcut is exact, i.e. its output is the bound (5.1) and not an upper bound of it.

Figure 5.3 illustrates the application of this shortcut with α = 0.1. In the left plot, this bound is applied to the set of 151 genes previously selected, and also to the two subsets consisting of genes with positive and negative log fold change. Note that this application to three sets is perfectly valid, because the bound in (5.8) is post hoc, ie it holds jointly an all sets of genes of interest. With probability larger than 0.9, we have the following guarantees: 79 of the 151 genes are truly differentially expressed; among the subset of 27 genes which are more expressed in BCR/ABL samples, at least 1 is a true discovery, while among the complementary subset of 123 genes which are less expressed in BCR/ABL samples, at least 62 are true discoveries. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -1. The right plot of Figure 5.3 illustrates the use of this shortcut to build a (1 -α)-level confidence envelope for the FDP among the most significant items. For example, we are 90% confident that the FDP of the 163 genes that were selected by the BH procedure at level 0.05 is less than 0.48. This statement can be obtained to the one that we obtained directly from FDX control (as a simple by-product of Markov's inequality) in Section 2.3.4, that is, FDX(R BH(0.05) ≥ 0.5) ≤ 0.1. The fundamental difference is that post hoc bounds yield confidence statements on arbitrary subsets of N m , not only S = R BH(0.05) .

Joint Error Rate

Definition and properties

Our proposed construction of post hoc inference procedures relies on a risk measure that we term Joint (Family-Wise) Error Rate (JER). In this section we give the definition of this risk measure and explain how in can be used to build post hoc confidence bounds.

Definition 5.1. Consider a family R = (R k , ζ k ) k∈N K , for some integer K ≥ 1, with R k ⊂ N m and ζ k ∈ N.
We define the Joint Family-Wise Error Rate (JER) of R as

JER(R) = P ∃k ∈ N K , |R k ∩ H 0 | > ζ k (5.2)
The family R is said to control JER at level α ∈ [0, 1] if:

P ∀k ∈ N K , |R k ∩ H 0 | ≤ ζ k ≥ 1 -α. (5.3)
In the above definition, we allow R to be a data-dependent family with R k (X), ζ k (X), but omit the dependence in X to ease notation. If we denote by

A(R) = {A ⊂ N m : ∀k, |R k ∩ A| ≤ ζ k }, then the functional V * R defined by V * R (S) = max{|S ∩ A| : A ∈ A(R)} , S ⊂ N m (5.4)
satisfies (PH α ) under (5.3). Moreover, it is optimal in the sense that it is the smallest upper bound on |S ∩ H 0 | which is valid for any S under (5.3), because H 0 could be any subset A of N m satisfying ∀k, |R k ∩ A| ≤ ζ k . However, it can be shown that the problem of computing V * R (S) given an arbitrary reference family R and S ⊂ N m , is NP-hard [J2, Proposition 2.3]. To overcome this computational complexity, we introduce the following coarser bound:

V R (S) = |S| ∧ min k∈N K |S ∩ R c k | + ζ k , S ⊂ N m . (5.5)
It is easy to see that V R satisfies (PH α ) when R controls JER. Indeed, we note that on the event of probability greater than 1 -α on which (5.3) holds, we have for any k ∈ N K and any S ⊂ N m ,

|S ∩ H 0 | = |S ∩ R c k ∩ H 0 | + |S ∩ R k ∩ H 0 | ≤ |S ∩ R c k | + |R k ∩ H 0 | ≤ |S ∩ R c k | + ζ k ,
where we have used (5.3) in the last inequality. Proposition 5.2 below formalizes the general link between JER control and the associated post hoc bounds:

Proposition 5.2 ([J2]). Let α ∈ [0, 1]. The family R controls JER at level α if and only if V R or V * R satisfies (PH α ).

Connection to confidence envelopes and FDX control

The JER framework makes it possible to draw interesting connections with the concepts on confidence envelopes [START_REF] Genovese | A Stochastic Process Approach to False Discovery Control[END_REF][START_REF] Genovese | Exceedance Control of the False Discovery Proportion[END_REF][START_REF] Meinshausen | False discovery control for multiple tests of association under general dependence[END_REF][START_REF] Meinshausen | Estimating the proportion of false null hypotheses among a large number of independently tested hypotheses[END_REF], and closed testing [START_REF] Goeman | Multiple testing for exploratory research[END_REF]. Some of these connections are also made in a recent paper by Katsevich and Ramdas [START_REF] Katsevich | Simultaneous high-probability bounds on the false discovery proportion in structured, regression, and online settings[END_REF]. This point is developed in detail in Section S-1 of the supplement of the paper [J2]. Here, we only summarize the main points regarding these connections.

Confidence envelopes and quantile bounding functions. Uniform upper confidence bounds on the empirical distribution function of null p-values are deduced by [START_REF] Genovese | A Stochastic Process Approach to False Discovery Control[END_REF][START_REF] Meinshausen | False discovery control for multiple tests of association under general dependence[END_REF][START_REF] Meinshausen | Estimating the proportion of false null hypotheses among a large number of independently tested hypotheses[END_REF][START_REF] Katsevich | Simultaneous high-probability bounds on the false discovery proportion in structured, regression, and online settings[END_REF] from probabilistic guarantees of the form

P ∀t ∈ [0, 1] : |H 0 ∩ R t | ≤ B(t) ≥ 1 -α ,
where R t = {i : p i ≤ t} denotes a p-value level set. The above guarantee assumption can be interpreted as a specific JER control, and the bounds derived in [START_REF] Meinshausen | Estimating the proportion of false null hypotheses among a large number of independently tested hypotheses[END_REF] can be obtained by V in (5.5) under natural assumptions.

Augmentation. Augmentation procedures [START_REF] Genovese | Exceedance Control of the False Discovery Proportion[END_REF] are based on the control of the k-FWER on a specific set R k , and an interpolation from this control to all possible rejection sets.

The augmentation bound of [START_REF] Genovese | Exceedance Control of the False Discovery Proportion[END_REF] may be written as

V aug (S) = min(|S ∩ R c k | + (k -1), |S|)
. This is a particular case of the bound V in (5.5), when the reference family only consists of a single element (R k , ζ k = k -1). As noted by [START_REF] Genovese | Exceedance Control of the False Discovery Proportion[END_REF], the augmentation approach of [START_REF] Van Der Laan | Augmentation procedures for control of the generalized family-wise error rate and tail probabilities for the proportion of false positives[END_REF] can itself be seen as the particular case when k = 1. Therefore, the bound V in (5.5) can be seen as an extension of the augmentation principle.

Inversion and closed testing. Inversion procedures [START_REF] Genovese | Exceedance Control of the False Discovery Proportion[END_REF] are based on family of local tests for all intersection hypotheses, which also form the basis over which [START_REF] Goeman | Multiple testing for exploratory research[END_REF] is built. The generic post hoc bound obtained in [START_REF] Goeman | Multiple testing for exploratory research[END_REF] is in fact equivalent to the inversion bound that can be derived from Genovese and Wasserman [START_REF] Genovese | Exceedance Control of the False Discovery Proportion[END_REF]Equation (10)], even though the latter relies only on local intersection tests, and not on closed testing.

Moreover, the inversion procedure of [START_REF] Genovese | Exceedance Control of the False Discovery Proportion[END_REF] can be seen as a particular case of the bound V * , for the choice ζ k = |R k | -1 which corresponds to a FWER guarantee in each R k provided by local intersection tests. Conversely, JER control based on a reference family can be embedded as a particular case of the local intersection test framework, by defining the local test of an intersection hypothesis H I as φ I = 0 if and only if

∀k ∈ N K , |R k ∩ I| ≤ ζ k .
To summarize these connections, with the formalism of JER control, the inversion procedure of [START_REF] Genovese | Exceedance Control of the False Discovery Proportion[END_REF], and therefore the post hoc bound of Goeman and Solari [START_REF] Goeman | Multiple testing for exploratory research[END_REF], can be seen as a particular case of the bound V * , while the bound V can be seen as an extension of the augmentation principle [START_REF] Genovese | Exceedance Control of the False Discovery Proportion[END_REF].

Structural assumption on the reference family

In the next sections we further study JER control, under additional assumptions on the reference family R:

• in Sections 5.3 and 5.4, we make two further assumptions: (i)

ζ k = k -1, (ii) R is nested, in the sense that for any k ≤ k , R k ⊂ R k .
In this case, it can be shown [J2, Proposition 2.5] that V R coincides with the optimal bound V * R . In the light of the above connections, this property can be seen as an extension of the equivalence between inversion and augmentation procedures stated in [START_REF] Genovese | Exceedance Control of the False Discovery Proportion[END_REF]Theorem 5].

• In Section 5.5 we take a drastically different point of view in order to address the case where the hypotheses are linearly ordered (e.g. in time or space). We assume that the (R k ) are given and deterministic, and that any two elements of the reference family are either disjoint or nested. The corresponding ζ k are then random and have to be calibrated from the data in order for (5.3) to be satisfied.

JER control via Simes inequality

We assume in this section that R is nested. Under this assumption, it can be shown that the bounds V and V * coincide, i.e. that the bound V is optimal. Moreover, we assume that

ζ k = k -1.
In this setting, it is natural to consider thresholding-based reference families based on p-value level sets, that is, families of the form

R k = {i ∈ N m : p i ≤ t k } , k ∈ {1, . . . , K}, (5.6)
where

t k ∈ R, 1 ≤ k ≤ K.
In view of Proposition 5.2, the main challenge in order to obtain post hoc bounds is to identify a suitable reference family R that controls JER at some prescribed level α.

Obtaining Simes' shortcut by a particular JER control

JER control is related to p (k:H0) , the k-th smallest value in the set {p i , i ∈ H 0 (P )}: JER(R) = P ∃k = 1, . . . , K ∧ m 0 , p (k:H0) ≤ t k .

(5.7)

The most natural example of a thresholding reference family is the Simes threshold family, defined by R k (α) = {i ∈ N m : p i ≤ αk/m} and ζ k = k-1 for k ∈ N m (here, K = m). A straightforward consequence of the Simes and Hommel inequalities (Proposition 2.3) is that the Simes family controls JER at level α under (PRDS(H 0 )), while the Hommel family R k (α/C(m)) controls JER at level α under general dependence. We recall that C(m) = i∈Nm i -1 . The corresponding post hoc bound according to (5.5) is given by

V R(α) (S) = |S| ∧ min 1≤k≤|S| i∈S 1 p i > αk m + k -1 (5.8)
for α ∈ (0, 1) and S ⊂ N m . Proposition 5.2 entails the following result:

Corollary 5.3. For any α ∈ (0, 1),

1. V R(α) satisfies (PH α ) under PRDS(H 0 ); 2. V R(α/C(m)) satisfies (PH α ) under general dependence.
This result illustrates the potential of JER control as a generic device to obtain post hoc bounds: we have obtained an easy-to-compute post-hoc bound as an elementary consequence of the Simes-Hommel inequality. It turns out that the bound (5.8) is the bound obtained by the Simes shortcut in [START_REF] Goeman | Multiple testing for exploratory research[END_REF]. We refer to [J2, Section S-1.4] for a proof.

In our construction, the problem of finding a JER controlling family replaces the twostep construction of [START_REF] Goeman | Multiple testing for exploratory research[END_REF]: (i) find a valid (and usually computationally prohibitive) post hoc bound under a specific probabilistic assumption, and (ii) identify shortcuts to calculate this bound efficiently. We find the formulation (5.8) to be simple and natural: for each k, adding k -1 to the number of p-values in S larger than αk/m (i.e., not rejected by the generalized Bonferroni procedure for k-FWER control) is a natural upper bound on the number of false positives in S.

Limitations of Simes-based post hoc bounds

We report a small simulation study in the same Gaussian ρ-equi-correlated model as in Section 2.3.4, which illustrates the possible conservativeness of the Simes inequality. We consider a "white" setting, that is, all null hypotheses are true, |H 0 | = m = 1, 000. In Table 5.1, we quantify the conservativeness of the Simes inequality as the ratio of its size (that is, the level actually achieved by the left-hand side in inequality (2.8)) to the target level α. Here, the size is estimated from 1, 000 replications. The fact that the size is very close to the target level for ρ = 0 illustrates the sharpness of the Simes inequality under independence. However, for ρ = 0.4, the achieved level is only 42% of the target level. Similarly, the Hommel inequality is sharp in the sense that there exists a worst-case p-value

Equi

| = m = 1, 000 and α = 0.2.
distribution such that it is an equality, but it is typically quite conservative when applied to a specific p-value distribution. By construction, the bounds obtained in Corollary 5.3 inherit the properties of the Simes and Hommel inequalities. In particular, the Simes bound (5.8) is sharp under independence and conservative under (PRDS(H 0 )). Instead of assuming a specific dependency setting, we propose in the next section a construction of JER controlling families (and associated post hoc bounds) that take into account the dependency between the tested hypotheses. This idea can be seen as a generalization of the notion of adaptivity to dependence by randomization for FWER control, that we reviewed in Section 2.2.2.

Adaptive JER control from a reference family

As in the preceding section, we assume in this section that R is nested, with ζ k = k -1, and we consider thresholding-based reference families.

General construction

We consider a reference family R(λ) of the form (5.6), based on thresholds t k (λ), 1 ≤ k ≤ K, for some functions t k : λ ∈ [0, 1] → t k (λ). Our goal is to choose λ = λ(α) in such a way that JER control (5.3) is satisfied. For example, in the above simulation example, we would like when ρ = 0.4 to use the Simes reference family with λ(α) = α/0.42, even if the datagenerating model was unknown. A template can be seen as a spectrum of curves, parametrized by λ. The theoretical results below are given in terms of a generic, fixed template

t k (λ), λ ∈ [0, 1], k ∈ N K .
When λ is fixed, we refer to t k (λ), k ∈ N K , as thresholds. In our examples and numerical illustrations we will work with the following two templates:

• Linear template: t k (λ) = λk/m, t -1 k (y) = ym/k; • Beta template: t k (λ) =λ-quantile of β(k, m -k + 1), t -1 k (y) = P(β(k, m -k + 1) ≤ y).
An illustration for these templates is provided in

k (λ), λ ∈ [0, 1], 1 ≤ k ≤ K, a functional λ(α, A), α ∈ (0, 1), A ⊂ N m , is called a λ-calibration if it is non-increasing in A and satisfies: ∀α ∈ (0, 1), JER(R(λ(α, H 0 ))) ≤ α.
The question of how to associate a valid λ-calibration to a given template is addressed in the next section. Given such a valid λ-calibration, λ(α, H 0 ) is still not accessible because H 0 is unknown. However, as λ is non-decreasing in A, we have that λ(α, H 0 ) ≥ λ(α, N m ), so that JER(R(λ(α, N m ))) ≤ α as well. Therefore, the template t k (λ(α, N m )) already gives us an operative reference family for JER control. In order to overcome the possible conservativeness this template due to the fact that H 0 N m , we define the following stepdown procedure:

Algorithm 1 General step-down algorithm j ← 0 A (0) ← N m repeat j ← j + 1 λ j ← λ(α, A (j-1) ) A (j) ← {i ∈ N m : p i (X) ≥ t 1 (λ j )} until A (j) = A (j-1) return  = A (j)
By construction, this algorithm leads to a tighter JER control than its single-step version based on t k (λ(α, H)). The proof that it yields JER control at the target level α can be obtained using the classical step-down methodology laid down by [START_REF] Romano | Exact and approximate stepdown methods for multiple hypothesis testing[END_REF].

Randomization-based JER control under general dependence

In [J2] we provide valid λ-calibrations both for the case where the observation follows a translation model with known dependence, and for the case of general dependence under the randomization assumption (Rand) introduced in Section 2.2. Here, we focus on the second case, which is more practically relevant. Under (Rand), there exists a group G of transformations acting on the observation set in such a way that the joint distribution of the transformed null p-values is invariant under the action of any element of G. This assumption is satisfied1 in particular in two important settings:

• location models of the form (2.2) with symmetric noise distribution:

G = {-1, 1}
n is the group of signs, which acts on the observation set by the element-wise product;

• two-sample multiple testing problems: G is the symmetric group of order n (where n is the sample size as in the Leukemia study), which acts on the observation set by permutation of the sample labels. This is the case in the differential gene expression study that we have been using as a working example throughout this part.

Given a G satisfying (Rand), we consider a (random) B-tuple (g 1 , g 2 , . . . , g B ) of G (for some B ≥ 2), where g 1 is the identity element of G and g 2 , . . . , g B have been drawn (independently of the other variables) as i.i.d. variables, each being uniformly distributed on G. Let us consider a deterministic template t k (•), 1 ≤ k ≤ K, and, for short, denote for all

A ⊂ N m , Ψ(X, A) = min 1≤k≤K∧|A| t -1 k p (k:A) (X) . Now introduce the (data-dependent) functional λ(α, A) = max λ ≥ 0 : B -1 B j=1
1 {Ψ(g j .X, A) < λ} ≤ α .

(5.9)

In practice, we can compute this functional easily as λ(α, A) = Ψ ( αB +1) where Ψ

(1) ≤ Ψ (2) ≤ • • • ≤ Ψ (B)
denote the ordered sample (Ψ(g j .X, A), 1 ≤ j ≤ B). We have the following result, where A denotes the output of Algorithm 5.4.1.

Theorem 5.6 (λ-calibration for unknown dependence). Consider any p-value family satisfying (Rand), a deterministic template and the associated reference family R(λ). Then the (data-dependent) functional λ(•, •) defined by (5.9) is a λ-calibration in the sense of Definition 5.5, and both R(λ(α, N m )) and R(λ(α, A)) control the JER at level α.

A related idea has been proposed independently by Hemerik, Solari, and Goeman [7] to build confidence envelopes for the False Discovery Proportion. To illustrate Theorem 5.6, we report (in Section 5.4.3) numerical experiments in a Gaussian location model where we use sign-flipping-based λ-calibration, and (in Section 5.4.4) the results of permutation-based λ-calibration in two-sample testing for the Leukemia data set.

Numerical experiments

We report experiments performed with the linear template in the two-sided Gaussian location model under equi-correlation, in the case of an unknown dependence. We let n be the sample size. The observations (X i,j ) i∈Nm ∈ R m , j ∈ N n are distributed as ρ-equi-correlated, and the test statistics for i ∈ N m is defined as T (X i,j , 1 ≤ j ≤ n) = n -1/2 n j=1 X i,j . We use sign-flipping to approximate the joint distribution of the test statistics under the null hypothesis. Specifically, the sign-flipped observation associated to a vector of signs s ∈ {-1, 1} n is defined as

(s.X) i,j = s j X i,j , i ∈ N m , j ∈ N n .
The location parameter is set to µ i = n -1/2 µ 1 {i ∈ H 1 }, where µ > 0 quantifies the signal-to-noise ratio. In Figure 5.5, the target JER level α and the level π 0 α are represented by horizontal lines. Each color corresponds to a different λ-calibration (or absence of calibration for the Simes reference family):

Simes Single step Step down Oracle α λ(α, N m ) λ(α, A) λ(α, H 0 )
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q ρ = 0 ρ = 0.2 ρ = 0.4 π 0 = 0.8 π 0 = 0.9 The JER is controlled at the target level α in all situations, as expected from Theorem 5.6.

Oracle calibration yields exact JER control, up to sampling fluctuations. As discussed above, the Simes reference family with parameter α yields JER equal to π 0 α under independence (ρ = 0), while it is more conservative under positive dependence ρ > 0. Single-step λcalibration addresses this conservativeness by adapting to the (unknown) dependence: it yields JER control at π 0 α in all settings considered. Finally, as the signal-to-noise ratio µ gets larger, the step-down λ-calibration yields a JER closer to the nominal level α in non-sparse situations (π 0 ∈ {0.8, 0.9}). In a sparse situation (π 0 = 0.99), corresponding to m 1 = 10 true alternative hypotheses, the single-step procedure is already quite sharp and essentially indistinguishable from its Oracle counterpart, so this setting has been omitted from Figure 5.5. at level 5%. The Simes bound ensures that the FDP of this subset is not larger than 0.48. As noted above concerning the BH procedure, we have a priori no guarantee that this bound is valid, because such multiple two-sample testing situations have not been shown to satisfy the PRDS assumption under which the Simes inequality is valid2 . In contrast, the λcalibrated bounds built by permutation are by construction valid here. Moreover, both are much sharper than the Simes bound while the λ-calibrated bound using the linear template is twice smaller, ensuring FDP< 0.23, and even smaller for the beta template with K = 50. The bound obtained by λ-calibration of the linear template is uniformly sharper that the original Simes bound (5.8), which corresponds to λ = α. This illustrates the adaptivity to dependence achieved by λ-calibration. The bound obtained from the beta template is less sharp for p-value level sets S k of cardinal less than k = 120, and then sharper. This is consistent with the shape of the threshold functions displayed in Figure 5.4.

Illustration on the Leukemia data set

To further illustrate the power of λ-calibration to obtain adaptivity to dependence in the classical problem of differential expression studies, Figure 5.7 illustrates the application of the bound derived from the linear template with (single-step) λ-calibration using B = 1, 000 random permutations, with α = 0.1 (that is, the same bound as the dashed red line in Figure 5.6). Again, the obtained guarantees are substantially more informative than the ones obtained from the Simes reference family with no λ-calibration (left plot of Figure 5.3). q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -1.

Spatially-structured hypotheses

In this section we describe a recent work [J3] that illustrates the flexibility of the JER framework introduced in Section 5.2. This contribution stems from a use-case described in [START_REF] Meijer | A region-based multiple testing method for hypotheses ordered in space or time[END_REF]. This paper focuses on applications where the hypotheses (H 0,i ) 1≤i≤m to be tested are linearly ordered (e.g. in time or space). In such situations, of particular interest are interval hypotheses, that is, intersection hypotheses of the form H i:j = i≤k≤j H 0,k . A natural idea in the JER framework is to include such sets in the reference family, in order for the associated post hoc bounds to be particularly sharp on such sets. These candidate sets are deterministic and given a priori. The corresponding (possibly random) ζ i:j achieving JER control have to be calibrated from the data or from probabilistic inequalities. This is in sharp contrast to Sections 5.3 and 5.4, where the ζ k are deterministic and given a priori, while the R k are random, and chosen in such a way that JER is controlled. The JER framework described in Section 5.2 accommodates both settings.

Theoretical results

A convenient assumption to accommodate the case of linearly ordered hypotheses is that the reference family satisfies a forest assumption, that is, that any two elements of {R k } k∈K are either disjoint or nested: Definition 5.7. A reference family R = (R k , ζ k ) k∈K is said to have a forest structure if following property is satisfied:

∀k, k ∈ K, R k ∩ R k ∈ {R k , R k , ∅}, (Forest)
If a reference family satisfies (Forest), it can be shown that there exists a partition of the original hypotheses into interval hypotheses called leaves, such that each element of the reference family may be written as the disjoint union of consecutive leaves [J3, Lemma 5].

When the reference family is not nested, we do not necessarily have V * = V . In particular, under (Forest), the bound V may not adequately capture the hierarchical structure of the forest. A natural extension of V to the non-nested case is given by the bounds:

V q R (S) = min Q⊂K,|Q|≤q k∈Q ζ k ∧ |S ∩ R k | + S \ k∈Q R k , 1 ≤ q ≤ K, S ⊂ N m . (5.10)
By construction, we have

V 1 R = V , V q
R is non-increasing in q, and V q R ≤ V * . The main result of [J3] (Theorem 3.6) is that under (Forest), we have

V * = V R ,
where is the number of leaves in a partition associated to the reference family. The proof of this result is constructive and therefore provides an algorithm to calculate V * . The complexity of the computation of V * (S) for a given S with this algorithm is in O(dm), where d is the depth of the forest. The above construction is generic and yields post hoc bounds tailored to the Forest structure as soon as the reference family controls JER. In particular, this is the case as soon as for each k ∈ N

K P |R k ∩ H 0 | > ζ k < α/K,
that is, as soon as ζ k is an over-estimator of the proportion of true null hypotheses in R k at level α/K. A simple way to achieve this without any dependence assumption is to calibrate ζ k using a Bonferroni test, i.e. to choose

ζ k = i∈R k 1 p i > α K |R k | .
While this calibration can be improved by using a Holm-Bonferroni test instead of a Bonferroni test, we can expect the associated post hoc bounds to be conservative because ζ k is deducted from a FWER control on the hypotheses in R k , which is more demanding than the mere estimation of the number of true nulls in R k .

In [J3], we consider an alternative calibration where ζ k is given by a procedure to directly estimate the H 0 ∪ R K . Specifically, we use an estimator inspired by the Storey-λ estimator (2.14). Formally, we define

ζ k = |R k | ∧ min t∈[0,1) C 2(1 -t) + C 2 4(1 -t) 2 + i∈R k 1{p i > t} 1 -t 1/2 2 , k ∈ K, (5.11)
where C = 1 2 log K α . The [START_REF] Dvoretzky | Asymptotic minimax character of the sample distribution function and of the classical multinomial estimator[END_REF] inequality with the optimal constant of [START_REF] Massart | The tight constant in the Dvoretzky-Kiefer-Wolfowitz inequality[END_REF] yields that the associated reference family controls JER under (indep) [J3, Proposition 4.1].

Application to the Leukemia data set

We illustrate these results with a preliminary application of the above bounds to the Leukemia data set. Our biological motivation is the fact that gene expression activity can be clustered along the genome. The m individual hypotheses are naturally partitioned into 23 subsets, each corresponding to a given chromosome. Within each chromosome, we consider sets of s = 10 successive genes. Hence, we focus on a reference family with the following elements

R c,k = {(k -1)s + 1, . . . , min(ks, m c )}, k ∈ N Kc , c ∈ {1, . . . , 23},
where, in chromosome c, m c denotes the number of genes, K c = m c /s the number of corresponding regions. In addition, for each (c, k) we obtain ζ c,k (X) by (5.11), with R k = R c,k and α = α c /K c . This choice accounts for a union bound over all chromosomes. In this genomic example, (indep) may not hold, so we have no formal guarantee that this bound is valid. Therefore, the results obtained below are merely illustrative of the approach and may not have biological relevance.

We report the results for chromosome c = 19, which contains m c = 626 genes. In this particular case, we obtain trivial bounds ζ c,k (X) = |R c,k | for all k ∈ N Kc . However, non-trivial bounds can be obtained by enriching the reference family by recursive binary aggregation of the neighboring R c,k , as described in detail in [J3]. The total number of elements in the enriched family is less than 2K c . In our example, it turns out that (5.10) combined with (5.11) yields 6 true discoveries in the interval R 17:24 and 1 true discovery in the interval R 53:54 , where we have denoted

R u:v = u≤k≤v R c,k .
This is illustrated by Figure 5.8 where the individual p-values are displayed (on the -log 10 scale) as a function of their order on chromosome 19. The sets R 17:24 and R 53:54 are highlighted in orange, with the corresponding number of true discoveries marked in each region. We obtain a non-trivial bound not because of the large effect of any individual q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 100 200 300 400 500 600 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Figure 5.8: Evidence of locally-structured signal on chromosome 19 detected by the bound (5.11).

gene, but because of the presence of sufficiently many moderate effects. In particular, in the rightmost orange region in Figure 5.8, the distribution of -log 10 (p) is shifted away from 0 when compared to the rest of chromosome 19. In comparison, we obtain trivial bounds V R (R 53:54 ) = |R 53:54 | = 2s and V R (R 17:24 ) = |R 17:24 | = 8s from (5.5) both for the linear or the beta template. These numerical results illustrate the interest of these bounds tailored to situations where the signal is expected to be spatially structured.

Part II

Inference from heterogeneous and ordered genomic data Chapter 6

Overview of contributions

This introductory chapter provides an overview of the contributions detailed in Chapters 7 to 10 and their scientific context. We also provide a short description of contributions that are not detailed further in this part. As announced in Chapter 1, these contributions can be categorized in two broad themes.

• Inference from heterogeneous genomic data: by "heterogeneous", we mean that several levels of biological information are available for the same variables (typically, genes). We focus on the question of identifying relevant features in situations where the sampling units are the same across levels of information or where prior knowledge on the dependency between variables is available in the form of gene networks;

• Inference from ordered genomic data: in this case, an important challenge is to summarize the data locally into biologically meaningful units. This can be addressed by segmenting a genome (or more precisely each of its chromosomes) into successive homogeneous regions.

A guiding principle for my research has been to try to take advantage of these data characteristics in order to build relevant trade-offs between biological interpretation and statistical and computational performance.

Remark: Throughout this part, the number of variables is denoted by p, whereas the number of hypotheses tested was denoted by m in Part I.

Inference from heterogeneous genomic data 6.1.1 Association between different levels of biological information

We focus on the case where different levels of biological information are observed for the same variables (e.g., genes) and the same observations (e.g. biological samples, or patients). For instance, within the Cancer Genome Atlas (TCGA) project, molecular profiling of each cancer sample is typically performed at least at 3 different biological levels: DNA copy numbers, DNA methylation, and gene expression. After summarizing each of these types of information at the gene scale (an important bioinformatic issue which is not covered in this document except in Section 9.2 for DNA copy numbers), the data is in the form of a 3 × n × p-dimensional array, where n is the number of observations and p the number of genes. A natural question in this setting is to identify variables for which the different levels of biological information are associated, in a statistical sense yet to be defined more formally for a specific biomedical question.

In Chapter 7, we tackle the question of identifying genes whose DNA copy number is associated with their expression level, accounting for DNA methylation. The approach proposed in that chapter is a marginal (i.e., gene by gene) feature selection method, which explicitly exploits some of the possible links between these levels of biological information. This approach provides statistical guarantees on the obtained estimates [J14, J8, S5]. This work was initiated while I was a post doc at UC Berkeley, jointly with Antoine Chambaz, then a visiting professor in the Biostatistics Department in the group of Mark van der Laan. It stems from the theory of Targeted Minimal Loss Estimation (TMLE), which has been initiated by van der Laan and Rubin [START_REF] Van Der Laan | Targeted maximum likelihood learning[END_REF].

Incorporating prior knowledge in the form of networks

Genes operate not individually but jointly via biological interaction networks, which are partly known. In order to better understand the mechanisms of cancer development, as well as the response of cancer cells to treatments, biomedical scientists need mathematical tools that exploit this knowledge of functional links between genes. Below we describe two contributions in this field. Only the first one is detailed in a specific chapter.

Chapter 8: Statistical tests on graphs [J16, S6]. One of the most classical application of multiple testing to genomics is the search for differentially expressed genes, as illustrated in Part I. Besides gene expression data, we have at our disposal prior knowledge on functional links between genes, for example via gene regulation networks. However, most existing methods either rely on marginal (that is, gene by gene) tests of differential expression, or reduce the biological information at hand to unstructured gene lists, disregarding the network information. In Chapter 8 we describe a method to perform multivariate tests of differential expression on graphs [J16, S6]. This method relies on a dimension reduction of gene expression data driven by the graph topology. This work was done while I was a post-doc at Berkeley, in collaboration with Laurent Jacob (also a post-doc at Berkeley at that time) and Sandrine Dudoit.

Identification of deregulated genes and potential regulators [J10, C1]. The mechanisms of regulation of normal cells are known to be altered in cancer cells. We have proposed a method to identify transcription factors involved in such gene deregulations [C1]. Transcription factors (TFs) are proteins that control or regulate gene expression, that is, the transcription of DNA into RNA. As any other protein, TFs are coded by genes, which justifies using their expression to quantify their activity. The method proposed in [C1] consists in three steps. First, a reference gene regulatory network that connects transcription factors to their downstream targets is inferred from gene expression data in a steady (or normal) state. This was done by adapting of the LICORN method [J25]1 . In a second step, the behavior of genes in tumor samples is then compared to this reference network in order to detect deregulated target genes [J10], as detailed in the next paragraph. Finally, the ability of each transcription factor to explain these deregulations is quantified via a linear model. The performance of this three-step strategy has been illustrated by numerical experiments on a TGCA breast cancer data set. These experiments show that the information about deregulation is complementary to the expression data, as the combination of the two improves the performance of the supervised classification of cancer samples.

The second step relies on a statistical method to identify deregulated genes from a reference regulation network, and gene expression data [J10]. We have proposed a model based on a regulatory process in which all genes are allowed to be deregulated, and the hidden variables correspond to the status (under/over/normally expressed) of the genes. The model parameters are estimated via a tailored Expectation-Maximization (EM) [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] algorithm. The resulting method infers posterior probabilities of gene deregulation for each (gene, observation) pair, whereas more classical approaches only produce a gene deregulation score common to all observations. This distinctive feature, which is relevant to biologists, comes at the price of a more tricky model inference step, due to the large number of possible states for the latent variables. This difficulty is bypassed by taking advantage of the fact that the Maximization (M) step of the EM algorithm only requires the marginal distributions of the latent variables to be known. We have shown by factorizing the model likelihood that these marginal distributions may be estimated at the Expectation (E) step of the EM by an algorithm of "Belief Propagation" [START_REF] Yedidia | Exploring Artificial Intelligence in the New Millennium[END_REF].

This work was done in the context of a collaboration with former colleagues at Évry: Julien Chiquet and Etienne Birmelé at LaMME for the statistical aspects, and Mohamed Elati (Institute of Systems & Synthetic Biology, iSSB). This collaboration started with the support of CNRS (PEPS BMI 2013, PI Etienne Birmelé), and continued thanks to the support of an INSERM grant (2015-2019, PI Mohamed Elati) called LIONS for "Largescale Integrative approach to unravel the complex relationships between differentiatiON and tumorigenesiS". The work was mainly carried out by Thomas Picchetti [J10] during his PhD thesis and then by Magali Champion during her post-doc [C1].

Inference from ordered genomic data

Genetic information is coded in long strings of DNA organized in chromosomes. Highthroughput sequencing such as RNAseq, DNAseq, ChipSeq and Hi-C makes it possible to study biological phenomena along the entire genome at a very high resolution [START_REF] Reuter | High-throughput sequencing technologies[END_REF]. In most cases, neighboring positions are expected to be statistically dependent. Using this a priori information is one way of addressing the complexity of genome-wide analyses. For instance, it is common practice to partition each chromosome into regions, because such regions hopefully correspond to biological relevant or interpretable units (such as genes or binding sites) and because statistical modelling and inference are simplified at the scale of an individual region. In simple cases, such regions are given (for example genes in gene expression analyses). However, in more complex cases as in Chapters 9 and 10, the regions of interest are unknown and need to be learned from the data.

Recovering the "best" partition of p loci for a given number of classes is a segmentation problem, also known as "multiple change point problem". Such segmentation problems are combinatorial in nature, as the number of possible segmentations of p loci for K change points (or breakpoints) for a given K = 1 . . . p -1 is p-1 K = O(p K ). When K = o(p), the best segmentation for all k = 1, . . . , K in terms of 2 loss can be recovered efficiently in quadratic time and space complexity using dynamic programming. However, this complexity is typically too large for genomic applications, where p ∼ 10 4 -10 6 . Therefore, an important challenge is to build statistically sound methods with sufficiently low algorithmic complexity and that provides biologically interpretable results.

Chapter 9: Statistical inference from DNA copy number data. The simplest scenario where the signals to be segmented are piecewise-constant is described in Chapter 9, which is dedicated to the analysis of DNA copy numbers in cancer studies. In this case, segmentation can be cast as a least squares minimization problem [START_REF] Picard | A statistical approach for array-CGH data analysis[END_REF][START_REF] Hocking | Learning smoothing models of copy number profiles using breakpoint annotations[END_REF]. I have started working on the statistical analysis of DNA copy numbers in cancers before my PhD thesis, when I got hired as a research engineer at the Curie Institute. At that time and during my PhD thesis, I developed data normalization methods [J28, S7], contributed to the development of data analysis pipelines [J27, s6], visualization methods [J26, s5] and to a review of clustering methods [J12]. I also worked on the development of a method to distinguish new primary tumors from true recurrences on the basis of copy-number profiles from a patient suffering from a second breast cancer [J23]. Copy numbers turned out to also have a central part in my postdoc with Terry Speed at UC Berkeley, thanks to a long-standing collaboration with Henrik Bengtsson (now an Associate Professor at UCSF). Together we developed a normalization method [J21, s3], and contributed to the development of another one [J18, s2] as well as segmentation method [J19, s1]. We also coauthored a book chapter on the statistical analysis of DNA copy numbers [B1], and participated to larger research projects in the context of the Cancer Genome Atlas project (TCGA) that funded my post doc [J20, J22, J15] I have continued working on this topic since I came back to France [J11, S4, S2], in particular with Guillem Rigaill and during the PhD thesis of Morgane Pierre-Jean (2013-2016) [START_REF] Pierre-Jean | Development of statistical methods for DNA copy number analysis in cancerology[END_REF], which I co-supervised with Catherine Matias.

Chapter 9 presents some statistical insights I gained from working with DNA copy number data. It is mostly written as a review chapter, with contributions posterior to my PhD thesis mentioned along the way in varying levels of detail: normalization methods in Section 9.2, a review of segmentation method in Section 9.3, and a performance evaluation framework in Section 9.4. Several software contributions are also highlighted.

Chapter 10: Adjacency-constrained clustering for Hi-C and GWAS. This chapter tackles the more general situation of data described via a similarity measure, as in the case for linkage disequilibrium in GWAS, and contact maps in Hi-C studies. The segmentation problem consists in finding a common partition of rows and columns of a matrix of similarity between objects, such that the signal is mostly concentrated in the diagonal blocks resulting from the partition. This type of segmentation problems can be tackled by kernel-based segmentation methods [START_REF] Harchaoui | Retrospective mutiple change-point estimation with kernels[END_REF][START_REF] Arlot | A kernel multiple change-point algorithm via model selection[END_REF]. However, the quadratic time complexity of these methods cannot be improved without making additional assumptions on the kernel [START_REF] Celisse | New efficient algorithms for multiple change-point detection with reproducing kernels[END_REF]. Indeed, for a generic kernel, even computing the loss (e.g., the least square error) of any given segmentation into a fixed number of segments has a computational cost of O(p 2 ).

Chapter 10 summarizes a series of works initiated by the PhD thesis of Alia Dehman (2012-2015) [START_REF] Dehman | Spatial clustering of linkage disequilibrium blocks for genome-wide association studies[END_REF], which I co-supervised with Christophe Ambroise. Motivated by the idea of taking into account the structure induced by linkage disequilibrium (LD) in GWAS, this thesis lead to the development of a method to detect blocks of LD associated to a given phenotype [J9] (Section 10.4). An important component of this method is an algorithm called adjclust dedicated to the above-described problem of segmenting a similarity matrix (Section 10.3), and its application to the detection of LD blocks. Since I moved to Toulouse in 2016, we started with Nathalie Vialaneix (INRA MIA-Toulouse) to work on the problem of detecting Topologically Associated Domains (TAD) from Hi-C data, for which this algorithm is also relevant [J5, S3]. This collaboration was supported by a grant from the Mission for Interdisciplinarity (MITI) at CNRS (SCALES project, 2017(SCALES project, -2019)). A major difference between GWAS and Hi-C from a statistical perspective is that while the LD similarity is a kernel, it is not necessarily the case for Hi-C similarity matrices, possibly leading to reversals in the constrained HAC. This point, which is briefly mentioned in Section 10.2 has been studied by Nathanaël Randriamihamison [J1] during the first year (2018-2019) of his PhD thesis. This thesis is co-supervised by Nathalie Vialaneix, Marie Chavent (Université de Bordeaux and Inria) and myself and funded by a joint INRA/Inria doctoral program.

Chapter 7

Targeted minimal loss estimation

Looking for genes whose DNA copy number is "associated with" their expression level in a cancer study can help pinpoint candidates implied in the disease and enhance our understanding of its molecular bases. Genomic covariates may play an important role in the biological process and should therefore be taken into account. For instance, DNA methylation is known to regulate gene expression. To quantify the association between DNA copy number and expression, accounting for such relevant genomic covariates, we propose to define a new parameter of interest, and build a method to infer it upon the targeted minimum loss-based inference principle (TMLE).

Considering associations between DNA copy numbers and expression levels in genes is not new [START_REF] Pollack | Microarray analysis reveals a major direct role of DNA copy number alteration in the transcriptional program of human breast tumors[END_REF][START_REF] Van Wieringen | Nonparametric Testing for DNA Copy Number Induced Differential mRNA Gene Expression[END_REF][START_REF] Andrews | Multi-Platform Whole-Genome Microarray Analyses Refine the Epigenetic Signature of Breast Cancer Metastasis with Gene Expression and Copy Number[END_REF][START_REF] Sun | Integrated analysis of gene expression, CpG island methylation, and gene copy number in breast cancer cells by deep sequencing[END_REF][START_REF] Louhimo | CNAmet: an R package for integrating copy number, methylation and expression data[END_REF]. In contrast to these earlier contributions, ours does explicitly exploit that DNA copy number measurements feature both a reference level and a continuum of other levels, instead of discretizing them or considering them as purely continuous. Moreover, we naturally handle multi-dimensional, continuous covariates without discretization. We do not need to assume that they are normally distributed, nor that their true effect of DNA copy number on gene expression is linear.

Defining the parameter of interest

Let O = (W, X, Y ) be a generic observation, where W , X and Y are respectively the covariates (e.g., DNA methylation), DNA copy number and expression of a gene of interest in a randomly picked biological sample. Let x 0 be a reference value for X, corresponding to the normal state of 2 DNA copies. We assume that the probability to observe X = x 0 is bounded away from 0 and 1. We assume without loss of generality that x 0 = 0. In the absence of additional solid knowledge regarding the law of O, we decide to focus on the following non-parametric variable importance measure Ψ. It is a mapping from M, the set of all laws compatible with the definition of O, to R given by

Ψ(P ) = arg min β∈R E P (Y -E P (Y |X = 0, W ) -βX) 2 . (7.1)
This parameter is a measure of the importance of X relative to Y , accounting for W . Moreover, it is a non-parametric measure as it is defined regardless of a semi-parametric model of the form Y = β(Xx 0 ) + η(W ) + U , with unspecified η and U such that E P (U |X, W ) = 0). Therefore, the parameter of interest, Ψ(P 0 ), is universally defined no matter what properties the unknown true data-generating distribution P 0 enjoys, or does not enjoy.

Interpretation in terms of risk excess. For all P ∈ M, Ψ(P ) may alternatively be written [J14, Proposition 1] as:

Ψ(P ) = E P {X(θ(P )(X, W ) -θ(P )(0, W ))} E P {X 2 } , (7.2)
where θ(P )(X, W ) = E P (Y |X, W ). Therefore, the functional Ψ may be interpreted as a generalization of the notion of risk excess to a continuous X. Indeed, when X ∈ {0, 1}, we have Ψ(P ) = E P {(θ(P )(x 1 , W ) -θ(P )(0, W ))h(P )(W )}, where h(P )(W ) = P (X = x 1 |W )/E P {X 2 }. That is, Ψ(P ) is a weighted version of the classical risk excess: E P {θ(P )(x 1 , W ) -θ(P )(0, W )}, where the weights are given by h.

Targeted Minimal Loss Estimation

General principle. As Ψ is known, a substitution estimator ψ n = Ψ(P n ) may be derived from any estimator P n of the law P 0 of the observations. We propose to construct one such estimator by relying on the theory of Targeted Minimal Loss Estimation (TMLE), which is described in detail in [START_REF] Van Der Laan | Targeted Learning: Causal Inference for Observational and Experimental Data[END_REF][START_REF] Van Der Laan | Targeted learning in data science: causal inference for complex longitudinal studies[END_REF]. This theory offers a generic principle of iterative estimation which may be summarized as follows. The first step consists in building a substitution estimator ψ 0 n = Ψ(P 0 n ) from an initial estimator P 0 n of P 0 . The estimator ψ 0 n is then used to build an updated estimator of P 0 , denoted by P 1 n , which itself induces an update of the substitution estimator: ψ 1 n = Ψ(P 1 n ). The updating step is at the core of the TMLE approach. It relies on the notion of efficient influence curve, which is pivotal in semi-parametric estimation [START_REF] Van Der | Cambridge Series in Statistical and Probabilistic Mathematics[END_REF]. It can be shown that the parameter Ψ defined in (7.1) is pathwise differentiable: for every P ∈ M, there exists a function ∇Ψ P of O such that for any bounded s : O → s(O) and all |ε| < s -1 ∞ , if we characterize P ε ∈ M by setting dP ε /dP = 1 + εs then

Ψ(P ε ) = Ψ(P ) + εE P {∇Ψ P (O)s(O)} + o(ε).
The map ∇Ψ is called the efficient influence curve of Ψ. Informally, and as it names suggests, ∇Ψ P quantifies the influence of the law P on the estimation of the parameter of interest, and indicates in which direction an estimator of P should be refined in order for the quadratic error of the associated substitution estimator to be improved. under additional conditions on the convergence rates of these features:

√ n(ψ * n -Ψ(P 0 )) N (0, Var P0 (∇Ψ(P 0 )(O))) (7.3) Again thanks to the double-robustness of ∇Ψ, the fact that either the estimator of θ or both the estimators of µ and σ 2 converge at rate n -1/2 is sufficient to guarantee the asymptotic normality of the TMLE. Finally, if the estimators of all features are consistent, then the TMLE is asymptotically efficient, and we can construct an estimate of its asymptotic variance.

Contrast to NP estimation. An obvious substitution estimator of Ψ(P 0 ) is

ψ NP n = arg min β∈R E Pn θn (X, W ) -θn (x 0 , W ) -β(X -x 0 ) 2 ,
an expression derived from (7.1) by substituting the empirical measure P n for P 0 and the Nadaraya-Watson estimator θn (X, W ) of E P0 (Y |X, W ) for it. Under regularity conditions of order on the true conditional expectation, the optimal bandwidth h n for θn satisfies h n = cn -1/(2 +d) where d = 2 is the dimension of (X, W ) [START_REF] Tsybakov | Introduction to Nonparametric Estimation[END_REF]. Now it is possible to characterize a

P 0 ∈ M such that E P0 {ψ NP n -Ψ(P 0 )} = c h n + o(h n ).
In particular, ψ NP n cannot achieve √ n-consistency under that P 0 . We see that ψ NP n suffers from the fact that the bias-variance trade-off, which is at the core of the construction of θn , is optimized for the sake of estimating the infinite-dimensional parameter E P0 (Y |X, W ) whereas we are eventually interested in estimating the one-dimensional parameter Ψ(P 0 ).

Simulations

Because association patterns between copy number, expression and methylation are generally non-linear, setting up a realistic simulation model is a difficult task. We design here a simulation strategy based on perturbations of real observed data structures. Specifically, we focus on the gene EGFR in a glioblastoma multiforme (GBM) data set from The Cancer Genome Atlas (TCGA) project [START_REF] Collins | Mapping the cancer genome[END_REF]111]. This gene is known to be altered in GBM. Figure 7.1(a) represents W = DNA methylation, X = DNA copy number, and Y = gene expression data for one particular gene, EGFR, which is known to be altered in GBM. For this gene, the association between copy number and expression is non-linear, and high methylation levels are associated with low expression levels.

Our simulation strategy implements the following constraints:

• There are generally up to three copy number classes: normal regions, and regions of copy number gains and losses;

• In normal regions, expression is negatively correlated with methylation;

• In regions of copy number alteration, copy number and expression are positively correlated. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q DNA copy number 0.87 0.02 0.04 0.06 0.08 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -1 0 1 2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q DNA copy number 0.80 0.00 0.02 0.04 0.06 0.08 0.10 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -2 -1 0 In both graphics, we represent kernel density estimates (diagonal panels), pairwise plots (lower panels), and report the pairwise Pearson correlation coefficients (upper panels). (a). Real dataset corresponding to the EGFR gene in 187 GBM tumor samples. For 130 among the 187 samples, only DNA copy number and gene expression data were available (circles in lower middle plot). (b). Simulated dataset consisting of n = 200 independent copies of the synthetic observed data structure described in the main text. Note that the constant O X 2 is added to each value of X so that graphics corresponding to real and simulated data can be more easily compared.

-
Monte-Carlo estimator of the asymptotic variance appearing in (7.3), which can be formed because ∇Ψ(P )(O) is available in closed form in our simulation (see [J14, Lemma 7]).

For each of B = 10 3 simulation runs, we record the parameter estimate obtained after k iterations of the TMLE procedure:

ψ k n,b = Ψ(P k n,b ), 1 ≤ b ≤ B .
The results of this simulation are summarized by Figure 7.2 in the case where the features are learned using algorithms based on generalized linear models. This figure provides the empirical distribution of the parameter estimate across the B simulation runs. These results illustrate some of the fundamental characteristics of the TMLE estimator and related confidence intervals: B above, which requires the knowledge of the simulation parameters, this estimator can be obtained only based on the n observations. While our theoretical results do not guarantee that it is safe to estimate the limit variance by (s k n,b ) 2 , the associated confidence intervals provide empirical coverage close to the nominal coverage.

An application to TCGA data

As an illustration, we study a breast cancer data set from The Cancer Genome Atlas (TCGA) Network [START_REF] Tcga) Network | Comprehensive molecular portraits of human breast tumours[END_REF]. We downloaded DNA methylation (W ), DNA copy number (X), and expression (Y ) of 11,314 genes for n = 463 patients2 . The dimension of W is the number of CpG loci in the gene's promoter region, which can vary from one gene to another. Conveniently, our implementation handles multi-dimensional covariates.

In order to quantify the influence of DNA methylation on the strength of association between DNA copy number and gene expression in this data set, we compare our proposed non-parametric variable importance measure Ψ(P ) to its counterpart neglecting W . The latter is a different mapping from M to R given by

F(P ) = arg min β∈R E P (Y -βX) 2 = E P {XY } E P {X 2 } .
A natural estimator of F(P ) is given by f

n = n i=1 X i Y i / n i=1 X 2 i .
Our theoretical results show that (f n , ψ n ) satisfies a central limit theorem. Furthermore, it is possible to estimate the corresponding asymptotic covariance matrix, hence the asymptotic variances of ψ n and of (ψ n -f n ).

We have performed the bilateral test of "Ψ(P ) = F(P )" against "Ψ(P ) = F(P )" for each of the 10,246 genes without missing data. Figure 7.3 presents the (-log 10 ) p-values of this test against the gene's position. A pattern emerges of regions featuring very small pvalues, among which chromosomes 1q, 8, 16q. The pattern is not correlated to the marginal Figure 7.3: Each dot corresponds to the genomic position and (-log 10 ) p-value of "Ψ(P ) = F(P )" against "Ψ(P ) = F(P )" for one of the 10,246 genes without missing data. The chromosomes are delimited by vertical grey lines. The background image represents, gene by gene, the proportions of the 463 samples for which X < 0 (blue), X > 0 (red), and X = 0 (white). distribution of X, represented in the background. We also compute the partial correlation of X and Y given W for each gene (not shown). No pattern emerges. This suggests that the approach we propose may be useful to identify novel regions worthy of interest.

For this data set, the typical run time of the method for a single gene with the default options of the package is 10 seconds on a standard laptop. Obviously, this analysis can be parallelized very easily, as each gene is treated independently of all the other ones.

Chapter 8

Graph-structured two sample tests

In this chapter we consider multivariate two-sample tests of means, where the location shift between the two populations is expected to be related to a known graph structure. An important application of such tests is the detection of differentially expressed genes between two patient populations, as shifts in expression levels are expected to be coherent with the structure of graphs reflecting gene properties such as biological process, molecular function, regulation, or metabolism. For a fixed graph of interest, we demonstrate that accounting for graph structure can yield more powerful tests under the assumption of smooth distribution shift on the graph. We also investigate the identification of non-homogeneous subgraphs of a given large graph, which poses both computational and multiple hypothesis testing problems.

The two-sample Hotelling test for multivariate data

Most approaches to the joint analysis of gene expression data and gene graph data involve two distinct steps. Firstly, tests of differential expression are performed separately for each gene. Then, these univariate (gene-level) testing results are extended to the level of gene sets, e.g., by assessing the over-representation of DE genes in each set based on p-values for Fisher's exact test1 (or a χ 2 approximation thereof) adjusted for multiple testing [START_REF] Beissbarth | GOstat: find statistically overrepresented Gene Ontologies within a group of genes[END_REF] or based on permutation adjusted p-values for weighted Kolmogorov-Smirnov-like statistics [START_REF] Subramanian | Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles[END_REF].

There are two major caveats to this general two-step approach:

• The first step is based on marginal tests of association between one gene and class labels, so the information of the joint distribution of gene expression across genes is lost for the second step;

• The sampling units for the tests in the second step are the genes (instead of the observations/patients), which are expected to have strongly correlated expression measures. As noted by [START_REF] Goeman | Analyzing gene expression data in terms of gene sets: methodological issues[END_REF], this renders the interpretation of the tests at the second step problematic and may lead to a large loss of power or Type I error control when sets of genes have correlated expression.

The first point is illustrated by Figure 8.1 where the gene expression levels of two genes (PIAS2 and UBE2C) are displayed for the 79 samples of the Leukemia data set. Although neither of these genes is differentially expressed (the p-value of the t-test for differential expression of PIAS2 and UBE2C are 0.11 and 0.12, respectively , the pair (PIAS2, UBE2C) can be called differentially expressed using a bi-variate Hotelling test (as formally defined below) with p-value 0.023. This example suggests that direct multivariate differential expression testing of gene sets could be more appropriate than posterior aggregation of marginal gene-level tests. We now give the definition of Hotelling's T 2 -test, the most classical multivariate test of location shift. Let us consider random samples (x 1 ) and (x 2 ) of size n 1 and n 2 drawn from two p-dimensional Gaussian distributions, N (µ i , Σ), i = 1, 2. Assuming that p < n 1 + n 2 -1 and Σ is invertible, Hotelling's T 2 -test statistic is defined by

T 2 = n 1 n 2 n 1 + n 2 (x 1 -x2 ) Σ-1 (x 1 -x2 ),
where xi , i = 1, 2 denote the sample means, and Σ the pooled sample covariance matrix. Up to the scaling factor n 1 n 2 /(n 1 + n 2 ), it corresponds to the squared Mahalanobis norm (∆ 2 (x, S) = x S -1 x) of the sample mean shift x1 -x2 associated to the empirical covariance matrix Σ). Under the null hypothesis H 0 :

µ 1 = µ 2 of equal means, N T 2 follows a (central) F -distribution F 0 (p, n 1 + n 2 -p -1), where N = n1+n2-p-1 (n1+n2-2)p . In general, N T 2 follows a non-central F -distribution F ( n1n2 n1+n2 ∆ 2 (δ, Σ); p, n 1 + n 2 -p -1)
, where the non-centrality parameter is, up to the same scaling factor as above, the Mahalanobis norm of the mean shift δ = µ 2 -µ 1 . We refer to ∆ 2 (δ, Σ) as the distribution shift. In the remainder of this chapter, unless otherwise specified, T 2 -statistics are assumed to follow the nominal F -distribution, e.g., for critical value and power calculations.

Hotelling's test is known to be uniformly most powerful invariant for multivariate normal distributions against global-shift alternatives. However, it is not directly applicable to multivariate differential expression testing of gene sets, for two main reasons:

• first, the typical size of gene sets is of the same order as the typical sample size (dozens of genes/samples). Hotelling's test it is only defined when p < n 1 + n 2 -1 and Σ is invertible. Even under these conditions, its power naturally decreases as p increases due to the increasingly poor conditioning of Σ-1 [START_REF] Bai | Effect of high dimension : by an example of a two sample problem[END_REF]. It is expected that genes from a particular gene set will have correlated expression levels. Continuing the example of Figure 8.1, the total expression level of the two genes (PIAS2 + UBE2C) has a much more significant association to the class labels (t-test p-value equal to 0.006) than the pair (PIAS2, UBE2C)

• second, such direct multivariate tests on unstructured gene sets do not take advantage of information on gene regulation or other relevant biological properties. An increasing number of regulation networks are becoming available, e.g., Gene Ontology (GO; http://www.geneontology.org), Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.genome.jp/kegg) or NCI Pathway Integration Database (NCI graphs; http://pid.nci.nih.gov). Such networks specify, for example, which genes activate or inhibit the expression of other genes. If it is known that a particular gene in a tested gene set activates the expression of another, then one expects the two genes to have coherent (differential) expression patterns, e.g., higher expression of the first gene in resistant patients should be accompanied by higher expression of the second gene in these patients.

In the next sections, we propose multivariate test statistics for identifying differential expression patterns (or, more generally, shifts in distribution) that are coherent with a given graph structure. In a nutshell, we propose to take advantage of the prior information encoded in gene networks or pathways to perform graph-structured dimension reduction, combined with a multivariate test in the low-dimensional space. The example of Figure 8.1 suggests that projecting the two-dimensional data (PIAS2, UBE2C) in the one-dimensional space (PIAS2 + UBE2C) may lead to increased differential expression power.

Graph-structured dimension reduction

Let us consider a network of p genes, represented by a graph G = (V, E), with |V| = p nodes and edge set E. Let δ ∈ R p denote the mean shift, i.e., the vector of differences in mean expression measures for these p genes between the two populations of interest. Suppose we expect the shift δ to be coherent with the graph G, in the sense that it has low energy E G (δ) for a particular energy function E G defined on G. Then, we wish to build a space of lower dimension k p capturing most of the low energy functions. To this end, we start by finding the function that has the lowest possible energy, then the function that has lowest possible energy in the orthogonal space of the first one, up to the kth function with lowest energy in the orthogonal subspace of the first k -1 functions. That is, for each i ≤ k, we define

u i =    arg min f ∈R p E G (f ) such that u i ⊥ u j , j < i. (8.1) If E G is a positive semi-definite quadratic form E G (δ) = δ Q G δ, for some positive semi- definite matrix Q G = U ΛU ,
where U is an orthogonal matrix and Λ a diagonal matrix with elements λ i , i = 1, . . . , p, then the solution to Equation (8.1) is given by the k eigenvectors of Q G corresponding to the smallest k eigenvalues. It is easy to check that these eigenvalues are the energies of the corresponding functions u i , i.e., E G (u i ) = λ i .

Any positive semi-definite matrix can be chosen for Q G , with different choices of Q G leading to different notions of coherence of the expression shift with the network. A classical choice is the graph Laplacian L. Suppose G is an undirected graph with adjacency matrix A, with a ij = 1 if and only if (i, j) ∈ E and a ij = 0 otherwise, and degree matrix D = Diag (A1), where 1 is a unit column-vector, Diag(x) is the diagonal matrix with diagonal x for any vector x, and D ii = d i . The Laplacian matrix of G is then typically defined as

L = D -A or L norm = I -D -1 2 AD -1 2
for the normalized version, leading to energies (ij)∈E (δ i -δ j )

2 and (i,j)∈E δ i / √ d i -δ j / d j 2 , respectively. Note that, in this case, the Laplacian matrix L, energy E, and basis functions u i extend the classical Fourier analysis of functions on Euclidean spaces to functions on graphs, by transferring the notions of Laplace operator, Dirichlet energy, and Fourier basis, respectively [START_REF] Evans | Partial Differential Equations[END_REF].

In the specific case of gene regulation networks, it may be relevant to model negative associations between genes by considering a signed version of the adjacency matrix, where a ij = 1 if gene i activates gene j, and -1 if it inhibits gene j. A signed version of the graph Laplacian is then L sign = D -A, where D = Diag (|A|1) is the degree matrix and |A| denotes the entry-wise absolute value of A. Note that the signed Laplacian is always positive definite, see e.g. [START_REF] Gallier | Spectral theory of unsigned and signed graphs[END_REF]Chapter 5]. As an example, let us consider a simple four-node graph whose signed adjacency matrix is given by

A =     0 1 0 0 1 0 1 -1 0 1 0 0 0 -1 0 0     , (8.2)
The eigenvectors of the corresponding signed Laplacian are represented in Figure 8.2. Following our principle to build a lower dimension space, we use the first few eigenvectors of Q G to obtain orthonormal functions with low energy. The first eigenvector, corresponding to the smallest energy (eigenvalue of zero), can be viewed as a "constant" function on the graph: its absolute value is identical for all nodes, but nodes connected by an edge with negative weight take on values of opposite sign. By contrast, the last eigenvector, corresponding to the highest energy, is such that nodes connected by positive edges take on values of opposite sign and nodes connected by negative edges take on values of the same sign.

Figure 8.3 illustrates the projection of two vectors onto the first two dimensions of the graph decomposition. The first vector is smooth along the graph, in the sense that its entries are more coherent with the graph structure. This is reflected by the magnitude of the coefficients of Therefore, it is essentially preserved by the projection in the first two dimensions. Conversely, the second vector is not smooth along the graph, so it is not preserved by the projection.

Smooth shift

Non-smooth shift Original vector 2d projection Original vector 2d projection

(1.45, -0.04, -0.21, 0.

Example of projection of two vectors along first two dimensions of the graph with adjacency matrix defined in (8.2). Left: a smooth vector along the graph; right: a non-smooth vector along the graph. Bottom: coefficients of the eigen decomposition of each vector.

While we have introduced the idea in the context of gene regulation networks and testing for differential expression, the same dimensionality reduction principle applies to any multivariate testing problem for which the variables have a known structure, as represented by a graph.

Graph-structured two-sample tests

In the remainder of this chapter, we denote by f = U f the coefficients of a vector f ∈ R |V| after projection on a basis U (typically the eigenvectors of a Q G matrix).

For any orthonormal basis U and, in particular, for our graph-based basis, direct calculation shows that

T 2 = T 2 ∆ = n 1 n 2 n 1 + n 2 (x 1 -x2 ) U U ΣU -1 U (x 1 -x2 ),
i.e., the statistic T 2 in the original space and the statistic T 2 in the new graph-based space are identical. More generally, for k ≤ p, the statistic in the original space after filtering out dimensions above k is the same as the statistic T 2 k restricted to the first k coefficients in the new space defined by U : Further numerical experiments reported in [J16] quantify empirically the robustness of our proposed method to a possible mis-specification of k 0 (as k 0 is typically unknown in practice), and to a possible mis-specification of the graph, where the graph used to perform dimension reduction and testing is a noisy version of the graph used for data generation.

T 2 k ∆ = n 1 n 2 n 1 + n 2 (x 1 -x2 ) U [k] U [k] ΣU [k] -1 U [k] (x 1 -x2 ) = n 1 n 2 n 1 + n 2 (x 1 -x2 ) U 1 k U U 1 k U ΣU 1 k U + U 1 k U (x 1 -x2 ), 0
F -distribution F n1n2 n1+n2 ∆ 2 k+l ; k + l, n 1 + n 2 -(k + l) -1 , for n 1 = n 2 =

Differential subgraph discovery

The methods described in this chapter are designed to test the differential expression of known graphs. An important related question is whether we can identify non-homogeneous subgraphs of a known graph. This poses a huge combinatorial problem even for moderately large graphs, as the number of (connected) subgraphs of size k of a graph of size p can be exponential in p and k. Exhaustive search is therefore not feasible in practice, especially for differential expression on gene networks, where p is typically in the dozens or hundreds of genes. To address this bottleneck, we describe in [J16] a branch and bound type algorithm which makes it possible to avoid testing all possible subgraphs. This algorithm relies on a pruning strategy based on an upper bound on the value of the test statistic for any subgraph containing a given set of nodes [J16, Supplement A, Lemma 2]. Two variants of this algorithm are proposed: an exact one, and a quicker, approximate one. The approximation in the latter consists in focusing on the subgraphs whose sample mean shift in the first k components in the graph space has a sufficiently large Euclidean norm. Just like for volcano plots, this approximation is justified in practice by the fact that significant subgraphs corresponding to small mean shifts are typically less relevant from a biological perspective. Finally, we have devised a strategy based on class label permutations to estimate the expected number of type I errors and thereby account for the multiple testing situation created by testing a large number of subgraphs.

Chapter 9

Statistical inference from DNA copy number data After a brief description of DNA copy number signals, this chapter reviews some statistical insights gained from working with DNA copy number data. We highlight the importance of data pre-processing, computational and statistical trade-offs for the design of segmentation methods, and the importance of performance evaluation and benchmarking studies in this context.

Copy-number signals

As noted in Section 1.3, an important issue in cancer research is to estimate the underlying copy number state at each position along the genome of a tumor sample. Formally, we define the copy number state of a tumor at a given genomic locus j as a pair of non-negative numbers (γ j , γ j ), where γ j ≤ γ j , which are respectively the smaller and the larger of the two parental copy numbers at this locus. Figure 9.1 illustrates the most common copy number states observed in cancer samples: normal, that is one copy from each parent (1,1), gain (1,2), deletion (0,1), copy-neutral LOH1 (0,2). By definition we have γ j ≤ γ j , and γ j = γ j + γ j is the total copy number. The quantities γ j and γ j are called minor and major copy numbers, respectively. Note that γ j , γ j , and γ j need not be whole numbers, especially because of the possible presence of normal cells in the tumor sample.

Minor and major copy numbers are not directly observed from genotyping microarray or sequencing assays, but they can be estimated from the data collected by these assays. These data can be summarized by a pair of vectors (θ A j , θ B j ) j where the index j refers to a genomic position along a chromosome. In Figure 9.1 four such positions are represented by horizontal segments. These positions correspond to (bi-allelic) Single Nucleotide Polymorphisms (SNPs), that is, genomic positions where the DNA sequence varies at a substantial rate across individuals of some population. For most SNPs only two letters (out of {A, C, G, T} ) variants are observed. These versions of the genetic sequence are called alleles and arbitrarily denoted by A and B in Figure 9.1. The values θ A j and θ B j correspond to the signal intensity for allele A and B at SNP j.

Preprocessing: allelic ratio normalization

These "raw" signals have to be preprocessed (or "normalized") in order to make them comparable (i) across samples for each locus, and (ii) across loci for each sample. While (i) is a common issue in all genomic analyses, we focus here on (ii), which is particularly crucial for copy-number data, where the locus of a given sample have to be segmented into genomic regions.

The TCGA project of molecular characterization of cancers has performed thousands of DNA copy number experiments for which both tumor sample and a normal sample (e.g. blood cells) from the same individual are available. For a given individual, we are therefore observing at each locus j (θ A ij , θ B ij ), where i ∈ {N, T } for the normal and the tumor sample, respectively. Let us define θ ij = θ A ij + θ B ij , and β ij = θ B ij /θ ij . The copy number data for one such tumor/normal pair are generally summarized using • total copy numbers in the tumor sample: c j = 2θ T j /θ N j ;

• allelic ratios (a.k.a. allele B fraction) in the normal sample: b N j = θ B N j /θ N j ;

• allelic ratios in the tumor sample: b T j = θ B T j /θ T j . Note that in the above definition of c j , the total intensity θ T j at locus j in the tumor sample is scaled by the corresponding quantity in the matching normal sample. The goal of this normalization is to correct for systematic and locus-specific biases in these intensities. These signals are displayed in the first three rows of Figure 9. the tumor allelic ratios after processing by the TumorBoost normalization method that we developed [J21, s3] in the context of the Aroma Project. As illustrated by the comparison between the (c) and (d) panels and between the (g) and (h) panels, this method provides a substantial denoising of the input data, with changes in the allelic ratio distribution being much more visible after normalization. The method takes advantage of two observations. First, the variability of allelic ratios at a given position is systematic, i.e. highly reproducible between the tumor and the paired normal sample. Second, one expects the true allelic ratios to be either 0, 1/2 or 1 in the normal samples. Therefore, the difference between the normal allelic ratio and the expected true one at a given locus provides a good estimate of the noise in the tumor allelic ratio, which can then be subtracted to estimate the tumor allelic ratio.

We have extended this method to the case where no paired normal is available in a collaboration with researchers at the University of Navarra [J18, s2]. We also stress the importance of the choice of a reference to estimated total copy numbers in this unpaired case, as illustrated in [B1]. While the methods described in this section are very simple from a mathematical point of view, their impact in the quality of the data is of major importance for downstream statistical analyses and results interpretation. Importantly, Figure 9.3 also illustrates the fact that change points occur at the same position in both signal dimensions. This is the case because changes in one of the parental copy numbers induce changes in both c and b (or d). In order to maximize the power of a change point detection method, it makes sense to require that a method can be applied jointly (rather than independently) to the two dimensions of the signals.

The problem of inferring the location of DNA copy number changes (also called copynumber change points) from locus-level estimates is an instance of the widely-studied changepoint detection problem. Two main statistical models have been used to model DNA copy number changes: change-point models and Hidden Markov Models (HMM). Before focusing on change point models, we give a brief overview of HMM. HMM assume that observed copy numbers at the locus level are emitted by an underlying Markov chain according to a small number of hidden true copy number states. HMM naturally incorporate and take advantage of the fact that different segments can have the same true copy number. Several HMM-based methods have been proposed for segmenting total copy numbers [START_REF] Fridlyand | Hidden Markov models approach to the analysis of array CGH data[END_REF][START_REF] Guha | Bayesian hidden Markov modeling of array CGH data[END_REF][START_REF] Lai | Stochastic segmentation models for array-based comparative genomic hybridization data analysis[END_REF]. These methods mainly differ in the assumptions that are made for the dynamics of the underlying Markov chain, and the approaches used for the estimation of the hidden states. Extensions of HMM to two-dimensional signals have also been proposed [START_REF] Chen | Estimation of parent specific DNA copy number in tumors using high-density genotyping arrays[END_REF].

We refer to [START_REF] Truong | A review of change point detection methods[END_REF] for a more comprehensive review on the general problem of change-point detection, which also covers the known statistical guarantees for the different methods. In the rest of this section we review change point models in the light of their application to copy number data. This section is adapted and updated from [B1]. In particular, we are interested in two key aspects:

• trade-offs between statistical accuracy and computational efficiency: sub-quadratic complexity in time or space is necessary for an algorithm to be applicable to copy number data;

• applicability to the joint segmentation of two-dimensional piecewise constant signals (c, d) which was not addressed in [B1].

A change-point model

A simple model for the observed DNA copy numbers c is:

c j = γ j + ε j , 1 ≤ j ≤ p (9.1)
where the errors (ε j ) 1≤j≤p are iid and the vector γ of true copy numbers is are assumed to be piecewise constant. For simplicity of exposition, (9.1) is formulated in terms of total copy numbers, but extending this model to two dimensions ((c, d) or c, c) is straightforward. The main difficulty with the two-dimensional version is of computational nature, as explained below.

We consider here the problem of minimizing the 2 loss, both for simplicity and for its relevance to DNA copy number segmentation. More general loss functions can be considered, see e.g. [START_REF] Rigaill | A pruned dynamic programming algorithm to recover the best segmentations with 1 to K max change-points[END_REF] which addresses the problem of minimizing the cost of a segmentation under the assumption that the cost is segment-additive. Under the 2 loss, our goal is to solve the following optimization problem:

min γ p j=1 (c j -γ j ) 2 s.t. Dγ 0 ≤ K (9.2)
where Dγ = (γ j+1 -γ j ) 1≤j≤p is the vector of first order differences associated to the true copy numbers γ, and • 0 is the 0 norm, that is, the number of non-null entries. The constraint Dγ 0 ≤ K encodes the fact that γ is piecewise constant, with at most K change points. Note that if we assume that the errors are Gaussian and homoscedastic, i.e. ε j ∼ N (0, σ 2 ), then maximizing the likelihood of the model (9.1) for a given number of change points K is equivalent to solving the optimization problem (9.2). Also, under this assumption, given change point locations, the optimal value of the true copy number between two change points is the average observed copy number in this region. Therefore, for a given K, solving (9.2) reduces to the combinatorial problem of finding the best possible change point locations according to (9.2). In practice, we are thus faced with two problems:

Combinatorial problem: (9.2) is non-convex, and the number of admissible change point locations is K p-1 , that is, O(p K ), so an exhaustive search of the best partition is practically infeasible in the case where p is large. This problem is discussed in Section 9.3.3;

Model selection problem: K is usually not known. As the segmentation models are nested with respect to K, the objective function in (9.2) is non-increasing in K. This model fit term has to be compensated by a penalty term that is increasing with the model size K:

min γ p j=1 (c j -γ j ) 2 + pen(K) . (9.3)
A variety of penalties have been proposed for the change point problem, see e.g. [START_REF] Truong | A review of change point detection methods[END_REF] for a general review and [START_REF] Hocking | Learning sparse penalties for changepoint detection using max margin interval regression[END_REF] for a focus on model section for DNA copy number change point detection. Formally, when the penalty is linear in K, 9.3 can be seen as the (Lagrangian) dual problem of (9.2). This is often the case, e.g. for the classical Akaike Information Criterion (AIC, [START_REF] Akaike | Information Theory and an Extension of the Maximum Likelihood Principle[END_REF]) or Bayesian Information Criterion (BIC, [START_REF] Schwarz | Estimating the dimension of a model[END_REF]).

Estimation of change point locations

The main difficulty to address the computational problem of estimating the change point locations is to find an appropriate balance between computational complexity, and statistical accuracy. We are interested in applications where K is of the order of 10 to 100, and p of the order of 10 5 to 10 6 . In such situations, a quadratic time or space complexity is already prohibitive. We review three types of approaches that we find both relevant to the biological problem, and statistically sound, and then discuss how these approaches may be combined.

Exact solutions. By taking advantage of the additivity of the objective function in the segments, it is possible to use a dynamic programming strategy to reduce the complexity of an exhaustive search from O(p K ) to O(K•p 2 ) [START_REF] Picard | A statistical approach for array-CGH data analysis[END_REF], at the price of increasing the space complexity from O(1) to O(p 2 ). This algorithm recovers the entire sequence of optimal segmentations in k segments for all k ≤ K, and it can be extended two multi-dimensional signals. However, a quadratic time and space complexity may be too high for recent microarray or sequencing copy number data. A pruned dynamic programming algorithm with linear space complexity has been proposed, that recovers the set of solutions faster [START_REF] Rigaill | A pruned dynamic programming algorithm to recover the best segmentations with 1 to K max change-points[END_REF]. In particular, although the worst case time complexity of the pDPA algorithm is still O(K•p 2 ), in practical situations it is almost linear in p. A related algorithm called PELT [START_REF] Killick | Optimal detection of changepoints with a linear computational cost[END_REF] provides the solution to the dual problem (9.3) with linear penalty pen(K) = λK in a linear time for a fixed value of λ. Although it does not provide the entire path of solutions (like pDPA does), it can be extended to multi-dimensional signals (unlike pDPA).

Binary segmentation. A widely used approach in various types of applications is binary segmentation [START_REF] Sen | On tests for detecting change in mean[END_REF], which recursively looks for the best partition of the data into two segments. This greedy algorithm has a very low (linear) time complexity of O(p log K). It has been adapted to the problem of copy-number segmentation by looking for the best partition in three segments at each step, with the constraint that the two extreme segments have the same true copy number [START_REF] Olshen | Circular binary segmentation for the analysis of array-based DNA copy number data[END_REF]. This method is called Circular Binary Segmentation (CBS). The depth of the recursion is determined by the estimated significance of the change points, which implicitly determines K. The original algorithm is quadratic in time due to the fast that two segment boundaries are looked for at each step. However, a pruned version which is also almost linear in time has been proposed [START_REF] Venkatraman | A faster circular binary segmentation algorithm for the analysis of array CGH data[END_REF]. We have proposed an extension of this method to the problem of segmenting both c and d [J19, s1]. It consists in a conditional segmentation method, where (i) each chromosome is segmented by CBS using the c statistic, and (ii) each segment is segmented by CBS using the d statistic. This extension inherits the quasi-linear complexity or CBS.

Convex relaxations. A classical approach in statistical machine learning is to replace a non-convex optimization problem with an approximate, but convex version of the problem, which can therefore be solved efficiently. Problem (9.2) can be convexified by replacing the 0 constraint on the number of change points by a 1 constraint. An adaptation of the fused lasso [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF] has been proposed [START_REF] Tibshirani | Spatial smoothing and hot spot detection for CGH data using the fused lasso[END_REF], which solves the constrained optimization problem The optimization problem (9.4) incorporates both a total variation constraint, that is, a constraint on the 1 norm of the jumps in γ, and sparsity constraint on γ -2, enforcing that most loci correspond to the normal copy number state. The complexity of the algorithm proposed in [START_REF] Tibshirani | Spatial smoothing and hot spot detection for CGH data using the fused lasso[END_REF] is (at best) O(p2 ), which as already discussed can be prohibitive for recent data sets. The second method [START_REF] Harchaoui | Multiple change-point estimation with a total variation penalty[END_REF] only keeps the total variation constraint, i.e. it simply convexifies Dγ 0 in (9.2) into Dγ 1 . This optimization problem can be written as a Lasso-type regression problem via the change of variable δ = Dγ. Therefore be solved in O(p•K 2 ) using a Least Angle Regression (LARS) algorithm [START_REF] Efron | Least angle regression[END_REF] to select the first K change points. An extension of this method to multi-dimensional signals into group-fused LARS has been proposed by [START_REF] Vert | Fast detection of multiple change-points shared by many signals using group LARS[END_REF][START_REF] Bleakley | The group fused Lasso for multiple change-point detection[END_REF], which inherits the linear complexity of the original method.

Combining approximate and exact methods. We also mention an interesting twostep strategy that consists in 1. performing a fast (yet approximate) search of K candidate change points, e.g. using binary segmentation or total variation penalty;

2. running dynamic programming only on this reduced set of candidates.

This two-step strategy can improve the exploration of all possible segmentations, while maintaining a low computational complexity. Indeed, when K p, the computational cost of the second step (O(K 3 )) can be negligible compared to the initial search. Such strategies have been proposed by [START_REF] Gey | Using CART to Detect Multiple Change Points in the Mean for Large Sample[END_REF] for binary segmentation, and by [START_REF] Harchaoui | Multiple change-point estimation with a total variation penalty[END_REF] for the total variation penalty. Importantly, this two-step approach can be used for the joint segmentation of multi-dimensional signals. These methods are implemented in the R jointseg [S2] which is available from CRAN 2 : in particular, the group-fused LARS method of [START_REF] Vert | Fast detection of multiple change-points shared by many signals using group LARS[END_REF][START_REF] Bleakley | The group fused Lasso for multiple change-point detection[END_REF] has been ported from Matlab to R, the CART/binary segmentation method of [START_REF] Gey | Using CART to Detect Multiple Change Points in the Mean for Large Sample[END_REF] has been implemented. The fast univariate implementation of pDPA [START_REF] Rigaill | A pruned dynamic programming algorithm to recover the best segmentations with 1 to K max change-points[END_REF] is also available in jointseg.

Performance evaluation

It is often the case that in absence of rich enough gold standard data sets, the numerical performance of statistical methods is assessed using unrealistic simulation studies, or based on small real data analyses with limited ground truth available. We have designed and implemented a framework to generate realistic DNA copy number profiles of cancer samples with known truth [J11, S4], that aims at combining the advantages and bypass some of the limitations of previous approaches [START_REF] Willenbrock | A comparison study: applying segmentation to array-CGH data for downstream analyses[END_REF][START_REF] Mosén-Ansorena | Comparison of methods to detect copy number alterations in cancer using simulated and real genotyping data[END_REF][START_REF] Lai | Comparative analysis of algorithms for identifying amplifications and deletions in array CGH data[END_REF][START_REF] Hocking | Learning smoothing models of copy number profiles using breakpoint annotations[END_REF]. In this section we describe this framework, and demonstrate how it can be applied to compare copy number segmentation methods. [START_REF] Rasmussen | Allele-specific copy number analysis of tumor samples with aneuploidy and tumor heterogeneity[END_REF] available from the Gene Expression Omnibus (GEO, [START_REF] Edgar | Gene Expression Omnibus: NCBI gene expression and hybridization array data repository[END_REF]). In this series, the cell line H1395 is mixed with a matched blood sample with several mixture proportions: 0, 30, 50 70 and 100% tumor cells. Only the last two mixture proportions are shown in Figure 9.4. The regions highlighted in purple were manually labelled as corresponding to two different copy number states (gain and loss of one DNA copy). Importantly, the labelling was performed on the 100% tumor sample where the signal to noise ratio is higher, but it is valid as well for the 70% tumor sample (and for the lower ones), where the signal to noise ratio is weaker. The regions highlighted in purple are then extracted and serve as a basis for a resampling-based data generation framework. These annotated data sets are available in the R package acnr [S4] which is available from CRAN3 .

Genomic regions with known copy-number state

Generating copy number profiles with known truth

A synthetic copy-number profile of length n with K change points can be generated in two steps from the annotated data described above: generation of truth: K change point positions (drawn uniformly out of the n -1 possible intervals between two successive loci), and K + 1 copy-number state labels for all K +1 regions between two consecutive change points, drawn from those of the existing annotated regions;

generation of signal: a n × 3 matrix of total copy numbers (c), allelic ratios (b), and genotypes. For each region of size n R between two change points, we sample n R data points from real copy-number data (such at the one displayed in Figure 9.1) corresponding to this type of region.

Using this framework, a variety of synthetic copy-number profiles can be generated. Figure 9.5 displays examples of four such profiles, which were generated from the same "truth" (change point position and region labels) as in Figure 9.3 with different tumor cell fractions (in rows), and from different annotated data sets (in column; the first column uses the data set described above, and displayed in Figure 9.4). Two attractive features of this Figure 9.5: Synthetic copy-number profiles that generated from the same "truth" as in Figure 9.3. Each block of two plots corresponds to total copy numbers (c) and allelic ratios (b) for one particular combination of fraction of tumor cells (in rows) and data set (in columns). Red vertical lines represent change points. Heterozygous SNPs are colored in black, and homozygous SNPs in gray. framework are the following:

• the distribution of copy number data in a given region matches a distribution actually observed in real data; therefore it does not rely on any probabilistic model assumption;

• the signal to noise ratio of these profiles can be tuned via the fraction of tumor cells, which has a clear biological interpretation.

Application to the joint segmentation of copy number profiles

In Section 9.3 we reviewed methods to segment copy number profiles with a focus on three main criteria: whether an algorithm solves the optimization problem (9.2) exactly or approximately, whether its time and space complexity is linear, quadratic, or worse, and whether it can be extended to two-dimensional signals (such as the ones displayed in Figure 9.3).

To complement this theoretical review, we empirically compare the performance of some of these algorithms on synthetic data generated as described in the preceding subsection. Several measures can be defined to assess the performance of segmentation methods, see e.g. [START_REF] Truong | A review of change point detection methods[END_REF]Section 3]. Here we focus on how well true change points are recovered, by casting the problem of change point detection as a binary classification problem. We define the number of true positives (TP) as the number of true change points for which at least one change point is detected closer than a given tolerance parameter. The number of false positives FP is FP = P-TP, where P is the total number of detected change points. With these definitions we are able for each method to plot a ROC curve parametrized by the number of change point found by the method. data points (n = 5, 000) was chosen in order for the dynamic programming (DP) approach to be applicable on the two-dimensional data (c, d), where it is quadratic as explained in Section 9.3.3. The methods compared are dynamic programming (DP, [START_REF] Picard | A statistical approach for array-CGH data analysis[END_REF][START_REF] Rigaill | A pruned dynamic programming algorithm to recover the best segmentations with 1 to K max change-points[END_REF]), (recursive) binary segmentation followed by dynamic programming (RBS + DP, [START_REF] Gey | Using CART to Detect Multiple Change Points in the Mean for Large Sample[END_REF]), circular binary segmentation (CBS [START_REF] Venkatraman | A faster circular binary segmentation algorithm for the analysis of array CGH data[END_REF] and PSCBS [J19]), and total variation penalty (LARS [START_REF] Harchaoui | Multiple change-point estimation with a total variation penalty[END_REF] and group-fused LARS [START_REF] Bleakley | The group fused Lasso for multiple change-point detection[END_REF]). For all these methods, we compared the 1d-segmentation of total copy numbers (c) alone to the 2d segmentation of (c, d). The main observations that may be drawn from these experiments are the following:

• The true positive rate seems to be bounded away from 1, meaning that on average, only 4 of the 5 segments could be identified.

• Most 2d methods (DP, RBS, CBS) generally perform better than their 1d counterpart, although 2d-DP is outperformed by 1d-DP at a high specificity (ie for a low number of false positives);

• The fact that the 2d total variation (GFLars) performs worse than is 1d counterpart (Lars) is an artifact of the implementation: in the implementation the data points for which one of the dimensions of the signal is missing are not considered, resulting in a substantial loss of power. Indeed, in this application only SNPs that were heterozygous in the germline have non-missing signals for d, and these SNP typically represent less than one third of the total number of SNPs.

• Although DP is the only method that exactly minimizes the objective function, it is not necessary the best performer. In the example of Figure 9.6 the 2d-RBS+DP method clearly outperforms DP at high specificities.

The above results illustrate the type of evaluation that can be performed using this framework. A general conclusion of this performance assessment is that no single method performs always better than the other ones. In particular, some methods are more robust than others to a decreasing tumor fraction, and the performance of the methods may also differ across types of change points (i.e. between a normal region and a region of gain).

An important point from a statistician's perspective is that most of the above observations would most probably not have been made if the performance of the methods had been assessed by classical simulations based on a probabilistic model for the noise. In particular, one would expect DP to have better statistical performance than its competitors, because it solves the original optimization problem (9.2) and not a relaxation. However, the statistical justification of considering this optimization problem relies on an assumption of homoskedastic Gaussian noise, which may not be appropriate for these data.

Chapter 10

Adjacency-constrained clustering

In the context of Genome Wide Association Studies (GWAS), region-scale approaches taking haplotype blocks into account can result in substantial statistical gains [START_REF] Gabriel | The structure of haplotype blocks in the human genome[END_REF]. Hi-C studies [START_REF] Dixon | Topological domains in mammalian genomes identified by analysis of chromatin interactions[END_REF] have demonstrated the existence of topological domains, which are megabase-sized local chromatin interaction domains correlating with regions of the genome that constrain the spread of heterochromatin. Both contexts raise the natural question of segmenting similarity matrices, either in the form of linkage disequilibrium (LD) matrices or of Hi-C contact maps.

In this chapter we study a particular segmentation method, hierarchical agglomerative clustering (HAC) with adjacency constraints. After defining the method and studying its conditions of applicability, we describe a quasi-linear algorithm to perform adjacencyconstrained HAC under an additional assumption derived from the biological context. Finally we describe how this method can be applied to GWAS studies, and combined with group lasso regression to detect LD blocks associated to a given phenotype.

HAC and adjacency-constrained HAC

HAC was initially described by Ward [START_REF] Ward | Hierarchical grouping to optimize an objective function[END_REF] for data in R d . Let Ω := {x 1 , • • • , x p } be a set of p objects to be clustered, with x i ∈ R d for i = 1, . . . , d. A cluster is a subset of Ω. The degree of inhomogeneity of a cluster G ⊂ Ω is quantified by the inertia (also known as Error Sum of Squares, ESS): (10.

I(G) = i∈G x i -xG 2 , ( 10
2)

The quantity δ is known as Ward's linkage and it is equal to the variation of within-cluster inertia (also called within-cluster sum of squares) after merging two clusters. It also corresponds to a scaled version of the squared distance between centers of gravity:

δ(G, G ) = |G||G | |G| + |G | xG -xG 2 . (10.3)
The HAC algorithm is described in Algorithm 2. Starting from the trivial partition P 1 = {{x 1 }, {x 2 }, • • • , {x p }} with p singletons, the HAC algorithm creates a sequence of partitions by successively merging the two clusters whose linkage δ is the smallest, until all objects have been merged into a single cluster.

Algorithm 2 (Contiguity-constrained) Hierarchical Agglomerative Clustering (HAC)

1: Initialization: P 1 = {{x 1 }, {x 2 }, • • • , {x p }} 2: for t = 1 to p -1 do 3:
Compute all pairwise linkage values between (contiguous) clusters of the current partition P t

4:

Merge the two (contiguous) clusters with minimal linkage value to obtain the next partition P t+1 5: end for 6: return {P 1 , P 2 , . . . , P p } Contiguity-constrained HAC is a simple modification of Algorithm 2, where only clusters that are adjacent are candidate mergers, as indicated by the word "contiguous" in line 3 and 4. The idea of incorporating such constraints was previously mentioned by Lebart [START_REF] Lebart | Programme d'agrégation avec contraintes[END_REF] to incorporate geographical (two-dimensional) constraints to cluster socio-economic data, and by [START_REF] Michel | A supervised clustering approach for fMRI-based inference of brain states[END_REF] to cluster functional Magnetic Resonance Imaging (fMRI) data into contiguous (three-dimensional) brain regions.

While contiguity-constrained HAC can be defined for any symmetric relation indicating which pairs of objects are considered as "contiguous", we focus in this chapter on the specific case of adjacency-constrained clustering. That is, we assume that the objects to cluster are ordered along a line, and only allow for adjacent objects or clusters of objects to be merged. This particular case has been studied by Grimm [START_REF] Grimm | CONISS: a fortran 77 program for stratigraphically constrained analysis by the method of incremental sum of squares[END_REF], and an R package implementing this algorithm, rioja [START_REF] Juggins | rioja: Analysis of Quaternary Science Data[END_REF], has been developed1 .

Formulation using pairwise Euclidean distances. The inertia of a cluster may be expressed only in function of sums of the entries of the pairwise distances ( x i -x j , 1 ≤ i, j ≤ p):

I(G) = ∆(G, G) 2|G| , (10.4)
where ∆ is defined by ∆(G, G ) = xi∈G,x i ∈G x i -x j 2 for any clusters G and G . Moreover, Ward's linkage between any two clusters G and G may be itself be written in function of these pairwise distances:

δ(G, G ) = 1 2 ∆(G ∪ G , G ∪ G ) |G ∪ G | - ∆(G, G) |G| - ∆(G , G ) |G | . ( 10
.5)

Extensions to possibly non-Euclidean settings

The HAC algorithm of Ward [START_REF] Ward | Hierarchical grouping to optimize an objective function[END_REF] has been designed to cluster elements of R d . In this chapter, we focus on the situation where the objects to be clustered are indirectly described by a matrix of pairwise similarities S = (s ij ) i,j=1,...,p . The motivation of this section is an important difference between GWAS and Hi-C from a statistical perspective. While the LD (r 2 ) similarity is a kernel (as explained below in Section 10.4), it is not necessarily the case for Hi-C similarity matrices, possibly leading to reversals in dendrograms obtained from constrained HAC as studied in detail in [J1].

Kernels. If S is positive definite, the theory of Reproducing Kernel Hilbert Spaces [START_REF] Aronszajn | Theory of reproducing kernels[END_REF] implies that the data can be embedded in an implicit Hilbert space. This allows to define the inertia of a cluster, and consequently to formulate Ward's linkage between any two clusters in terms of the similarity using the so-called "kernel trick":

∀ C, C ⊂ {1, . . . , p}, δ(C, C ) = S(C) |C| + S(C ) |C | - S(C ∪ C ) |C ∪ C | , (10.6)
where S(C) = (i,j)∈C 2 s ij only depends on S and not on the embedding. This expression shows that Ward's Linkage also has a natural interpretation as the decrease in average intracluster similarity after merging two clusters. Equation (10.6) is derived from (10.3) using the formalism of kernels in Section S1.1 of the Supplementary material of [J5]. Another way to prove it is to use Equation (10.5) with the distance associated to S by the kernel mapping, d 2 ij = s ii + s jj -2s ij , and noticing that the diagonal terms cancel out.

Similarities. An extension of this approach to the case of a general (that is, possibly nonpositive definite) similarity matrix has been studied by Miyamoto, Abe, Endo, and Takeshita [START_REF] Miyamoto | Ward method of hierarchical clustering for non-Euclidean similarity measures[END_REF]. Noting that (i) for a large enough λ, the matrix S λ = S + λI p is positive definite and that (ii) Ward's linkage associated to S λ is simply given by δ λ (C, C ) = δ(C, C )+λ, applying Ward's HAC to S and S λ yields the exact same hierarchy, only shifting the linkage values by +λ [START_REF] Miyamoto | Ward method of hierarchical clustering for non-Euclidean similarity measures[END_REF]Theorem 1]. This result, which a fortiori holds for adjacency-constrained Ward's HAC, justifies the use of Equation (10.6) in the case of a general similarity matrix.

Dissimilarities. For completeness, we also describe a generic construction of HAC for arbitrary dissimilarity data, proposed by [START_REF] Chavent | ClustGeo2: an R package for hierarchical clustering with spatial constraints[END_REF]. This construction is based on an analogy between distance and dissimilarity. If the p objects are described by a dissimilarity matrix D, that is, a matrix satisfying: for all i, j ∈ {1, . . . , p}, d ij ≥ 0; d ii = 0; d ij = d ji , then one can simply define the inertia of a cluster by (10.4), with (sums of squared) distances replaced by (sums of squared) dissimilarities. The corresponding HAC is then formally obtained as the output of Algorithm 2. This construction generalizes the (less intrinsic) one proposed by [START_REF] Székely | Hierarchical clustering via joint between-within distances: extending Ward's minimum variance method[END_REF] and, later, [START_REF] Strauss | Generalising Ward's method for use with Manhattan distances[END_REF], for specific dissimilarities: in those papers, instead of starting from Equation (10.4), the linkage between two clusters is defined by a reformulation of (10.5) using the fact that ∆ inserting an element in an ordered list is linear in p, again leading to a quadratic time complexity. Instead, we propose to store the candidate fusions in a partially-ordered data structure called a min heap [START_REF] Williams | Algorithm 232 -heapsort[END_REF]. This type of structure achieves an appropriate trade-off between the cost of maintaining the structure and the cost of finding the minimum element at each iteration. A min heap is a binary tree such that the value of each node is smaller than the value of its two children. The advantage of this structure is that all the operations required in Algorithm 2 to create and maintain the list of candidate fusions can be done in O(log p). A detailed description of the algorithm, which is implemented in the adjclust package [S3], is given in [J5].

(G ∪ G , G ∪ G ) = ∆(G, G) + ∆(G , G ) + 2∆(G, G ). P (j, k) P (i, k) k h S(C) i i j j P (j, h) P (i, h) k h S (

Application to Genome-Wide Association Studies

As explained in Chapter 1, genome-Wide Association Studies (GWAS) aim at identifying markers associated with a phenotype of interest. Given that an individual's genotype is characterized by millions of markers (known as SNPs) this approach yields a large multiple testing problem. Due to recombination phenomena, the hypotheses corresponding to SNPs that are close to each other along the genome are statistically dependent. This dependence is usually quantified by the linkage disequilibrium (LD), as displayed in the right panel of Figure 1.2. A widely used measure of LD in the context of GWAS is the r2 coefficient, which can be estimated directly from genotypes measured by genotyping array or sequencing data using standard methods [START_REF] Clayton | snpStats: SnpMatrix and XSnpMatrix classes and methods[END_REF]. The similarity S = (r 2 ij ) i,j induced by LD can be shown to be a kernel [J5, Supplementary Materials] 2 .

In this section, we show how the above-described algorithm for adjacency-constrained clustering may be used in this context. In Subsection 10.4.1 we illustrate the relevance of the band similarity assumption introduced above, and the computational gains offered by the algorithm. In Subsection 10.4.2 we show how this algorithm can be used to perform tests of association at the level of LD blocks instead of SNPs. The numerical experiments below have been performed on a SNP dataset coming from a GWA study on HIV [START_REF] Dalmasso | Distinct genetic loci control plasma HIV-RNA and cellular HIV-DNA levels in HIV-1 infection: the ANRS Genome Wide Association 01 study[END_REF] mentioned in Chapter 1.

Linkage disequilibrium block inference in GWAS

We used genotype data corresponding to five chromosomes that span the typical number of SNPs per chromosome observed on the 317k Illumina genotyping microarray (5, 000 ≤ p ≤ 25, 000).

Quality of the band approximation. We compared the dendrogram obtained with h < p to the reference dendrogram obtained with the full bandwidth (h = p) by recording the index t of the last clustering step (among p -1) for which all the preceding fusions in the two dendrograms are identical. The quantity t/(p -1) can then be interpreted as a measure of similarity between dendrograms, ranging from 0 (the first fusions are different) to 1 (the dendrograms are identical). at h = 1000, the dendrograms differ from the reference dendrogram only in the last 0.5% of the clustering step. We obtained similar results when using other criteria for evaluating the quality of the band approximation, including Baker's Gamma correlation coefficient [START_REF] Baker | Stability of two hierarchical grouping techniques case I: sensitivity to data errors[END_REF], which corresponds to the Spearman correlation between the ranks of fusion between all pairs of objects. Importantly, the influence of the bandwidth parameter is the same across chromosomes, that is, across values of p. Therefore, it makes sense to assume that h does not depend on p and that the time and space complexity of our proposed algorithm, which depends on h, is indeed quasi-linear in p.

Scalability and computation times. Figure 10.2 displays the computation time for the LD matrix (dotted lines) and for the adjacency-constrained HAC with respect to the size of the chromosome (x axis), both for rioja (dashed line) and adjclust (solid lines). As expected, the computation time for rioja did not depend on the bandwidth h, so we only represented h = p. For adjclust, the results for varying bandwidths are represented by different colors.

First, the computation times of rioja are several orders of magnitude larger than those of adjclust, even when h = p where both methods implement the exact same algorithm. As expected, the complexity of adjclust with h = p is quadratic in p, while it is essentially linear in p for fixed values of h < p. For large values of p the gain of the band approximation is substantial. We also note that regardless of the value of h, the total time needed for the clustering is of the order of (and generally lower than) the time needed for the computation of the LD. This implies that runnning the adjclust algorithm in this case comes essentially for free from a computational point of view.

Block-wise GWAS

To overcome the intrinsic limitations of classical Genome-wide association studies (GWAS) involving univariate tests between statistically dependent markers, we have proposed in [J9] the following approach: (i) infer LD blocks using adjclust, (ii) estimate the number of LD blocks using a model selection criterion based on the Gap statistic [START_REF] Tibshirani | Estimating the number of clusters in a data set via the gap statistic[END_REF], and (iii) perform Group Lasso regression to identify which of the inferred LD blocks are associated with the response.

We have investigated numerically the efficiency of this approach compared to state-ofthe art regression methods: haplotype association tests, single marker analyses, and Lasso and Elastic-Net regressions. Our numerical experiments show that the proposed method outperforms state-of-the-art approaches (the haplotype association module of the PLINK genome association analysis tool [START_REF] Purcell | PLINK: a tool set for wholegenome association and population-based linkage analyses[END_REF], single marker analyses and Lasso and Elastic-Net regressions) as soon as the number of causal SNPs within a LD block exceeds 2. This method has also been applied to the analysis of the HIV data set introduced in Chapter 1, where one goal is to identify LD blocks associated to the viral load (that is, the abundance of the virus in blood cells). This is illustrated by Figure 10.3, where a small part of the LD similarity matrix is displayed together with the results of the proposed method (left panel) and PLINK (right panel). The region displayed corresponds to 68 contiguous markers located in the major histocompatibility complex (MHC), of which 8 (marked with a red star *) had previously been identified as (marginally) associated to the viral load [START_REF] Dalmasso | Distinct genetic loci control plasma HIV-RNA and cellular HIV-DNA levels in HIV-1 infection: the ANRS Genome Wide Association 01 study[END_REF]. The proposed method (Figure 10.3, left plot) is able to identify LD blocks (delimited with a red contour line) larger than PLINK. within the structure inferred by the proposed method, the blocks selected by the Group Lasso are delimited with a red contour line. The SNPs selected by single marker analysis are plotted with a red star (*), and the SNPs missed by Lasso with a blue dash (-). Right: within the structure inferred by the haplotype association method, the blocks selected by the competing method are delimited with a red contour line.

Emerging challenges in post hoc inference 11.1.1 Improved post hoc bounds

New JER controlling procedures. The generic framework based on JER control introduced in Section 5 makes it possible to construct several post hoc bounds based on different reference families. One natural extension is to look for new probabilistic inequalities that yield JER control, or to extend the validity of existing one. In particular, it is conjectured in the multiple testing field that the Simes inequality holds for Gaussian tests statistics with arbitrary covariance matrix. If this conjecture was proved, then the validity of the Simes post hoc bound would be automatically extended as well. Another interesting research direction is to obtain JER control for dependent tests by building reference families based on harmonic mean p-values [START_REF] Wilson | The harmonic mean p-value for combining dependent tests[END_REF].

Aggregation of reference families. Bounds based on distinct reference families may perform differently on different subsets of selected hypotheses. For example, Figure 5.6 in Chapter 5, suggests (consistently with our theoretical results) that the Linear bound performs better than the Beta bound for smaller subsets, and worse for larger subsets. Similar observations have been made in the case of locally-structured hypotheses [J3]. An interesting perspective is to use aggregation techniques to build bounds performing almost as well as the best possible bound for any size or shape of subset.

Toward Post Selection Inference and online multiple testing

Connection to post selection inference. The post hoc bounds originally proposed by [START_REF] Goeman | Multiple testing for exploratory research[END_REF] and our JER-based bounds [J2] intrinsically rely on testing marginal hypotheses. Valid inference methods for multiple regression after arbitrary model selection have been proposed by [START_REF] Berk | Valid post-selection inference[END_REF][START_REF] Bachoc | Uniformly valid confidence intervals post-model-selection[END_REF]. These methods involve inflating standard confidence intervals by a constant factor called "Post Selection Inference (PoSI) constant". However, this approach cannot be used in high-dimensional settings due to the exponential time complexity (in the number p of variables) required to calculate the PoSI constant for a given design matrix. Kuchibhotla, Brown, Buja, Cai, George, and Zhao [4] have recently introduced confidence intervals enjoying two remarkable properties: they have a much reduced computational cost of O(p 2 ), and their volume scales with the size of the selected model and not with the largest possible model size considered. An important perspective is to derive post hoc bounds from [4], and to compare these bounds to our "marginal" bounds.

Post hoc bounds for online tests. Online multiple testing consists in performing sequence of tests, when the number of tests is unknown an possibly infinite. Motivated by applications to A/B testing, much progress has been recently made for online FDR and FDP control, see [START_REF] Javanmard | Online rules for control of false discovery rate and false discovery exceedance[END_REF] and the works of Aaditya Ramdas on the subject, e.g. [START_REF] Katsevich | Simultaneous high-probability bounds on the false discovery proportion in structured, regression, and online settings[END_REF]10]). A natural question is the construction of online post hoc bounds.

Inference for structured signals

The problem of differential analysis along the genome has been studied for DNA methylation data [START_REF] Robinson | Statistical methods for detecting differentially methylated loci and regions[END_REF] and ChIP seq data [START_REF] Lun | De novo detection of differentially bound regions for ChIP-seq data using peaks and windows: controlling error rates correctly[END_REF] among other types of genomic data. Just like differential expression analyses, the goal is to identify genomic features whose "activity" is significantly different between two biological or clinical conditions, based on the observation of several activity profiles. However, one is generally more interested in detecting differentially active regions (differentially methylated regions, or differentially bound regions in the above two examples) than differentially active individual loci. The main statistical difficulty is that the candidate regions are not known a priori and thus have to be constructed from the data. We are interested in two specific practical instances of such structured differential analysis problems. The first one is the case of differential DNA copy number analysis in cancers between cancer subtypes. The second one is differential analysis of Hi-C contact maps. For the latter, at least two "regional" scales should be considered, which are known as Topologically Associated Domains (TAD) and A/B compartments. These scales correspond to two levels of chromosomal organization that can be detected from Hi-C contact maps from a single DNA sample [START_REF] Dixon | Topological domains in mammalian genomes identified by analysis of chromatin interactions[END_REF].

In both contexts, it is possible to first perform differential analysis at the most local level, using methods derived from differential expression studies proposed in [START_REF] Robinson | Statistical methods for detecting differentially methylated loci and regions[END_REF][START_REF] Lun | De novo detection of differentially bound regions for ChIP-seq data using peaks and windows: controlling error rates correctly[END_REF][START_REF] Lun | diffHic: a Bioconductor package to detect differential genomic interactions in Hi-C data[END_REF], and then to aggregate this information into a regional signal. Using the strategy outlined in [J3] it is possible to derive post hoc bounds for any candidate region. Depending on the constraints on the shape of the regions of interest, finding efficient ways of scanning all relevant regions raises interesting computational challenges.

From a statistical perspective, following the argument developed in Chapter 8 to motivate the method proposed in [J16], aggregating marginal tests can be under-powered compared to performing multivariate tests. An additional difficulty in the above applications is that contrary to the case of differential expression of pathways, the regions of interest are generally not known in advance. Accounting for this feature will rely on testing strategies tailored to each problem. For differential copy number analysis, one direction to investigate consists in fitting a Hidden Markov Model in order to obtain locus-level probabilities of differential expression that account for spatial dependence, and to derive a tailored JER control. Another relevant direction consist in first performing a joint segmentation of the genome as in Chapter 9, and then design tailored tests of differential copy number analysis of such data-driven segments. For differential Hi-C analysis, one possibility is to elaborate on the TAD detection method presented in Chapter 10. We expect that these methodological developments will also raise interesting statistical and computational challenges.

Broader scientific challenges

In this section I discuss broader scientific challenges that I believe as a statistician to be fundamental for the research community to work in a sound and efficient manner. While I will certainly not solve these challenges by myself, I hope to be able to make some contributions to some of the points they raise in future years.

Interdisciplinarity: methods to connect theory and practice

An important context element for method development in genomic data analysis is the rapid evolution of genomic assays, which has several consequences for statisticians. First, the increasing complexity of these assays makes it difficult for non-specialists to gain enough understanding to propose relevant models for statistical and experimental variability. Then, the number of proposed bioinformatic methods makes it difficult to identify state-of-theart methods, and the fact that these methods are typically pipelines made of several tools hampers the statistical understanding of the strengths and weaknesses of such methods. Finally, and as a result, the time needed to finely address the statistical issues raised by the analysis of a given technology can be much longer than the half-life of this technology. Thus, the breadth of the interface between biology/genomics applications and statistical theory is expanding. An important point for future years is the ability for the scientific community to strengthen the links between theory and applications.

My personal experience is that it is difficult to foster scientific interactions between biology and mathematics (or even biology and statistics) at a macroscopic scale. I have found it more effective and more stimulating to organize seminars, reading groups or workshops dedicated to the statistical, bioinformatic and biological challenges raised by the analysis of one particular type of data, or one particular biological question. In the context of the PhD thesis of Nathanaël Randriamihamison on Hi-C data analysis, we have created in 2018 in Toulouse a very active reading group called chrocogen (for chromatin Conformation and gene expression), which gathers statisticians, biostatisticians, bioinformaticians and biologists. We have also organized a successful workshop on Hi-C data analysis in December 2019.

Reproducible research and performance evaluation

Reproducible research and the role of statisticians. The ability to reproduce or replicate scientific experiments is at the very core of the scientific method. Almost fifteen years ago, a paper raised awareness by claiming that an important proportion of published biomedical studies could not be replicated [START_REF] Ioannidis | Why most published research findings are false[END_REF]. This issue is not specific to biomedical sciences1 . While the actual proportion of "wrong" papers in the scientific literature can be discussed (see e.g. [START_REF] Leek | Statistics: P values are just the tip of the iceberg[END_REF][START_REF] Jager | An estimate of the science-wise false discovery rate and application to the top medical literature[END_REF]), there is consensus that this proportion is unacceptably large. Most results in biomedical studies are obtained from some data analysis. A first and obvious necessary condition for reproducing such computational results is the presence of all the required data and code. A number of journals now favor methods reproducibility by imposing to publish code and data allowing to reproduce the results of a paper; see also the special collection "Challenges in irreproducible research" from the Nature journals. The question of "reproducibility" of research findings has several important aspects from a statistician's perspective:

1. Part of non-reproducibility can simply be attributed to the fact that a number of studies still disregard or poorly address the issue of selection or multiple testing.

Even when all of the hypotheses tested are explicitly reported in a study, it is not uncommon that only those whose p-value is less than 0.05 are retained and mentioned in the abstract. This issue can readily be addressed by "classical" multiple testing procedures, and it is worrisome that not all scientific journals seem to be able to catch these errors before publication.

2. Scientific results may be non-reproducible because of selection biases in the study, such as silently adding/removing some cases in a study, reporting only the results of one particular combination of method and tuning parameters whereas many such combinations have also been tried. This issue of fishing for significance can partly be addressed by pre-registration of the methods to be used: such practices are highly relevant to confirmatory analyses, clinical trials, but less so in exploratory contexts. We believe that developments in post hoc inference could help address this issue.

3. Similar issues can be found in the statistical methodology literature itself. In papers presenting new statistical methods, it is not uncommon that the performance of several methods are tested on several synthetic or experimental data sets and/or performance criteria, but only a subset including those favoring the method proposed by the authors of the paper are reported. This issue has been studied in depth in [START_REF] Boulesteix | Publication bias in methodological computational research[END_REF], and practical guidelines for benchmarking studies have been proposed in [START_REF] Weber | Essential guidelines for computational method benchmarking[END_REF].

Importance of performance evaluation in scientific literature. The bias in the methodological literature discussed in the last item partly comes from the current publication system itself. Indeed, in order for a statistical method to be published, most journals require evidence of the superiority of this method on existing ones in simulations or data set analyses, which encourages the above-discussed selective reporting of numerical experiments. This lack of objectivity is problematic per se (leading for example to very contrasted performance assessments for a given method in different papers), but also because the literature on statistical methods often provides only limited insight into when and why a given method is appropriate or not, whereas these questions are actually those that contribute to a global understanding and improvement of the available methods. While the need for neutral comparison studies (also known as benchmarking studies) is widely acknowledged by us as statisticians, it is not always (or even often not) applied in our own methods publications, as noted in [START_REF] Boulesteix | Towards evidence-based computational statistics: lessons from clinical research on the role and design of real-data benchmark studies[END_REF]. The question of designing such benchmarking studies in several fields has itself become a subject of research, see e.g. [START_REF] Boulesteix | On the necessity and design of studies comparing statistical methods[END_REF]8,[START_REF] Robinson | Benchmarking comes of age[END_REF][START_REF] Weber | Essential guidelines for computational method benchmarking[END_REF] for recent contributions and the STRATOS Initiative2 in clinical biostatistics. Some statistical journals including Briefings in Bioinformatics or Genome Biology explicitly welcome neutral comparison studies. These stimulating questions provide food for thought for building novel forms of scientific publications that could be specifically designed for these types of studies.

Toward truly applicable methods

Although the applied statistics community has made tremendous efforts to make newly developed methods available for the bioinformatics or biology community, I believe that there is still much room for improvement in this direction, and that such improvement is crucial to consolidate the links between theory and applications. R packages. The past decades have undergone massive developments of open-source implementations of data analysis methods, especially in R and python. The current standard for the implementation of statistical methods is in the form of R packages. From a statistician's perspective, there is a clear tension between the collective interest of putting together and maintaining well-documented packages, and the corresponding individual time requirements -in particular, the more a package is used, the more costly it is to maintain it. As a result, the sweet spot for developers of statistical methods is currently to put together a (github) repository where their code is available, generally in the form of an R package3 . However, an important caveat of this package-centered model is that several (and possibly dozens of) packages using different inputs and outputs can be devoted to similar statistical problems. This makes it difficult for end users to know which implementation to use for a particular task 4 .

The computational statistics community should continue to foster the development of open-source implementations of new statistical methods. The scikit-learn Python library for machine learning is based on a different model, relying on a growing user and developer community to contribute to the code and documentation of a single module, which is a well-recognized entry point. This project has recently been awarded a prestigious prize from the Académie des Sciences, but it has required a huge time investment from the main contributors 5 . The time saved (collectively) by having a readily usable implementation of recent statistical methods is invaluable, not only for potential users, but also for developers of competing methods. I believe that it is important for the efficiency and quality of current and future scientific research to increase the (individual) support and recognition of this type of contributions.

Interactive visualization and inference. Interactive graphical user interfaces (GUI) are an important tool for the diffusion of statistical methods toward end users with little or no programming skills. Post hoc inference methods play particularly nicely with interactive visualization tools because these methods by construction allow the user to perform data snooping without compromising the statistical guarantees offered by the methods, thereby making statistical analysis truly interactive. The development of such tools has been made easier for R users with the advent of Shiny Applications. For example, we have quickly developed a proof of concept application as a supplementary material for the book chapter [P1], in order to illustrate the interest of post hoc bounds for localized signals 6 . I plan to make future contributions available in the form of such applications, in particular for the interactive visualization and analysis of DNA copy number profiles in cancers, and for the differential analysis of Hi-C data. Finally, such developments are a not only a way to bridge the gap between theory and practice: developing such applications together with end users can also yield a better understanding of which type of signal is relevant, therefore possibly opening the way for new statistical developments to specifically detect this particular type of signal.

Contents 1 . 1 1 . 2 1 . 3 1 . 4 1 . 5

 1112131415 Foreword: scientific path . . . . . . . . . . . . . . . . . . . . . . 2 From hypothesis-driven to data-driven research . . . . . . . . 2 Instances of biomedical questions . . . . . . . . . . . . . . . . . 3 Statistical challenges and opportunities . . . . . . . . . . . . . 5 Overview of contributions . . . . . . . . . . . . . . . . . . . . . 7 1

Figure 1 . 1 :

 11 Figure 1.1: DNA copy numbers for the lung tumor cell line NCI-H1395-4W. Left: karyotype. Right: DNA copy number profile of Chromosome 5.

Figure 1 . 2 :

 12 Figure 1.2: Example of a Manhattan plot (left panel) and a linkage disequilibrium plot (right panel) obtained from a GWA study on HIV. Credits: [103, Figure 1 and 2].

  Let us denote by (p (1) ≤ p (2) ≤ . . . ≤ p (m) ) the ordered p-values. Definition 2.1 (Step-up and step-down multiple testing procedures). Let c = (c i ) i∈Nm be a non-decreasing set of values in [0, 1].

Figure 2 . 1 :

 21 Figure 2.1: FDR control is not FDP control.

Contents 3 . 1 FDP

 31 as a stochastic process of a random threshold . . . . . . 28 3.2 Results for PI-0 procedures . . . . . . . . . . . . . . . . . . . . 31 3.3 Results for PI-1 procedures . . . . . . . . . . . . . . . . . . . . 33 3.4 Extensions to other dependency settings . . . . . . . . . . . . 34

. 2 )

 2

2 ;

 2

Figure 3 . 1 :

 31 Figure 3.1: Dual interpretations of the BH threshold. Left: ordered p-values; right: empirical cumulative distribution function of the p-values.

( 3 . 2 )Figure 3 . 2 :

 3232 Figure 3.2: Critical value of the BH procedure: illustration of the results of[START_REF] Chi | On the performance of FDR control: constraints and a partial solution[END_REF] in the Gaussian (left; no criticality, α = 0) and Laplace/double exponential (right: criticality; α > 0) location models.

4. 1

 1 Settings 4.1.1 Random effects model We consider the random effects model described in Section 3.1, with additional assumptions on the distribution of the test statistics. In particular, we study the specific case of Subbotin location and scale models where the density d of the test statistics conditionally on θ i = 0 is the ζ-Subbotin density defined in Section 2.1.1: d(x) = (L ζ ) -1 e -|x| ζ /ζ , for ζ ≥ 1, where the normalization constant L ζ is defined in Equation (2.3). The density of the test statistics conditionally on θ i = 1 is either a shift or a scaling of d. In the Subbotin location model, we have d 1,m (x) = d(xµ m ), for some (unknown) location parameter µ m > 0. The corresponding p-values are obtained using the transformation p i = D(X i ), and the cumulative distribution function of the p-values conditionally on

Theorem 4 . 1

 41 ([J17], Corollary 4.4). Take ζ > 1, γ = 1 -ζ -1 for the location case and ζ ≥ 1, γ = 1 for the scale case. Consider a ζ-Subbotin density (2.3) in the sparsity regime τ m = m β , 0 < β ≤ 1 and assume that C m is bounded away from 0 and 1. Letting r m = α m + (log(α -1 m /(log m) γ )) + /(log m) γ , the following holds: (i) The pFDR threshold t m (α m ) is asymptotically optimal if and only if α m → 0 and log α m = o ((log m) γ ) , (4.3) in which case it is asymptotically optimal at rate r m . (ii) The FDR threshold tF DR m at a level α m satisfying (4.3) is asymptotically optimal at rate r m . (iii) Choosing α m ∝ 1/(log m) γ , pFDR and FDR thresholding are both asymptotically optimal at rate r m = 1/(log m) γ .

Figure 4 .

 4 Figure 4.1 compares the relative excess risks of the Bayes procedure and the BH procedure for several choices of α m in the Gaussian location model. For each procedure (in rows) and each value of m (in columns), the behavior of the relative excess risk is studied as the (unknown) true model parameters (β, C m ) vary in [0, 1] × [0, 1]. The threshold of the Bayes procedure depends on the model parameters (β, C m ); the results reported in the first row "Bayes0" correspond to the choice (β, C m ) = (β 0 , C 0 ). For FDR thresholding, we report the results for the choice α m = α opt m (β 0 , C 0 ). The values for β 0 and C 0 are arbitrarily chosen as the midpoints of the corresponding intervals, i.e. (β 0 , C 0 ) = (1/2, 1/2).Bayes0 performs well when the sparsity parameter β is correctly specified, and its performance is fairly robust to C m . However, it performs poorly when β is misspecified, and increasingly so as m increases. The results is markedly different for FDR: FDR thresholding is less adaptive to C m than Bayes0, but much more adaptive to the sparsity parameter β, as illustrated by the fact that the configurations with low relative excess risk span the whole range of β. For FDR thresholding at α m = α opt m (β 0 , C 0 ), the fraction of configurations (β, C m ) for which E m ≤ 0.1 increases as m increases. This illustrates the asymptotic optimality of FDR thresholding for classification under sparsity in the Gaussian location model. Similar conclusions hold for the Gaussian scale and the Laplace scale models, see the Supplementary Materials of[J17]).

  Figure 4.1 compares the relative excess risks of the Bayes procedure and the BH procedure for several choices of α m in the Gaussian location model. For each procedure (in rows) and each value of m (in columns), the behavior of the relative excess risk is studied as the (unknown) true model parameters (β, C m ) vary in [0, 1] × [0, 1]. The threshold of the Bayes procedure depends on the model parameters (β, C m ); the results reported in the first row "Bayes0" correspond to the choice (β, C m ) = (β 0 , C 0 ). For FDR thresholding, we report the results for the choice α m = α opt m (β 0 , C 0 ). The values for β 0 and C 0 are arbitrarily chosen as the midpoints of the corresponding intervals, i.e. (β 0 , C 0 ) = (1/2, 1/2).Bayes0 performs well when the sparsity parameter β is correctly specified, and its performance is fairly robust to C m . However, it performs poorly when β is misspecified, and increasingly so as m increases. The results is markedly different for FDR: FDR thresholding is less adaptive to C m than Bayes0, but much more adaptive to the sparsity parameter β, as illustrated by the fact that the configurations with low relative excess risk span the whole range of β. For FDR thresholding at α m = α opt m (β 0 , C 0 ), the fraction of configurations (β, C m ) for which E m ≤ 0.1 increases as m increases. This illustrates the asymptotic optimality of FDR thresholding for classification under sparsity in the Gaussian location model. Similar conclusions hold for the Gaussian scale and the Laplace scale models, see the Supplementary Materials of[J17]).

Figure 4 . 1 :

 41 Figure 4.1: Adaptation to sparsity by FDR thresholding in the Gaussian location model. relative excess risks E m for various thresholding procedures (rows) and different values of m (columns). In each panel, the corresponding risk is plotted as a function of β ∈ [0, 1] (horizontal axis) and C m ∈ [0, 1] (vertical axis). Colors range from white (low risk) to dark red (high risk), as indicated by the color bar at the bottom. Black lines represent the level set E m = 0.1. The point (β, C m ) = (β 0 , C 0 ) is marked by "+". We chose β 0 = 1/2 and C 0 = 1/2. See main text for details.

Figure 5 . 1 :

 51 Figure 5.1: Volcano plot.

  of[START_REF] Goeman | Multiple testing for exploratory research[END_REF] in Figure5.2, with m = 3 and where the hypotheses rejected by the local test are crossed. In this particular example, all of the hypotheses rejected by the local test are also rejected by the closed testing procedure. The rejection of H {2,3} = H 0,2 ∩ H 0,3 is non-consonant, because neither H {2} = H 0,2 or H {3} = H 0,3 is rejected by the closed testing procedure. The bound V GS α tells us not only that V GS α ({1}) = 0 (which we already knew from the classical FWER control), but also that V GS α ({2, 3}) = 1, that is, at least H 2 or H 3 is a true positive.

Figure 5 . 2 :

 52 Figure 5.2: [81, Figure 1]: "Intersection hypotheses formed by elementary hypotheses H 1 , H 2 and H 3 . Rejected hypotheses have been marked with a cross. The rejection of H 2 ∩ H 3 is a non-consonant rejection."

Figure 5 . 3 :

 53 Figure 5.3: Left: volcano plot using post hoc inference; right: corresponding FDP confidence envelope.

Definition 5 . 4 (

 54 Threshold template). A template is a family of functions t k (λ), λ ∈ [0, 1], k ∈ N K , such that K ∈ N mand for all k ∈ N K , t k (0) = 0 and t k (•) is non-decreasing and left-continuous on [0, 1]. The parameter K is called the size of the template.

Figure 5 Figure 5 . 4 :

 554 Figure 5.4: Curves k → t k (λ) for a wide range of λ values. Left: linear template. Right: beta template.

Figure 5 . 5 :

 55 Figure 5.5: JER control based on the linear template for equi-correlated test statistics: λ-calibration by sign-flipping in the location model.

Figure 5 .Figure 5 . 6 :

 556 Figure 5.6 compares confidence envelopes obtained from the Simes reference family (longdashed purple curve), to the permutation-based λ-calibration derived from Theorem 5.6 using B = 1, 000 permutation of the sample labels for the linear template with K = m (dashed red curve) and the beta template with K = 50 (solid green curve). Note that Assumption (Rand) holds in this two-sample framework. While (1 -α)-lower confidence bounds on the number of true positives of the form k, |S k | -V (S k ) : k ∈ N m are displayed in the left panel, (1-α)-upper confidence bounds on the proportion of false positives k, V (S k )/ |S k | : k ∈ N m are shown in the right panel. The vertical line in Figure 5.6 corresponds to the 163 genes selected by the BH procedure

Figure 5 . 7 :

 57 Figure 5.7: Volcano plot using adaptive post hoc inference.

Figure 7 . 64 (

 764 Figure 7.1(b) summarizes a simulation run with n = 200 independent copies of the synthetic observed data structure based on two specific GBM samples for the EGFR gene.Such a simulation study provides an arguably realistic use case where closed-form expressions for the features of interest θ(P ), µ(P ), g(P ), and σ 2 (P ) can be derived, and the value of Ψ(P ) can be estimated accurately. We estimated this value by ψ B (P ) using B = 10 5 Monte-Carlo samples. Gray rectangles represent 95%-accuracy intervals [ψ B (P ) ± ξ 0.975 s B (P )/√ n] and [ψ B (P ) ± ξ 0.975 s B (P )/ √ B] for the true parameter Ψ(P ) based on n observed data structures (light gray) and B = 10 5 observed data structures (dark gray). Here, ξ a is the 1 -a quantile of the standard Gaussian distribution, and s B (P ) = Var P B (∇Ψ(P )(O)) is a

Figure 7 . 1 :

 71 Figure 7.1: Illustrating DNA methylation, DNA copy number, and gene expression data. In both graphics, we represent kernel density estimates (diagonal panels), pairwise plots (lower panels), and report the pairwise Pearson correlation coefficients (upper panels). (a). Real dataset corresponding to the EGFR gene in 187 GBM tumor samples. For 130 among the 187 samples, only DNA copy number and gene expression data were available (circles in lower middle plot). (b). Simulated dataset consisting of n = 200 independent copies of the synthetic observed data structure described in the main text. Note that the constant O X

  Convergence and robustness. The substantial bias in the initial estimation is diminished (if not perfectly corrected) at the first updating step of the TMLE procedure, illustrating the robustness of the targeted estimator. The empirical distributions of {ψ k n,b : b ≤ B } for k = 1, 2, 3 are not (visually) markedly different, an empirical indication that the TMLE procedure converges quickly. Asymptotic normality, and coverage. The asymptotic normality of the TMLE estimator was confirmed by Lilliefors and Kolmogorov-Smirnov tests. We used the asymptotic variance estimate (s k n,b ) 2 = Var P n,b ∇Ψ(P k n,b )(O) to form confidence intervals. Contrary to s 2

Figure 7 . 2 :

 72 Figure 7.2: Empirical distribution of {ψ k n,b : b ≤ B } based on n = 200 independent observed data structures for k = 0 (initial estimator) and k iterations of the updating procedure (k = 1, 2, 3), as obtained from B = 10 3 independent replications of the simulation study. Gray rectangles represent 95%-accuracy intervals for the true parameter based on n observed data structures (light gray) and B = 10 5 observed data structures (dark gray).

Figure 8 . 1 :

 81 Figure 8.1: Expression measurements of the PIAS2 and UBE2C genes in two patient populations (BCR/ABL and NEG) from the Leukemia data set. The marginal distributions are similar across groups for its genes, but the scatter plot reveals differential expression in two dimensions.

Figure 8 . 2 :

 82 Figure 8.2: Eigenvectors of the graph with adjacency matrix defined in (8.2).

2 Figure 8 . 5 :

 285 Figure 8.4: Left: Increase in distribution shift required for Hotelling's T 2 -test to maintain a given power when increasing the number of tested new coefficients: ∆2 k+l -∆ 2 k vs. l such that β α,k+l (∆ 2 k+l ) = β α,k (∆ 2 k ). Power β α,k+l (∆ 2 k+l ) computed under the non-central F -distribution F n1n2 n1+n2 ∆ 2 k+l ; k + l, n 1 + n 2 -(k + l) -1 , for n 1 = n 2 =100 observations, k = 5, and α = 10 -2 . Each line corresponds to the fixed shift ∆ 2 k and power β α,k (∆ 2 k ) pair indicated in the legend. Right: Zoom on the first 30 dimensions.

Figure 9 . 1 :

 91 Figure 9.1: Illustration of the main DNA copy-number states in three cancer samples: normal, gain, deletion, copy-neutral LOH. Illustration by Henrik Bengtsson.

Figure 9 . 2 :

 92 Figure 9.2: Copy-number signals in two genomic regions (left and right columns). Orange box: allelic signals after Tumorboost normalization (last row) have a much higher signal to noise ratio than before Tumorboost normalization (third row).

9. 3 Figure 9 . 3 :

 393 Figure 9.3: Copy-number profile of a tumor sample. Top: total copy numbers (c): the same panel is displayed twice. Bottom: allelic ratios (b) on the left, and decrease in heterozygosity (d) on the right. Red vertical bars indicate the presence of change points. Gray and black dots correspond to SNP that were called homozygous and heterozygous, respectively, in a paired normal sample.

  γ j ) 2 s.t. Dγ 1 ≤ v and γ -2 1 ≤ u . (9.4)

Figure 9 .

 9 Figure 9.4 displays microarray copy-number profiles (c, b) for two chromosomes (in columns) in two samples (in rows) from the tumor cell line H1395, whose karyotype is displayed in Figure1.1. These samples are part of a dilution series[START_REF] Rasmussen | Allele-specific copy number analysis of tumor samples with aneuploidy and tumor heterogeneity[END_REF] available from the Gene

Figure 9 . 4 :

 94 Figure 9.4: (c, b) signals for two annotated regions from tumor cell line NCI-H1395-4W. Left: gain on chromosome 5; right: loss on chromosome 6. Top: sample with pure tumor; Bottom: sample with a mixture of 70% tumor cells and 30% normal cells.

Figure 9 . 6 (Figure 9 . 6 :

 9696 Figure 9.6: Comparing the accuracy of change point detection of different change point methods. Left: ROC curves; right: an example of sample illustrating the simulation scenario. 2d methods are in solid lines, 1d methods are in dashed lines.

. 1 )

 1 where • is the Euclidean norm in R d , xG = |G| -1 i∈G x i is the center of gravity of G and |G| denotes the cardinal of G. The loss of information when merging two disjoint clusters G and G into G ∪ G is quantified by : δ(G, G ) := I(G ∪ G ) -I(G) -I(G ).

Figure 10 . 1 :

 101 Figure 10.1: Example of forward pencils (in yellow and green) and backward pencils (in green and blue), and illustration of Equation (10.7) for cluster C = {i, . . . , j -1}. Left: cluster smaller than bandwidth (k ≤ h); right: cluster larger than bandwidth k ≥ h.

Figure 10 . 2 :

 102 Figure 10.2: Left: quality of the band approximation as a function of the bandwidth h for five different chromosomes. Right: computation times versus p for LD matrices, and clustering with rioja and adjclust.

Figure 10 . 3 :

 103 Figure 10.3: A linkage disequilibrium (r 2 ) plot with the inferred block structures (black and red contour lines) for a set of 68 contiguous SNPs located on the MHC region. Left:within the structure inferred by the proposed method, the blocks selected by the Group Lasso are delimited with a red contour line. The SNPs selected by single marker analysis are plotted with a red star (*), and the SNPs missed by Lasso with a blue dash (-). Right: within the structure inferred by the haplotype association method, the blocks selected by the competing method are delimited with a red contour line.

  Contents 2.1 Statistical setting . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2 Family-Wise Error Rate control . . . . . . . . . . . . . . . . . . 15 2.3 False Discovery Rate control . . . . . . . . . . . . . . . . . . . . 18 2.4 From selective inference to simultaneous inference . . . . . . 21 2.5 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

Table 2

 2

			.1.
		H 0 not rejected H 0 rejected
	H 0 true	true negative	false positive
	H 0 false	false negative	true positive

Table 2 .

 2

1: Qualification of the four possible outcomes of a statistical test.

Table 2 .

 2

	Setting	Procedures
	General dependence Bonferroni	Holm
	PRDS(H 0 )	Hochberg	Hommel
	indep	Šidák

2: Summary of FWER-controlling procedures

Table 5 .

 5 1: Conservativeness of Simes inequality in the Gaussian equi-correlated model. Here, |H 0

	-correlation level: ρ	0	0.1	0.2	0.4	0.8
	Size (achieved level)	0.99 0.85 0.72 0.42 0.39

  Figure 8.4: Left: Increase in distribution shift required for Hotelling's T 2 -test to maintain a given power when increasing the number of tested new coefficients: ∆ 2 k+l -∆ 2 k vs. l such that β α,k+l (∆ 2 k+l ) = β α,k (∆ 2 k ). Power β α,k+l (∆ 2 k+l ) computed under the non-central

	Increase in ∆ 2 required to maintain power	0 1 2 3 4	50 ∆ 2 =0.2 (power=0.42) 100 ∆ 2 =0.15 (power=0.28) ∆ 2 =0.1 (power=0.16) ∆ 2 =0.05 (power=0.06)	150	Increase in ∆ 2 required to maintain power	0.00 0.05 0.10 0.15 0.20 0.25	0	5 ∆ 2 =0.2 (power=0.42) 10 15 ∆ 2 =0.15 (power=0.28) ∆ 2 =0.1 (power=0.16) ∆ 2 =0.05 (power=0.06)	20	25	30
			Number of graph-based dimensions added				Number of graph-based dimensions added

Source: https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/ tcga, retrieved on October 8,

2 See e.g. the discussion Collect Data First, Ask Questions Later from the podcast "The Effort Report".

This data set will be used in Chapters 2, 5 and 8.

These data will be used in Chapter 9

According to scholar.google.com these papers have been cited

[START_REF] Goeman | Multiple hypothesis testing in genomics[END_REF],000 and 30,000 times, respectively, as of December 2019. Both are among the most cited statistics paper, and Yoav Benjamini has been awarded the 2019 Karl Pearson Prize from the International Statistics Institute for the FDR paper.

For notation simplicity, we also follow[START_REF] Casella | Statistical inference[END_REF] in using the same letter P for probability measures on the observation space (X , X) and on the implicitly defined domain of definition of X.

The results stated in this chapter under PRDS in fact hold for a slightly larger class of distribution called weak PDRS. We refer to[START_REF] Roquain | Contributions to multiple testing theory for high-dimensional data[END_REF] for references on this distinction, which is not essential for our purpose.

However, both of them can also be seen as conservative, as discussed and illustrated in Section 5.4.

Formally we should write ∀i ∈ H, {i} ∈ H but we identify i with {i} in order to alleviate notation.

The closed testing principle may also be used to prove the FWER control of the Hommel procedure[START_REF] Hommel | A stagewise rejective multiple test procedure based on a modified Bonferroni test[END_REF] under (PRDS(H 0 )). This procedure is slightly more powerful but also more complicated to define than the Hochberg procedure, see[START_REF] Dickhaus | Simultaneous statistical inference[END_REF] Definition 5.4].

Chapter 3 provides theoretical results supporting this observation.

In the context of analysis of variance Scheffé[START_REF] Scheffé | The analysis of variance[END_REF], a contrast is simply a linear combination of the group means such that the coefficients of the combination sum to 0. This naturally covers the test for equality of means of two subgroups.

The exact translation of the definition in(2.14) incorporates an additional 1/m term in the numerator of Π λ . This term is required to ensure finite sample FDR control, but the two associated FDR controlling procedures are asymptotically equivalent, in a sense formally defined in[J24].

Here, we consider Hadamard-differentiability on the set of càdlàg functions on [0,1], tangentially to the set of continuous functions on [0, 1]. We refer to[START_REF] Van Der | Cambridge Series in Statistical and Probabilistic Mathematics[END_REF] for a formal definition of Hadamard differentiability.

We refer to the "vanish-second-order" condition in[START_REF] Delattre | On empirical distribution function of high-dimensional Gaussian vector components with an application to multiple testing[END_REF] for a formal definition of this condition.

See [J2, Section

2.1 and Appendix S-4] for proofs.

In this particular case, λ-calibration with the linear template yields λ(α) > α, which a posteriori implies that the Simes inequality was indeed valid.

LICORN has been developed by Mohamed Elati, and I participated to this development during my PhD thesis.

http://cran.r-project.org/package=tmle.npvi

The data are available from the TCGA at: https://tcga-data.nci.nih.gov/docs/publications/brca 2012.

Sometimes referred to as hypergeometric test in the bioinformatics literature.

Loss of heterozygosity (LOH) is the loss of the contribution of one parent in a genomic region. It includes the case of a deletion, but also of copy-neutral LOH, corresponding to two copies from the same parent.

http://cran.r-project.org/package=jointseg

http://cran.r-project.org/package=acnr

available on CRAN at https://cran.r-project.org/package=rioja.

The proof is a consequence of the fact that r 2 can be formulated as the squared correlation between (latent) variables that depend on haplotypes.

"Is There a Reproducibility Crisis in Science?". Nature Video, Scientific American.

May 2016.

STRATOS stands for STRengthening Analytical Thinking for Observational Studies, see http://www. stratos-initiative.org/.

This is a clear improvement with respect to the standards of 10-15 years ago, where implementations were generally "available upon request" to the authors. This progress is a consequence of (i) the more and more frequent requirements of journals to make code truly available, (ii) the development of efficient packages for package development[START_REF] Wickham | devtools: Tools to Make Developing R Packages Easier[END_REF], testing[START_REF] Wickham | testthat: Get Started with Testing[END_REF], and documentation[START_REF] Wickham | roxygen2: In-Line Documentation for R[END_REF], and (iii) the recognition of R packages as a contribution in the applied statistic community, which can even lead to dedicated publications[J8].

To address this difficulty, the Comprehensive R Archive Network (CRAN) offers Task Views, which give a list of packages (curated by specialists) pertaining to a given topic. For example, the recent task view on missing data lists more than 130 CRAN packages related to this particular topic. Also, the developers of Bioconductor packages, which are specifically dedicated to computational biology, are encouraged to use standard methods and classes.

See the blog post Getting a big scientific prize for open-source software by Gaël Varoquaux.

This app is available at https://pneuvial.shinyapps.io/posthoc-bounds_ordered-hypotheses/.

Remerciements

Chapter 13

Bibliographic references

Algorithm. A closed-form expression of ∇Ψ can be derived [J14, Proposition 1] for the parameter Ψ defined in (7.1). The expression of ∇Ψ involves finite-and infinite-dimensional features of P , notably including:

Moreover, it can be shown that ∇Ψ is double-robust, in the sense that for any (P, P ) ∈ M 2 , we have Ψ(P ) = Ψ(P ) as soon as either (µ(P ) = µ(P ) and g(P ) = g(P )) or θ(P )(0, •) = θ(P )(0, •).

Let us assume that we observe n independent random variables O 1 , . . . , O n drawn from P 0 . By (7.2), building an initial substitution estimator Ψ(P 0 n ) of Ψ(P 0 ) requires the estimation of θ P0 , of σ 2 P0 , and of the marginal distribution of (W, X) under P 0 . It can be shown that the latter can itself be obtained from estimates of g P0 and µ P0 , and the marginal distribution of W under P 0 n , which itself is simply estimated by its empirical counterpart. Therefore, using Monte-Carlo estimation, we can obtain an initial estimator of Ψ(P 0 n ) from initial estimates (θ 0 n , g 0 n , µ 0 n , σ 0 n ). The k th update step (for k ≥ 0) consists in estimating the efficient influence curve ∇Ψ P k n based on the features θ k n , µ k n , g k n and σ k n . This estimation step defines a one-dimensional model

This model is called a fluctuation of the law P k n . An updated estimator of P 0 is then defined by

where ε k n is the maximum likelihood estimator of the fluctuation model. Finally, updated estimations of the features and of the parameter: ψ k+1 n are obtained similarly to the initialization step, i.e. via Monte-Carlo estimation of the marginal law of (X, W ). The series of updates is interrupted if the total variation distance d T V (P k n , P k+1 n

Finally, the TMLE ψ * n is defined as ψ * n = lim k→∞ ψ k n , assuming that the limit exists, or more generally as ψ kn n for a conveniently chosen sequence {k n } n≥0 .

This estimation method has been implemented in the R package tmle.npvi [S5] which is available from CRAN 1 . The implementation is described in detail in [J8]. This implementation incorporates estimators of the features (θ, g, µ, σ) based on classical generalized linear models, as well as estimators relying on super-learning, a generic method of aggregation of estimators proposed by [START_REF] Van Der Laan | Super Learner[END_REF].

Inference

Convergence of the iterative procedure. Lemmas 3 and 4 in [J14] ensure that if the sequence (parametrized by the iteration index k) P n ∇Ψ(P k n ) 2 is bounded away from 0, then the above iterative procedure converges, in the sense that the sequence (ε k n ) k>0 tends to 0. If moreover the series k |ε k n | is convergent, then the sequence of estimators ψ k n = Ψ(P k n ) is also convergent.

Asymptotic properties of the proposed estimator.

The TMLE is consistent as soon as the estimators of the features θ, µ, σ and g are convergent, and that one of the estimator of θ(0, •) or the estimators of µ and of g is consistent [J14, Proposition 2]. This remarkable property is inherited from the double robustness of the efficient influence curve ∇Ψ, combined with the fact that the TMLE solves the efficient influence curve equation (i.e., P n ∇Ψ(P kn n ) ≈ 0). The TMLE satisfies a central limit theorem [J14, Proposition 3] where A + denotes the generalized inverse of a matrix A, the p × k matrix U [k] denotes the restriction of U to its first k columns, and 1 k is a p × p diagonal matrix, with ith diagonal element equal to one if i ≤ k and zero otherwise. Note that, as retaining the first k dimensions corresponds to a non-invertible transformation, this filtering indeed has an effect on the test statistic, that is, we have T 2 k = T 2 in general. Lemma 1 stated below shows that gains in power can be achieved by filtering, under the assumption of a smooth shift along the graph. We let δ = U δ and Σ = U ΣU denote, respectively, the mean shift and covariance matrix in the new space. Given

the distribution shift restricted to the first k dimensions of δ and Σ, i.e., based on only the first k elements of δ, (δ i : i ≤ k), and the first k × k diagonal block of Σ, (σ ij : i, j ≤ k).

Lemma 8.1. For any level α and any 1 < l ≤ pk, there exists η(α, k, l) > 0 such that

where

. Under the assumption that the distribution shift is smooth, i.e., lies mostly in the first few graph-based coefficients, so that ∆ 2 k ( δ, Σ) is nearly maximal for a small value of k, Lemma 1 states that performing Hotelling's test in the new space restricted to its first k components yields more power than testing in the entire new space. Equivalently, the test is more powerful in the original space after filtering than in the original unfiltered space. The increase in shift η(α, k, l) required to maintain power when increasing dimension can be evaluated numerically for any (α, k, l). Note that this result holds because retaining the first k new components is a non-invertible transformation.

Corollary 8.2 states that if the distribution shift lies in the first k new coefficients, then testing in this subspace yields strictly more power than using additional coefficients:

In particular, if there exists k < p such that δj = 0 ∀ j > k (i.e., the mean shift is smooth) and Σ is block-diagonal such that σij = 0 ∀ i < k, j > k, then gains in power are obtained by testing in the first k new components. Although non-necessary, this condition is plausible when the mean shift lies at the beginning of the spectrum (i.e., has low energy), as the coefficients which do not contain the shift are not expected to be correlated with the ones that do contain it. Figure 8.4 illustrates, under different settings, the increase in distribution shift necessary to maintain a given power level against the number of added coefficients. Under the assumption of block-diagonal covariance, it is also possible to directly relate the energy of the mean shift vector to the gain in power, see Corollary 2 in [J16, Supplement A].

Numerical experiments

The empirical performance of the graph-structured test has been assessed in numerical experiments where the distribution shift ∆ 2 satisfies the above-defined smoothness assumptions. More specifically, we assume that this shift lies in the first k 0 -dimensional eigen space of the graph. As expected in this favorable setting, Figure 8.5 shows (in the case of a block diagonal covariance structure and for k 0 = 3) that our proposed T 2 -statistic in the first k 0 graph-based coefficients (dashed red lines) compares favorably to several alternatives: the standard Hotelling T 2 -statistic in the original space and T 2 -statistic in the first k 0 principal components (left panel), the statistics of [START_REF] Bai | Effect of high dimension : by an example of a two sample problem[END_REF] (BS), [START_REF] Chen | A two-sample test for high-dimensional data with applications to gene-set testing[END_REF] (CQ), and [START_REF] Srivastava | A test for the mean vector with fewer observations than the dimension[END_REF] (SD) (middle panel), and the Adaptive Neyman statistics of [START_REF] Fan | Test of Significance when data are curves[END_REF] (right panel).

Fast segmentation of a band similarity matrix

When the p objects to be clustered belong to R d , with d < p, the computation of Ward's linkage between two clusters can be done in O(d) by exploiting its explicit alternative formulation as the distance between centers of gravity (10.3). In such cases, it is possible to obtain unconstrained HAC in O(p 2 log 2 p) in time [START_REF] Eppstein | Fast hierarchical clustering and other applications of dynamic closest pairs[END_REF], and lower complexities could possibly be achieved for adjacency-constrained HAC. However, we focus in this chapter in the situation where the input objects are represented by pairwise similarities. In such cases there is generally no explicit or finite-dimensional representation of the centers of gravity, and the time complexity of adjacency-constrained HAC is intrinsically quadratic in p because all of the p 2 similarities are used to compute all of the required linkage values required by Algorithm 2. This is the case for the implementation of adjacency-constrained HAD provided in the CRAN/R package rioja [START_REF] Juggins | rioja: Analysis of Quaternary Science Data[END_REF].

Band similarity assumption. In applications where adjacency-constrained clustering is relevant, such as Hi-C and GWAS data analysis, this quadratic time complexity is a major practical bottleneck because p is typically of the order of 10 4 to 10 5 for each chromosome.

Fortunately, in such applications it also makes sense to assume that the similarity between physically distant objects is small. Specifically, we assume that S is a band matrix of bandwidth h + 1, where h ∈ {1 . . . p}: s ij = 0 for |i -j| ≥ h. This assumption is not restrictive, as it is always fulfilled for h = p. However, we will be mostly interested in the case where h p. We now describe an algorithm proposed in [J5] to perform adjacency-constrained HAC under this assumption. This algorithm relies on (i) constant-time calculation of each of the Ward linkages involved in Algorithm 2 using Equation (10.6), and (ii) storage of the candidate fusions in a min-heap.

Ward's linkage as a function of pre-calculated sums. We have shown in [J5] that the sum of all similarities in any cluster C = {i, . . . , j -1} of size k = ji can easily be obtained from sums of elements in the first min(h, k) sub-diagonals of S: Lemma 10.1. For 1 ≤ r, l ≤ p, let P (r, l) = {s ij : 1 ≤ i, j ≤ r, |i -j| < l} and P (r, l) = P (p + 1 -r, l). Letting h k := min(h, k), we have

P (p, h k ) is the "full" pencil of bandwidth h k (which also corresponds to P (1, h k )). Lemma 10.1 is illustrated in Figure 10.1, with r ∈ {i, j}, for l = k ≤ h in the left panel, and l = h ≤ k in the right panel. In both panels, P (j, min(h, k)) is the sum of elements in the yellow and green regions, while P (i, min(h, k)) is the sum of elements in the green and blue regions. Because P and P are sums of elements in pencil-shaped areas, we call P (r, l) a forward pencil and P (r, l) a backward pencil.

Lemma 10.1 makes it possible to compute δ(C, C ) in constant time from the pencil sums using Equation (10.6). By construction, all the bandwidths of the pencils involved are less than h. Therefore, only pencils P (r, l) and P (r, l) with 1 ≤ r ≤ p and 1 ≤ l ≤ h have to be pre-computed, so that the total number of pencils to compute and store is less than 2ph. These computations can be performed recursively in a O(ph) time complexity.

Storing candidate fusions in a min-heap. Iteration t of Algorithm 2 consists in finding the minimum of pt elements, corresponding to the candidate fusions between the pt + 1 clusters in C t-1 , and merging the corresponding clusters. Storing the candidate fusions in an unordered array and calculating the minimum at each step would mean a quadratic time complexity. One intuitive strategy would be to make use of the fact that all but 2 to 3 candidate fusions at step t are still candidate fusions at step t -1. However, maintaining a totally-ordered list of candidate fusions is not efficient because the cost of deleting and Chapter 11

Directions for future research

The previous chapters summarize some attempts to address the complexity of genomic data by considering their sparsity, heterogeneity and structure as constraints that alleviate the curse of dimensionality, and guide statistical methods toward solutions that are more plausible from a biological standpoint. My recent research contributions, particularly in the context of the SCALES project (2017-2019) funded by the CNRS MITI and the SansSouci project (2016-2020) funded by ANR have also raised a number of new research questions, that cover a broad spectrum from theory to applications. In this chapter we describe some research directions toward exploiting the multiscale nature of biological problems, and providing inferential guarantees for exploratory findings. We will pay particular attention to the evaluation of the statistical and computational "price to pay" for accounting for these important features. Some open questions in the field of post hoc inference are discussed in Section 11.1. In Section 11.2 specific methodological challenges for inference for signals that are structured along the genome or in three dimensions are described. Section 11.3 is an opening to broader scientific challenges which I consider to be of primary importance, especially in the context of interdisciplinary research. Journal papers