Biological Context

All biological traits are influenced by both genetic and environmental factors. Quantifying these two contributions for a particular trait is a fundamental and challenging question in biology. The concept of heritability refers to the part of the variability of an observed trait (or phenotype) which can be attributed to genetic causes. Several missconceptions regarding heritability are due to the use of the term in the common language, which differs from the technical definition in the genetic field. For instance, a frequent missconception would be that heritability is the proportion of a phenotype that is transmitted to the next generation. Firstly, genes are passed on from parents to offspring but phenotypes are not. Secondly, if half of the genetic effects are indeed transmitted from each parent, this particular half is specific to each offspring. [START_REF] Visscher | Heritability in the genomics era -concepts and misconceptions[END_REF] gathered these frequent questions and mistakes regarding heritability. The concept of heritability as it is used in the field of genetics is presented in the following section.

Definition of heritability

As elegantly explained by [START_REF] Visscher | Heritability in the genomics era -concepts and misconceptions[END_REF], we consider the simple modeling where a phenotype of interest is the result of genetic and environmental effects considered as independent: Phenotype (P) = Genotype (G) + Environment (E).

The variance of the observable phenotypes (σ 2 P ) can then be expressed as a sum of unobserved underlying variances (σ 2 G and σ 2 E ):

σ 2 P = σ 2 G + σ 2 E .
Heritability (H 2 ) is defined as a ratio of variances and expresses the proportion of the phenotypic variance that can be attributed to genetic factors:

H 2 = σ 2 G σ 2 P .
The genetic variability may be partitioned into variances coming from different sources, in particular the variance σ 2 A of additive genetic effects. Such additive effects are characterized by the impact of single nucleotide polymorphisms (SNPs), which are DNA sequence differences at the positions of the genome where there exists considerable variability in the population. These positions are actually not frequent compared to the totality of the genome: the human genome is indeed composed of approximately 3 billions of base pairs, a very large fraction of which are identical for all humans. In the sequel, we will consider the "narrow sense heritability" which is the proportion of variability explained only by additive genetic effects, defined by

h 2 = σ 2 A σ 2 P .
Since the access to the genotype of thousands of individuals has been made possible by the spectacular decreased cost of DNA sequencing, the heritability of quantitative traits and pathologies has become widely studied. [START_REF] Yang | Common snps explain a large proportion of the heritability for human height[END_REF] estimated for instance that around 45 % of human height was explained only by the most frequent SNPs.

Heritability in human genetics

We will here motivate the estimation of heritability for human traits. It is indeed a step toward the understanding of complex diseases, which have often multiple causes. We refer in particular to diseases which are not caused by a single affected gene but are nevertheless suspected to have a strong genetic component, probably split among different genes. For instance, the causes of psychiatric disorders, such as autism or schizophrenia, remain vague. A genetic component has been suggested by the results of monozygotic and dizygotic twin studies (monozygotic twins have identical genomes while dizygotic twins share around 50% of their genomes). These studies show that if one twin is affected by autistic disorders, the other one is also affected in 82 to 92 % of cases for monozygotic twins [START_REF] Bailey | Autism as a strongly genetic disorder: evidence from a british twin study[END_REF]), or in 20% of cases for dizygotic twins [START_REF] Hallmayer | Genetic heritability and shared environmental factors among twin pairs with autism[END_REF]. Moreover, for a family which already has an autistic child, the risk of having another one is evaluated at 20 % against 1% in the general population. These studies describe autism as the psychiatric disease with the most important genetic component. However, the severity of the autistic traits (language and interaction disorders, intellectual disabilities...) can be very different for two patients with similar causes, for instance the same mutation. It would thus seem, as it is also the case for other genetic diseases, that the genetic background modulates the effect of a causal mutation and renders an individual more or less sensitive to developing autistic traits. Furthermore, all studies show that despite their identical genetic patrimony, the concordance of symptoms of monozygotic twins is never total, which confirms an epigenetic and/or environmental component. However, quantifying these different possible causes and potential interactions between them remains a challenging issue.

Determining a significant genetic component of a disease also constitutes a strong argument to refute some popular beliefs about causes of some illnesses. For instance, an important wave of anti-vaccine movement has been fueled by a presumed connection between the hepatitis B vaccine and multiple sclerosis. Similarly, the measles vaccine has been accused to cause autism [START_REF] Uno | The combined measles, mumps, and rubella vaccines and the total number of vaccines are not associated with development of autism spectrum disorder: the first case-control study in asia[END_REF]. Although no link has ever been demonstrated [START_REF] Poland | Understanding those who do not understand: a brief review of the anti-vaccine movement[END_REF], the consequences of the fact that many parents refuse to vaccinate their children remains a major public health issue. Indeed, a recent study [START_REF] Uno | The combined measles, mumps, and rubella vaccines and the total number of vaccines are not associated with development of autism spectrum disorder: the first case-control study in asia[END_REF] showed that more than 25% of parents in the US refused to vaccinate their children against mortal diseases like measles. Regarding other proposed causes of autistic disorders, the "refrigerator mother theory" was developed by the psychiatrist Leo Kanner who claimed to observe a "genuine lack of maternal warmth" among his patients' mothers. Even though this theory has since been discarded, the mothers of autistic patients have suffered severe and unwarranted accusations for several decades.

Heritability in vegetal and animal genetics

In the field of vegetal and animal genetics, heritability estimation is the first step to the selection of traits of interest, generally related to the yield of a valuable resource. We can mention the examples of the optimization of the yield of milk, [START_REF] Visscher | Genetic parameters for milk yield, survival, workability, and type traits for australian dairy cattle[END_REF] or wheat [START_REF] Eid | Estimation of heritability and genetic advance of yield traits in wheat (Triticum aestivum L.) under drought condition[END_REF]. The goal in [START_REF] Eid | Estimation of heritability and genetic advance of yield traits in wheat (Triticum aestivum L.) under drought condition[END_REF] is to determine strongly heritable traits related to the yield and then to obtain an optimal genotype. This genotype was even selected to be the most resistant to extreme environmental conditions like water deprivation, which is currently a fundamental issue.

If this kind of practice is generally accepted in animal genetics, it creates a controversy on possible consequences of heritability estimations of human traits. Several studies estimated the IQ heritability [START_REF] Toro | Genomic architecture of human neuroanatomical diversity[END_REF] and [START_REF] Davies | Genome-wide association studies establish that human intelligence is highly heritable and polygenic[END_REF] even announced that "Genome-wide association studies establish that human intelligence is highly heritable and polygenic". The controversy about IQ heritability is discussed in [START_REF] Visscher | Heritability in the genomics era -concepts and misconceptions[END_REF], who enumerates reasons for the polemic nature of this issue. These include the very controversial definition of IQ as a measure of intelligence as well as historical abuses related to eugenics. We will not further discuss this controversy here, we just mention it to illustrate a frequent issue when dealing with heritability of human features.

Having briefly argued for the general interest of estimating heritability, we will now present the statistical modeling used to provide these estimations.

1.2 Heritability estimations in high dimensional linear mixed models 1.2.1 State of the art Linear Mixed Models (LMMs) have been widely used in several fields, especially in medicine and genetics. [START_REF] Yang | Common snps explain a large proportion of the heritability for human height[END_REF] proposed to estimate the heritability of human height using a classical LMM defined as follows:

Y = Xβ + Zu + e (1.1)
where Y = (Y 1 , . . . , Y n ) is the vector of observations of a phenotype of interest, X is a n × p matrix of predictors (or fixed effects), β is a p × 1 vector containing the unknown linear effects of the predictors, and u and e correspond to the Gaussian random effects with variances respectively equal to σ 2 u and σ 2 e . Moreover, Z is a n × N matrix which contains the genetic information. More precisely, the Z i,j 's are normalized random variables in the following sense: they are defined from a matrix W = (W i,j ) 1≤i≤n, 1≤j≤N by Z i,j = W i,j -W j s j , i = 1, . . . , n, j = 1, . . . , N ,

where

W j = 1 n n i=1 W i,j , s 2 j = 1 n n i=1
(W i,j -W j ) 2 , j = 1, . . . , N .

(1.3)

In (1.2) and (1.3) the W i,j 's are such that for each j in {1, . . . , N } the (W i,j ) 1≤i≤n are independent and identically distributed random variables and such that the columns of W are independent. In genetic applications, the matrix W contains all the genetic information about all the individuals in the study. With this definition the columns of Z are empirically centered and have an empirical variance equal to 1.

The LMM appears to be an intuitive modeling to describe the biological concept of heritability as a ratio of genetic and phenotypic variances. [START_REF] Yang | Common snps explain a large proportion of the heritability for human height[END_REF] and [START_REF] Pirinen | Efficient computation with a linear mixed model on large-scale data sets with applications to genetic studies[END_REF] proposed to estimate the parameter .4) commonly considered as the mathematical definition for heritability since it determines how the variance is shared between u and e.

η = N σ 2 u N σ 2 u + σ 2 e , ( 1
In Model (1.1), the log-likelihood conditionnaly to Z is given by:

L(β, σ 2 u , σ 2 e ) = - n 2 log(2π) - 1 2 log(|ZZ σ 2 u + σ 2 e Id R n | - 1 2 (Y -Xβ) (ZZ σ 2 u + σ 2 e Id R n ) -1 (Y -Xβ).
(1.5) [START_REF] Searle | Variance Components. Wiley Series in Probability and Statistics[END_REF] gathered plenty of optimization techniques to estimate the parameters β, σ 2 u and σ 2 e , among which we can quote Henderson equations or iterative methods like Fisher-Scoring and Newton-Raphson. A natural idea to estimate heritability is to estimate the variance parameters σ 2 u and σ 2 e in order to obtain an estimator as the ratio:

N σu 2 /(N σu 2 + σe 2 ). [START_REF] Pirinen | Efficient computation with a linear mixed model on large-scale data sets with applications to genetic studies[END_REF] noticed that the model defined in (1.1) could be reparameterized with β, η and σ 2 = N σ 2 u + σ 2 e as new parameters. More precisely,

Y ∼ N Xβ, η σ 2 R + (1 -η )σ 2 Id R n ,
where R = ZZ /N . Let U be the orthogonal matrix (U U = UU = Id R n ) such that URU = diag(λ 1 , . . . , λ n ) is a diagonal matrix having its diagonal entries equal to λ 1 , . . . , λ n . Hence, Y = U Y is a zero-mean Gaussian vector having a covariance matrix equal to diag(η σ 2 λ 1 + (1 -η )σ 2 , . . . , η σ 2 λ n + (1 -η )σ 2 ), where the λ i 's are the eigenvalues of R. Let us also denote X = U X. Finally they computed and maximized the log-likelihood:

L n (β, σ 2 , η) = - n 2 log(σ 2 ) - 1 2 n i=1 log(η(λ i -1) + 1) - 1 2σ 2 n i=1
( Y i -Xβ) 2 η(λ i -1) + 1 -n 2 log(2π), (1.6) where Y = ( Y 1 , ..., Y n ).

The aforementioned approaches raise two main concerns: firstly, they all have been validated in the framework where N is fixed and n goes to infinity. Indeed, using classical results of the LMM, we can obtain properties of consistency and asymptotic normality for the maximum likelihood estimator of heritability. However, since in practice the number N of SNPs is widely greater than the number of individuals n, it would be more appropriate to validate these methods in the framework where n and N go to infinity, with n/N going to a ∈ (0, +∞). Moreover, they all have been developed in a non sparse Gaussian framework, which would imply that all the available genetic information would impact the observed phenotype. This unlikely hypothesis has been discussed in particular by [START_REF] Jiang | High-dimensional genome-wide association study and misspecified mixed model analysis[END_REF], who studied the potiential error caused by non impacting SNPs in the model when considering a maximum likelihood approach from both theoretical and numerical points of view.

We have only mentioned heritability estimation in linear mixed models, but there exist other ways to define and estimate heritability. Indeed, important theoretical results on heritability estimation have been proven in the framework where n and N go to infinity, with n/N going to a ∈ (0, +∞), in the linear model

Y = Xβ + ε (1.7)
where the random component comes from the residual vector ε which is assumed to be a zeromean Gaussian vector with variance σ 2 and from the "SNP matrix" X which columns are assumed to be independent and identically distributed Gaussian variables. The heritability in this model is defined as the ratio

η = ||β|| 2 2 σ 2 + ||β|| 2 2 .
(1.8)

An advantage of this model is that there is no assumption on the distribution of β, in particular on its sparsity. However, strong assumptions are required on the stucture of the matrix X. Several methods were proposed to estimate heritability in Model (1.7). [START_REF] Dicker | Variance estimation in high-dimensional linear models[END_REF] proposed a method-of-moments estimator which is asymptotically normal when n, N → +∞ and n/N → a ∈ (0, +∞). [START_REF] Janson | Eigenprism: Inference for high-dimensional signal-tonoise ratios[END_REF] developed the Eigenprism procedure to build accurate confidence intervals for the heritability in finite sample size and also studied the asymptotic behavior of their estimator when n, N → +∞ and n/N → a ∈ (0, +∞). [START_REF] Dicker | Maximum likelihood for variance estimation in highdimensional linear models[END_REF] studied the properties of the maximum likelihood estimator and conducted a numerical comparison of the aforementioned methods which showed that the maximum likelihood estimator had a smaller empirical variance than the two others. [START_REF] Dicker | Maximum likelihood for variance estimation in highdimensional linear models[END_REF] showed the consistency and the asymptotic normality of the maximum likelihood estimator and computed as well an explicit form of the asymptotic variance.

In the same model, [START_REF] Verzelen | Adaptive estimation of high-dimensional signal-to-noise ratios[END_REF] studied the optimality of different procedures depending on the sparsity. Indeed, [START_REF] Verzelen | Adaptive estimation of high-dimensional signal-to-noise ratios[END_REF] compared the performances of an approach with variable selection (Gauss-LASSO estimator) and without selection (dense estimator) in different sparsity regimes. They computed for each range of sparsity values the minimax risk and proposed an adaptive estimator which achieves the minimax risk in all sparsity regimes.

Contribution

Our first contribution was to propose an estimator for heritability in the context where n and N go to infinity, with n/N goes to a ∈ (0, +∞) and to establish its theoretical properties. This work is developed in Chapter 2 of this manuscript and has been published in the Electronic Journal of Statistics. We studied a model as the one defined in (1.1) except that we assumed that the random effects could be sparse, that is that only a proportion q of the components of u were non-zero:

u i i.i.d.
∼ (1 -q)δ 0 + qN (0, σ u 2 ) , for all 1 ≤ i ≤ N and e ∼ N 0, σ e 2 Id R n , (1.9)

where Id R n denotes the n × n identity matrix, q is in (0, 1], and δ 0 is the point mass at 0. Up to considering the projection of Y onto the orthogonal of the image of X and for notational simplicity, we studied the following model

Y = Zu + e .
(1.10) Moreover, since in our case we are only interested in estimating η , we plugged in L n defined in (1.6) an estimator of σ 2 , that is

σ2 = 1 n n i=1 Y i 2 η(λ i -1) + 1 .
We implemented an estimator of η as the maximizer of this likelihood function depending only on parameter η:

L n (η) = -log 1 n n i=1 Y 2 i η(λ i -1) + 1 - 1 n n i=1
log (η(λ i -1) + 1) , (1.11)

We obtained two main results in the framework where n and N go to infinity, with n/N going to a ∈ (0, +∞): first, we proved that our estimator was √ n-consistent despite the presence of null components in the random effects. This result was obtained under mild assumptions on the matrix W and for any unknown sparsity q. Then we established a central limit theorem under the additional assumption that for all i and j, Z i,j were Gaussian variables with zero mean and unit variance. We computed a closed-form expression for the asymptotic variance, given by τ 2 (a, η , q) = 2 γ 2 (a, η ) + 3 a 2 η 2 γ 4 (a, η ) 1 q -1 S(a, η ) (1.12)

where

γ 2 (a, η ) = λ -1 η (λ -1) + 1 2 dµ a (λ) - λ -1 η (λ -1) + 1 dµ a (λ)
2 and S(a, η ) = λ(λ -1) (η (λ -1) + 1) 2 dµ a (λ) -λ (η (λ -1) + 1) dµ a (λ) λ -1 (η (λ -1) + 1) dµ a (λ)

2 .

In the previous expression dµ a (λ) is the density of Marchenko-Pastur, which is the distribution of the eigenvalues of ZZ /N . This distribution obtained by [START_REF] Marchenko | Distribution of eigenvalues for some sets of random matrices[END_REF] was a key element to establish the proof of our results. We implemented this approach in the R package HiLMM, which is available on the CRAN. We also conducted a simulation study with finite sample size corresponding to realistic practical studies. We showed that although the asymptotic variance defined in (1.12) was theoretically depending on the sparsity q, its influence was barely noticeable in practice. However, the asymptotic variance was shown to be very sensitive to the parameter a = n/N : more precisely, when the number of observations is very small compared to the size of the random effects (which is often the case in genetic studies), the variance of the heritability estimator increases substantially. This numerical result motivated the idea of developing a variable selection approach in order to reduce the size of the random effects and to improve the accuracy of heritability estimation.

Variable selection in the random effects of a high dimensional sparse linear mixed model

Motivated by the numerical performance of our estimator described in the previous section, it appeared to be a good idea to include a variable selection step in our method. The aim of this variable selection step is to recover the support of the random effects, which means in practice that we want to find the SNPs involved in the phenotypic variations. We would then consider only the matrix of SNPs reduced to these relevant SNPs and estimate the heritability with smaller standard error than we would have obtained with the whole matrix of SNPs. Let us first present the existing methods and results regarding variable selection in the random effects of sparse linear mixed models.

State of the art

Although the case of linear mixed models has received less attention than the linear model, there exist several methods to perform variable selection in linear mixed models.

Several works focus on selecting variable in the fixed effects of sparse LMMs, as for instance [START_REF] Schelldorfer | Estimation for high-dimensional linear mixedeffects models using 1 -penalization[END_REF]. For a complete review of these methods, we refer the reader to the work of [START_REF] Müller | Model selection in linear mixed models[END_REF]. Regarding selection in the random effects, we are only aware of the work of [START_REF] Fan | Variable selection in mixed effects models[END_REF] and [START_REF] Bondell | Joint variable selection for fixed and random effects in linear mixed-effects models[END_REF]. [START_REF] Bondell | Joint variable selection for fixed and random effects in linear mixed-effects models[END_REF] proposed indeed a method to select jointly fixed effects and random effects based on a EM algorithm. [START_REF] Fan | Variable selection in mixed effects models[END_REF] proposed a penalized criterion with a particular penalty named SCAD (Smoothly Clipped Absolute Deviation) which combines L1 and L2 penalties. Both methods can be computationally very demanding in high dimension: on the one hand, the EM algorithm, on the other hand, the cross validation to choose the two regularization parameters.

Variable selection in such high dimensional frameworks as those we are interested in can be very tricky, as proven by [START_REF] Verzelen | Minimax risks for sparse regressions: Ultra-high dimensional phenomenons[END_REF] who studied the case of the random linear model defined in Equation (1.7). [START_REF] Verzelen | Minimax risks for sparse regressions: Ultra-high dimensional phenomenons[END_REF] indeed established that if the condition N q log(1/q) >> n holds, namely when the number of causal SNPs (that is the number of non null components in the random effects) is larger than the number of individuals, the support cannot be fully recovered.

Regarding heritability estimation, the idea of introducing a variable selection step beforehand was already proposed by [START_REF] Guan | Bayesian variable selection regression for genome-wide association studies and other large-scale problems[END_REF] in a Bayesian framework. [START_REF] Guan | Bayesian variable selection regression for genome-wide association studies and other large-scale problems[END_REF] proposed indeed an approach, named BVSR (Bayesian Variable Selection Regression), that is very efficient to estimate heritability in a very sparse framework but which is biased when the number of causal SNPs is high. [START_REF] Zhou | Polygenic modeling with bayesian sparse linear mixed models[END_REF] then proposed a practical approach, called BSLMM (Bayesian Sparse Linear Mixed Model) defined as an hybrid estimator between BVSR and a classical maximum likelihood approach (without selection). This hybrid estimator behaves closely to BVSR in very sparse frameworks and like the maximum likelihood estimator (no selection) otherwise. These numerical observations of [START_REF] Zhou | Polygenic modeling with bayesian sparse linear mixed models[END_REF] are consistent with the theoretical grounds established by [START_REF] Verzelen | Adaptive estimation of high-dimensional signal-to-noise ratios[END_REF] in the linear model and described in Section 1.2.1.

Contribution

Methodology

We proposed a practical variable selection method to improve the accuracy of heritability estimation. This work has been submitted for publication and is contained in Chapter 3 of this manuscript. Our method is implemented in the R package EstHer available on the CRAN.

Our approach has two main features: firstly, it is very efficient from a statistical point of view since it provides confidence intervals considerably smaller than those obtained with methods without variable selection. Secondly, its very low computational burden makes it usable on very large data sets coming from quantitative genetics. Our method can handle ultra high dimension scenarios by using as a first step the Sure Independence Screening developed by [START_REF] Ji | UPS delivers optimal phase diagram in high-dimensional variable selection[END_REF]. Then we apply a LASSO criterion [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] combined with the stability selection [START_REF] Meinshausen | Stability selection[END_REF]. We also propose a methodology to compute confidence intervals based on a non parametric bootstrap approach and validated on synthetic data. In the course of the numerical study, we observed similar conclusions to those obtained by [START_REF] Zhou | Polygenic modeling with bayesian sparse linear mixed models[END_REF] in the Bayesian framework: in very sparse scenarios (namely, less than 200 causal SNPs out of 100 000), the estimator which includes a variable selection step is unbiased and its variance is substantially smaller than the variance of the ML estimator. However, when the number of causal SNPs is high, the selection step is not efficient and the corresponding estimator can severely underestimate the heritability. We developed a criterion based on the data in order to have an idea of the sparsity regime and whether we should apply a variable selection technique or not. We developed a hybrid estimator able to adapt according to the sparsity and we showed on synthetic data that this procedure allows us to reduce substantially the confidence intervals of the heritability estimations compared to a classical maximum likelihood estimator in very sparse scenarios. Otherwise, if the number of causal SNPs is too high, our hybrid estimator behaves like the maximum likelihood estimator, which was expected after introducing the decision criterion we proposed. The benefit of our method compared to the Bayesian approach developed by [START_REF] Zhou | Polygenic modeling with bayesian sparse linear mixed models[END_REF] lies mainly in substantially smaller computational times than for MCMC procedures, and also we do not have to deal with the settings of the different parameters in order to ensure the convergence of the algorithms.

Applications in human neuroanatomy and in animal genetics

We applied this method to two different datasets.

The first one comes from the European project IMAGEN, which is a study on teenagers' mental health. We estimated the heritability of the brain volume and the volumes of the different subcortical regions. Six phenotypes out of nine did not pass the criterion so we can suspect that a large number of SNPs are involved in their variations, and we obtained similar results to those obtained with a classical maximum likelihood approach, such as those obtained by [START_REF] Toro | Genomic architecture of human neuroanatomical diversity[END_REF] who studied the same data thanks to the software GCTA developed by [START_REF] Yang | GCTA: A tool for genome-wide complex trait analysis[END_REF]. However, for the other three phenotypes, we obtained heritability estimations with very small standard errors as well as a list of potential causal SNPs, the relevance of which could be analyzed from a biological point of view. This application to neuroanatomical data is described after the description and the validation of our method on synthetic data in Chapter 3 of this manuscript.

The second application is the study of a trout species named Salmo trutta. This brown trout, which lives in fresh water, may or may not, during its life, decide to leave fresh water to migrate to the sea. This migration has a major impact in the trout conservation, and we aim to understand the reasons of this decision. It appears that growth during the freshwater phase could potentially be a critical factor determining the fate of individuals as brown trout (remained in fresh water) or sea trout; indeed, if a fish is growing fast in fresh water then there is no real need to go to the sea where survival rate is much lower. However, if it is struggling to grow in fresh water, then the benefit of going to sea and having better growth prospects might compensate the higher predation risk. Hence we aim to investigate the proportion of genetics and environment effects in length variations, and if possible we would like to determine which SNPs and environmental variables are associated with growth patterns. The size of this data set was substantially smaller than in the previous application: we had indeed access to the length of 192 trouts, the genotype of which is described by 4069 SNPs. We noticed, according to the results of a numerical study, that our R package EstHer, which was dedicated to the analysis of very large datasets, was not efficient in this case. However, when removing the first step of our method, the Sure Independence Screening, which was specific to ultra high dimension frameworks, we obtained satisfactory results again. This application is described in Chapter 4.

Heritability estimation for binary traits

We are interested in the extension of the previous methods to the heritability estimation of a disease, where the observations are categorical (patient or control). We found in the literature different models used to define and estimate heritability for binary data.

Generalized Linear Mixed Model and Liability Model

An intuitive generalization of the previous work for estimating the heritability of a binary trait would be to consider the following Generalized Linear Mixed Model:

Y i ∼ B(q i ),
(1.13) with q i = g(l i ) where g is a link function and l i is defined as (1.14) with u ∼ N (0, σ 2 u ) and e ∼ N (0, σ 2 e ), as in the classical LMM defined in Section 1.2. A classical choice of link function in the case of binary data is for instance

l = Zu + e,
g(x) = exp(x) 1 + exp(x) ,
which ensures that q i ∈ (0, 1).

The heritability can then be defined at the liability scale, that is the heritability of the continuous variable l which is identical to the definition of heritability in the previous sections when considering a Gaussian phenotype:

η = N σ 2 u N σ 2 u + σ 2 e .
(1.15)

Another modeling and definition for heritability of a binary trait was proposed by [START_REF] Falconer | The inheritance of liability to certain diseases, estimated from the incidence among relatives[END_REF], who assumed that the binary observations could be seen as an indicator function of a Gaussian variable exceeding a certain threshold t:

Y i = 1 {l i >t} (1.16)
with l i defined by the same expression (1.14) as in the previous model.

The unobserved Gaussian variable l is also called the liability in this modeling, which is usually called the "liability model" [START_REF] Falconer | The inheritance of liability to certain diseases, estimated from the incidence among relatives[END_REF], [START_REF] Lee | Estimating missing heritability for disease from genome-wide association studies[END_REF][START_REF] Tenesa | The heritability of human disease: estimation, uses and abuses[END_REF]). The heritability is then also defined as the heritability at the liability scale as written in Equation (1.15).

Existing methods for heritability estimation in the liability model

In the literature specific to the heritability of binary phenotypes, we found methods for heritability estimation based on each of the previously described models. For the first modeling described in Equation (1.13), [START_REF] Villemereuil | Comparing parent-offspring regression with frequentist and Bayesian animal models to estimate heritability in wild populations: a simulation study for Gaussian and binary traits[END_REF] proposed to estimate the variance of the random effects σ 2 u by using MCMC methods developed by [START_REF] Hadfield | MCMC methods for multi-response generalized linear mixed models: The MCMCglmm R package[END_REF] and then to estimate the heritability as η = σ2

u σ2 u + 1 + 1
, where the first 1 in the denominator stands for the residual variance and the second 1 for the distribution-specific variance of a probit-link function [START_REF] Nakagawa | Repeatability for gaussian and non-gaussian data: a practical guide for biologists[END_REF]. The residual variance is indeed set to 1 because the binary data do not provide enough information to infer both variances σ 2 u and σ 2 e . Since the expression of the likelihood is not possible to optimize directly, [START_REF] Breslow | Approximate inference in generalized linear mixed models[END_REF] proposed to maximize a penalized quasi-likelihood, using a Laplace approximation of the likelihood. This method has been shown to underestimate the variance parameters, for instance in the numerical comparative study performed by [START_REF] Villemereuil | Comparing parent-offspring regression with frequentist and Bayesian animal models to estimate heritability in wild populations: a simulation study for Gaussian and binary traits[END_REF].

Regarding the procedures based on the second modeling defined in Equations (1.16) and (1.14), [START_REF] Lee | Estimating missing heritability for disease from genome-wide association studies[END_REF] proposed to use a maximum likelihood approach as if the binary traits were Gaussian, and then to apply a multiplicative factor to correct this approximation. [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF] showed that this heritability estimator was strongly biased in several realistic scenarios, in particular it was very sensitive to the prevalence of the disease (when the disease is rarer, the bias increases). The estimator also underestimates the heritability when the real heritability is high. [START_REF] Weissbrod | Accurate liability estimation improves power in ascertained case-control studies[END_REF] presented a different methodology to estimate heritability also in the liability model. They proposed a maximum likelihood based strategy to rebuild the underlying liability before estimating the heritability. More precisely, they differentiated the likelihood with respect to u and e and maximized the obtained functions for chosen values of σ 2 u and σ 2 e . It gave them an estimator û of u and ê of e from which they obtained an estimator of the liability as l = Zû + ê. Then they used this liability to provide another estimation of σ 2 u and σ 2 e and they repeated this procedure until convergence.

A key element of the method proposed by [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF] was to notice the particular shape of the data, which is the oversampling of cases in the case-control study. In a medical study, the number of patients is indeed similar to the number of controls even though the studied disease might be rare, which means that the proportion of cases in the study does not reflect the proportion of cases in the population. The method developed by [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF] takes into account such oversampling of the cases and, as far as we are aware, it is the only method to do so. They proposed a method of moments to estimate heritability.

More precisely, [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF] considered a simplified version of Model (1.14), where the liability is given by

l = g + e,
where g is a genetic random effect, which can be correlated across individuals, and e is the environmental random effect, which is assumed to be independent of the genetic effect. Both effects are assumed to be Gaussian: e has a variance equal to (1-η )Id R n and g has a covariance matrix, the diagonal entries of which are equal to η and the non diagonal term (i, j) is equal to η G i,j . For 1 ≤ i = j ≤ n, the covariance matrix of (l i , l j ) is given by

1 G i,j η G i,j η 1 .
They defined the variable

p i = Y i -P P (1 -P ) ,
where P is the proportion of cases in the study and the event {S = 1} happens if both individuals i and j are selected in the study. The heritability estimator proposed by [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF] is a least square estimator obtained as the minimizer of i =j

(p i p j -E(p i p j |S = 1)) 2 .
Since E(p i p j |S = 1) has no explicit formula, [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF] suggested to take advantage of the fact that the correlations G i,j were small and proposed an approximation based on Taylor developments around G i,j .

Contribution

Since the method of [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF] seemed very efficient according to their numerical results and since it is the only method we have seen which considered the specificity of the data in a case-control study, we decided to establish the theoretical properties of their estimator in the framework: n and N go to infinity, and n/N goes to a ∈ (0, +∞). We considered the model defined in Equations (1.16) and (1.14), and we assumed that Z was a random matrix with centered and normalized columns as defined in Section 1.2. In this model, the covariance matrix of (l i , l j ) can be written as

Σ (N ) = 1 + η (G N (i, i) -1) η G N (i, j) η G N (i, j) 1 + η (G N (i, i) -1),
where for all 1 ≤ i, j ≤ n,

G N (i, j) = 1 N N k=1 Z i,k Z j,k .
(1.17)

The main idea is to notice that the quantities G N (i, j), G N (i, i) -1 and G N (j, j) -1 are small, which implies that the matrix Σ (N ) is close to identity.

Inspired by the method developed by [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF], we proposed a first and second order approximation of

E[p i p j |Z, S = 1] based on Taylor developments in G N (i, j), G N (i, i) -1 and G N (j, j) -1.
Despite differences in the models we considered, we found the same first order approximation than obtained by [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF] but we have several differences in the second order approximation.

First, we investigated the theoretical properties of the estimator found with the first order approximation: we showed indeed that it was consistent under very mild assumptions on the matrix of SNPs. Then, we compared the numerical performances of the estimators obtained from the first and second order approximations, from both a statistical and computational point of view. We highlighted that the computation of the estimator obtained with the second order approximation was slower and did not improve the results compared to the first order approximation estimator. These results are contained in Chapter 5 of this manuscript.

Chapter 2

Heritability estimation in high dimensional sparse linear mixed models

The content of this chapter is contained in the article: A.

Abstract

Motivated by applications in genetic fields, we propose to estimate the heritability in highdimensional sparse linear mixed models. The heritability determines how the variance is shared between the different random components of a linear mixed model. The main novelty of our approach is to consider that the random effects can be sparse, that is may contain null components, but we do not know either their proportion or their positions. The estimator that we consider is strongly inspired by the one proposed by [START_REF] Pirinen | Efficient computation with a linear mixed model on large-scale data sets with applications to genetic studies[END_REF], and is based on a maximum likelihood approach. We also study the theoretical properties of our estimator, namely we establish that our estimator of the heritability is √ n-consistent when both the number of observations n and the number of random effects N tend to infinity under mild assumptions. We also prove that our estimator of the heritability satisfies a central limit theorem which gives as a byproduct a confidence interval for the heritability. Some Monte-Carlo experiments are also conducted in order to show the finite sample performances of our estimator.

Introduction

Linear mixed models (LMMs) have been widely used in various fields such as agriculture, biology, medicine and genetics. In quantitative genetics, LMMs have been used for estimating the heritability of traits and breeding values as explained for instance by [START_REF] Lynch | Personal genomes: The case of the missing heritability[END_REF]. In Genome Wide Association Studies (GWAS), which is the application field that inspired our work, [START_REF] Yang | GCTA: A tool for genome-wide complex trait analysis[END_REF] suggested the use of linear mixed models to measure genotypes at a large number of single nucleotide polymorphisms (SNPs) in large samples of individuals in order to identify genetic variants that explain variations in phenotypes.

The model that we shall study in this paper is a LMM defined as

Y = Xβ + Zu + e , (2.1)
where Y = (Y 1 , . . . , Y n ) is the vector of observations, X is a n × p matrix of predictors, β is a p × 1 vector containing the unknown linear effects of the predictors, and u and e correspond to the random effects. Moreover, in (2.1), Z is a n × N random matrix which will be further described in Section 2.2. We shall assume that the random effects can be sparse, that is only a proportion q of the components of u are non-zero:

u i i.i.d. ∼ (1 -q)δ 0 + qN (0, σ u 2 ) , for all 1 ≤ i ≤ N and e ∼ N 0, σ e 2 Id R n , (2.2)
where Id R n denotes the n × n identity matrix, q is in (0, 1], and δ 0 is the point mass at 0. Note that this corresponds to a more general situation than the usual assumption of (non-sparse) Gaussian random effects which is recovered when q = 1. The use of linear mixed models to estimate heritability has been proposed by [START_REF] Yang | GCTA: A tool for genome-wide complex trait analysis[END_REF] as an alternative to the regression models usually used in GWAS. The goal is to consider the joint effect of all SNPs on a phenotype, and the heritability corresponds to the proportion of phenotypic variance explained by all SNPs.

In the GWAS framework, Z is thus a matrix having a number of rows equal to the number of individuals in the experiment that is n ≈ 1000 and a number of columns equal to the number of SNPs taken into account in the experiment, namely N ≈ 500, 000. This application motivated the framework that we chose where n and N tend to infinity.

The major difference between the framework of [START_REF] Yang | GCTA: A tool for genome-wide complex trait analysis[END_REF] and ours is that they consider that the random effects are Gaussian while we consider a mixture model between a Chapter 2 -Heritability estimation in high dimensional sparse linear mixed models point mass at 0 and a Gaussian distribution. With this modeling, we assume that all SNPs are not necessarily causal, that is that all SNPs do not explain a given phenotype.

Our main goal in this paper is to propose an estimator for the heritability in this possibly sparse framework and to establish its theoretical properties in the non standard theoretical context where n and N tend to infinity.

In this paper, we prove that using a strategy close to the one proposed by [START_REF] Pirinen | Efficient computation with a linear mixed model on large-scale data sets with applications to genetic studies[END_REF], which has been devised in the case q = 1, provides consistent estimators even in the case where q < 1. Moreover, we prove that this estimator is √ n-consistent in the following asymptotic framework: n → ∞ and N → ∞ such as n/N → a > 0 and satisfies under mild assumptions a central limit theorem in both cases q = 1 and q < 1. It has to be noticed that the classical results that exist in linear mixed models are established only in the case where q = 1, n tends to infinity and N is constant.

The paper is organized as follows. Section 2.2 provides a detailed description of the model and the heritability estimator that we propose. Section 2.3 reviews existing methods for heritability estimation. Section 2.4 is dedicated to the theoretical properties of our estimator. The numerical results are presented in Section 2.5. They have been obtained thanks to the R package HiLMM that we have developed and which is available from the Comprehensive R Archive Network (CRAN). In Section 2.6, we provide some additional comments on our work as well as some prospects such as the estimation of the proportion q of non null components in the random effects. Finally, the proofs are given in Section 2.7.

Model and heritability estimator 2.2.1 Model

In the sequel, up to considering the projection of Y onto the orthogonal of the image of X and for notational simplicity, we shall focus on the following model

Y = Zu + e , (2.3)
where Y = (Y 1 , . . . , Y n ) is the vector of observations, u and e correspond to the random effects, which are defined in (2.2). Moreover, Z is a n × N random matrix such that the Z i,j are normalized random variables in the following sense: they are defined from a matrix W = (W i,j ) 1≤i≤n, 1≤j≤N by

Z i,j = W i,j -W j s j , i = 1, . . . , n, j = 1, . . . , N , (2.4)
where

W j = 1 n n i=1 W i,j , s 2 j = 1 n n i=1
(W i,j -W j ) 2 , j = 1, . . . , N .

(2.5)

In (2.4) and (2.5) the W i,j 's are such that for each j in {1, . . . , N } the (W i,j ) 1≤i≤n are independent and identically distributed random variables and such that the columns of W are independent. With this definition the columns of Z are empirically centered and normalized.

In genetic applications, the matrix W contains all the genetic information about all the individuals in the study. More precisely, for each j, the (W i,j ) 1≤i≤n are i.i.d binomial random

Chapter 2 -Heritability estimation in high dimensional sparse linear mixed models variables with parameters 2 and p j . W i,j = 0 (resp. 1, resp. 2) if the genotype of the ith individual at locus j is qq (resp. Qq, resp. QQ) where p j is the frequency of Q allele at locus j.

In Model (2.1) with (2.4), (2.5), (2.2), one can observe that Var(Y|Z) = N qσ u 2 R + σ e 2 Id R n , where R = ZZ N and q is defined in (2.2) .

Inspired by [START_REF] Pirinen | Efficient computation with a linear mixed model on large-scale data sets with applications to genetic studies[END_REF], Model (2.1) can be rewritten by using the following parameters:

σ 2 = N qσ u 2 + σ e 2 and η = N qσ u 2 N qσ u 2 + σ e 2 .
(2.6)

Thus, Var(Y|Z) = η σ 2 R + (1 -η )σ 2 Id R n .
The parameter η which belongs to [0, 1] is commonly called the heritability in the case where q = 1, see for instance [START_REF] Yang | Common snps explain a large proportion of the heritability for human height[END_REF], and determines how the variance is shared between u and e when all the components of u are non zero. We propose in (2.6) to extend this definition to the case where u may contain null components and q is in (0, 1]. The parameter q actually corresponds to the proportion of non null components in u that is to the proportion of causal SNPs. Then, the heritability defined by η in (2.6) corresponds to the proportion of phenotypic variance explained by the causal variants.

Heritability estimator

In the case where q = 1, observe that

Y|Z ∼ N 0, η σ 2 R + (1 -η )σ 2 Id R n ,
where η and σ are defined in (2.6).

Let U as the orthogonal matrix (U U = UU = Id R n ) such that URU = diag(λ 1 , . . . , λ n ) is a diagonal matrix having its diagonal entries equal to λ 1 , . . . , λ n . Hence, in the case where q = 1 and conditionally to Z, Y = U Y is a zero-mean Gaussian vector having a covariance matrix equal to diag(η σ 2 λ 1 + (1 -η )σ 2 , . . . , η σ 2 λ n + (1 -η )σ 2 ), where the λ i 's are the eigenvalues of R.

The method proposed by [START_REF] Pirinen | Efficient computation with a linear mixed model on large-scale data sets with applications to genetic studies[END_REF] consists in computing the log-likelihood

L n (σ 2 , η) = - n 2 log(σ 2 ) - 1 2 n i=1 log(η(λ i -1) + 1) - 1 2σ 2 n i=1 Y i 2 η(λ i -1) + 1 - n 2 log(2π)
and to maximize this function of two variables by iterative optimization techniques. Since in our case we are only interested in estimating η , we plugged an estimator of σ 2 that is

σ2 = 1 n n i=1 Y i 2 η(λ i -1) + 1 in L n .
Thus, in the case q = 1, the maximum likelihood strategy would lead to estimate η , assumed to be in [0, 1 -δ] with δ > 0, by η defined as a maximizer of

L n (η) = -log 1 n n i=1 Y 2 i η(λ i -1) + 1 - 1 n n i=1 log (η(λ i -1) + 1) , (2.7)
where the Y i 's are the components of the vector Y = U Y. We shall establish in Theorem 2, which is proved in Section 2.7, that this strategy produces √ n-consistent estimators of η in both cases: q = 1 and q < 1 and also that this estimator satisfies a central limit theorem which provides as a by-product confidence intervals for η .

Existing methods for heritability estimation

Several approaches can be used for estimating the heritability in the case where q = 1 but to the best of our knowledge, no theoretical results concerning the estimation of the heritability or the estimation of σ u , σ e have been established in the framework where both n and N tend to infinity. This is one of the contributions of our paper. Among these approaches, we can quote the REML (REstricted Maximum Likelihood) approach, originally proposed by [START_REF] Patterson | Recovery of inter-block information when block sizes are unequal[END_REF] and described for instance in [START_REF] Searle | Variance Components. Wiley Series in Probability and Statistics[END_REF], which consists in estimating first σ u and σ e and then to estimate η as the ratio η = N σu 2 /(N σu 2 + σe 2 ). However, this type of approach is based on iterative procedures that require expensive matrix operations. Hence, several approximations have been proposed such as the AI algorithm [START_REF] Gilmour | Average information reml: An efficient algorithm for variance parameter estimation in linear mixed models[END_REF]) which is used for instance in the software GCTA (Genome-wide Complex Trait Analysis) described in [START_REF] Yang | GCTA: A tool for genome-wide complex trait analysis[END_REF]. Other approximations have also been proposed in the EMMA algorithm [START_REF] Kang | Efficient control of population structure in model organism association mapping[END_REF]). For further details on the different approximations that could be used we refer the reader to [START_REF] Pirinen | Efficient computation with a linear mixed model on large-scale data sets with applications to genetic studies[END_REF]. The latter paper proposes another methodology for estimating the heritability which consists in rewriting Model (2.1) with the parameters (2.6) and in using an eigenvalue decomposition of the matrix R. Further details on their methodology are given hereafter. According to the numerical experiments conducted in [START_REF] Pirinen | Efficient computation with a linear mixed model on large-scale data sets with applications to genetic studies[END_REF] their approach has the lowest computational burden among the available algorithms.

In the case of sparse high dimensional framework, most of the papers studied the case of linear models. Among them, we can quote: [START_REF] Meinshausen | Stability selection[END_REF] and [START_REF] Beinrucker | Extensions of stability selection using subsamples of observations and covariates[END_REF]. The high dimensional linear mixed model where u is sparse, that is the case where q < 1, which is the most realistic case for the applications that we have in view, has received little attention. It has been studied according two directions: detection and estimation. Concerning the detection field in this framework, we are only aware of the work of [START_REF] Arias-Castro | Global testing under sparse alternatives: Anova, multiple comparisons and the higher criticism[END_REF] in which a testing procedure for detecting a sparse vector in high dimensional linear sparse regression model is also proposed and compared with the one proposed by [START_REF] Ingster | Detection boundary in sparse regression[END_REF]. As for the procedures dedicated to the heritability estimation, there exist, to the best of our knowledge, only three approaches: the approach of [START_REF] Yang | Common snps explain a large proportion of the heritability for human height[END_REF] who propose to approximate the genetic correlation between every pairs of individuals across the set of causal SNPs by the genetic correlation across the set of all SNPs, the approach of [START_REF] Golan | Accurate estimation of heritability in genome wide studies using random effects models[END_REF] who propose a methodology based on MCEM (Monte-Carlo expectation-maximization) developed by [START_REF] Wei | A monte carlo implementation of the em algorithm and the poor man's data augmentation algorithms[END_REF] and the Bayesian approaches of [START_REF] Guan | Bayesian variable selection regression for genome-wide association studies and other large-scale problems[END_REF] and [START_REF] Zhou | Polygenic modeling with bayesian sparse linear mixed models[END_REF]. However, as far as the estimation issue in the high dimensional linear mixed model is concerned, the authors of these papers did not establish the theoretical properties of their estimators in the framework where both n and N tend to infinity.

Theoretical results

Observe that another way of writing Model (2.3) with the parameters defined in (2.6) consists in writing

Y = 1 √ N Zt + σ 1 -η ε , (2.8)
where ε is a n × 1 Gaussian vector having a covariance matrix equal to identity and t = (t 1 , . . . , t N ) is a random vector such that

t i = σ √ η √ q w i π i ,
where the w i 's and the π i 's are independent, w = (w 1 , . . . , w N ) is a Gaussian vector with a covariance matrix equal to identity and the π i 's are i.i.d Bernoulli random variables such that P(π 1 = 1) = q.

The estimator η is defined as a maximizer of L n (η) for η ∈ [0, 1 -δ] for some small δ > 0, L n being given in (2.7). We shall study the asymptotic properties of η as n and N tend to infinity in a comparable way, that is when n/N → a > 0. To understand the asymptotic behavior of η, we shall first prove its consistency, then use a Taylor expansion of the derivative of L n around η in the usual way. The computations as can be seen in (2.7) involve empirical means of functions of the eigenvalues λ i of R = ZZ N . Using Theorem 1.1 of [START_REF] Bai | Large sample covariance matrices without independence structures in columns[END_REF], we shall prove the almost sure convergence of such empirical quantities under a weak assumption denoted by Assumption 1 as follows.

Assumption 1. Let Z and W be two matrices defined by (2.4) and (2.5). Recall that for each j in {1, . . . , N } the (W i,j ) 1≤i≤n are independent and identically distributed random variables and such that the columns of W are independent (but not necessarily identically distributed). Assume that the entries W i,j of W are uniformly bounded, and have variance uniformly lower bounded, that is: there exist W M < ∞ and κ > 0 such that 0 ≤ W i,j ≤ W M and σ 2 j = Var(W i,j ) ≥ κ, for all j.

The following lemma ensures that the result of [START_REF] Marchenko | Distribution of eigenvalues for some sets of random matrices[END_REF] which gives the empirical spectral distribution of sample covariance matrices ZZ /N holds even when the entries Z i,j of the matrix Z are not i.i.d. random variables but when Z is obtained by empirical standardization of a matrix W satisfying Assumption 1.

Lemma 1. Under Assumption 1, as n, N → ∞ such that n/N → a > 0, the empirical spec- tral distribution of R N = ZZ /N : F R N (x) = n -1 n
k=1 1 {λ k ≤x} tends almost surely to the Marchenko-Pastur distribution defined as the distribution function of µ a where, for any measurable set A,

µ a (A) = 1 -1 a 1 0∈A + ν a (A) if a > 1 ν a (A) if a ≤ 1 with dν a (λ) = 1 2π (a + -λ)(λ -a -) aλ 1 [a -,a + ] (x)dx, a ± = (1 ± √ a) 2 .
(2.9)

In F R N (x), the λ k 's denote the eigenvalues of R N .

Our first main result is the √ n-consistency of the estimator η.

Chapter 2 -Heritability estimation in high dimensional sparse linear mixed models

Theorem 1. Let Y = (Y 1 , . . . , Y n ) satisfy Model (2.8) with η > 0 and the entries W i,j of W satisfy Assumption 1. Then, for all q in (0, 1], as n, N → ∞ such that n/N → a ∈ (0, 1],

√ n(η -η ) = O P (1).
Such a result is a theoretical cornerstone to legitimate the use of an estimator. However, statistical inference has to be based on confidence sets. The next step is thus to find the asymptotic distribution of √ n(η -η ). Define for any η ∈ [0, 1] and λ ≥ 0

g(η, λ) = λ -1 η(λ -1) + 1 .
Define also

γ 2 n =    1 n n i=1 g(η, λ i ) 2 - 1 n n i=1 g(η, λ i ) 2    and γ 2 (a, η ) = g(η, λ) 2 dµ a (λ) - g(η, λ)dµ a (λ) 2 .
(2.10)

We are now ready to state our second main result about the asymptotic distribution of √ n(η-η ). For general q, the result only holds when the entries of Z, that is the random variables Z i,j are i.i.d. standard Gaussian. Indeed, as may be seen when computing the variances, we need to be able to find the asymptotic behavior of empirical means of functions of the eigenvalues together with the eigenvectors of the matrix R = ZZ /N .

Theorem 2. Let Y = (Y 1 , . . . , Y n ) satisfy Model (2.8) with η > 0 and assume that the random variables Z i,j are i.i.d. N (0, 1). Then for any q ∈ (0, 1], as n, N → ∞ such that n/N → a > 0, √ n(η -η )

converges in distribution to a centered Gaussian random variable with variance

τ 2 (a, η , q) = 2 γ 2 (a, η ) + 3 a 2 η 2 γ 4 (a, η ) 1 q -1 S(a, η ) where S(a, η ) = λ(λ-1) (η (λ-1)+1) 2 dµ a (λ) - λ (η (λ-1)+1) dµ a (λ) λ-1 (η (λ-1)+1) dµ a (λ) 2 .
In the case where q = 1, the result holds in the general situation described in Assumption 1, and allows us to propose confidence sets with precise asymptotic confidence level.

Theorem 3. Let Y = (Y 1 , . . . , Y n ) satisfy Model (2.8) with q = 1 and with η > 0. Assume also that the entries W i,j of W satisfy Assumption 1 then, as n, N → ∞ such that n/N → a > 0,

γ n n 2 (η -η )
converges in distribution to N (0, 1).

Let us now give some additional comments on the previous results. Firstly, it has to be noticed that none of the limiting variance depends on σ . Secondly, Theorem 2 is proved here only in the case where the Z i,j are i.i.d. Gaussian. This is because we used several times that the matrix of eigenvectors of ZZ /N is independent of the eigenvalues, and uniformly distributed on the set of orthonormal matrices. We think that the result of Theorem 2 is also valid when the Z i,j are defined from the W i,j satisfying Assumption 1, as suggested by the numerical results obtained in Section 2.5. To prove it requires new results in an active research topic of the random matrix theory field. We can observe in the expression of τ 2 (a, η ) given in Theorem 2 that the presence of q is counterbalanced by the presence of a 2 . This will be confirmed by the results obtained in the numerical results given in Section 2.5. Finally, we can observe that 2/(nγ 2 n ) corresponds to the usual inverse of the Fisher information associated to η. This result is classical in the case where N is fixed and n tends to infinity but did not exist in the framework where both n and N tend to infinity even if it was already used in biological applied papers for deriving standard errors and confidence intervals. Theorem 3 proves that this result still holds even in the case where both n and N tend to infinity.

To the best of our knowledge, the effect of the presence of null components in the random effects has never been taken into account for computing the asymptotic variance of an estimator of the heritability. This is the contribution of Theorem 2. This theorem shows that the asymptotic variance contains an additional term which increases its value in the case q < 1 with respect to the case q = 1. It is shown in Section 3.3 how the computation of the asymptotic variance can be altered if this additional term is neglected. In practical situations, computing the standard error given by Theorem 2 requires the knowledge of q which is in general unknown. However, if an estimation of q is available for any practical reasons, the result of Theorem 2 can be used for computing confidence intervals and standard errors, see Section 2.6 for further details.

Numerical experiments

In this section, we first explain how to implement our method and then we illustrate the theoretical results of Section 2.4 on finite sample size observations for both cases: q = 1 and q < 1. We also compare the results obtained with our approach to those obtained by the GCTA software described in [START_REF] Yang | Common snps explain a large proportion of the heritability for human height[END_REF] and [START_REF] Yang | GCTA: A tool for genome-wide complex trait analysis[END_REF] which is a reference in quantitative genetics.

Implementation

In order to obtain η, we used a Newton-Raphson approach which is based on the following recursion: starting from an initial value η (0) ,

η (k+1) = η (k) - L n (η (k) ) L n (η (k) ) , k ≥ 1 ,
where L n and L n denote the first and second derivatives of L n defined in (2.7), respectively. The closed form expression of these quantities are given in (2.13) and (2.25), respectively. In practice, this approach converges after at most 20 iterations and is not very sensitive to the initialization, namely to the value of η (0) . However, in particular cases, the value of the initialization can have an influence on the estimation of η . This is the case, for instance, when the real value η is close to 1. In these situations, our algorithm can provide an estimation bigger than 1 and we constrained our method to return a value equal to 0.99. Figure 2.1 shows the estimations obtained on 100 replications when a = 0.1 and η = 0.8. From this figure, we can see that the estimation of η does not depend in general on the initialization, except in some cases. Moreover, the best choice for η (0) is not constant from one replication to another. In order to limit the effect of the initialization, our algorithm uses several values for η (0) and whenever the estimations differ, it keeps the estimation which is the farthest away from the boundaries. 1.0 q qq q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q 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 2.1 -Estimation of η obtained in the case a = 0.1 and η = 0.8 for different values of initialization: η (0) = 0.1 (dots), η (0) = 0.5 (triangles) and η (0) = 0.9 (crosses). The plain line displays the estimations obtained with our method to select the best initialization value and the x-axis is the replication number.

Results in Model (2.2) when q = 1

We shall first consider the performance of the estimator η when q = 1 for η in {0.3, 0.5, 0.7}, n = 1000, σ u = 0.1 and for a in {0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1}, where a = n/N . We generated 500 data sets according to Model (2.1) using these parameters and Z as defined in (2.4) where the W i,j are binomial random variables with parameters 2 and p j . In our experiments the p j 's are uniformly drawn in [0.1, 0.5]. The corresponding boxplots of η are displayed in Figure 2.2. We can see from this figure that our approach provides unbiased estimators of η and that the smaller the a the larger the empirical variance.

In order to illustrate the central limit theorem given in Theorem 3, we first display in Figure 2.3 the histograms of γ n (n/2) 1/2 (η -η ) along with the p.d.f of a standard Gaussian random variable for η = 0.5 and different values of a. We can see that the Gaussian p.d.f fits well the data in all the considered cases. We also display in Figure 2.4 the values of n -1/2 2γ -2 n and q q q q 0.01 0.02 0.05 0.1 0.2 0.5 1 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q 0.01 0.02 0.05 0.1 0.2 0.5 1 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q q q q q q q q q q q q q q 0.01 the empirical standard deviation of (η -η ) averaged over all the experiments. As shown in Theorem 3, we also observe empirically that both quantities are very close.

In practice, the value of γ -1 n (n/2) -1/2 can be used for deriving confidence intervals for η . As we can see from Figure 2.4, our approach leads to very accurate confidence intervals for a larger than 0.1 even in finite sample size cases.

Let us now compare our results with those obtained with the software GCTA. As we can see from Figure 2.5 which displays the boxplots of η for different values of a when η = 0.7, the results found by our approach and GCTA are very close. In both cases, we observe that when a is close to 1 the estimations of η are very accurate but when a is small the standard error becomes very high.

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.1 0.2 0.3 0.4 a q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.1 0.2 0.3 0.4 a q q q q q q q Figure 2.4 -Values of n -1/2 2γ -2 n ("•"
) and the empirical standard deviation of (η -η ) (plain line) for several values of η (0.3 (left), 0.5 (right)).

q q q q q q q q q q 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.01 0.01 0.05 0.05 0.1 0.1 0.5 0.5

1 1 0.0 0.2 0.4 0.6 0.8 1.0 Figure 2
.5 -Boxplots of η for different values of a, using our method (dark gray) and GCTA (light gray). The whiskers of each boxplot are the first and third quartiles.

Results in Model (2.2) when q < 1

This section is dedicated to the study of the performance of η when q < 1. We generated 500 data sets according to Model (2.1) for η = 0.7, a ∈ {0.05, 0.1, 0.5, 1}, different values of q and Z defined in (2.4) where the W i,j are binomial random variables with parameters 2 and p j . In our experiments the p j 's are uniformly drawn in [0.1, 0.5].

Figure 2.6 displays the boxplots of η for these parameters. We can see from this figure that for small values of a, the estimators of η have the same behavior for q = 1 and q < 1. However, when a = 1 or a = 0.5, we can see from this figure that the presence of null components strongly alter the performance of the estimator of η . Since in typical GWAS experiments, a = 0.01 or even smaller, the results of Figure 2.6 could lead to conclude that considering the case q < 1 is not necessary for such values of the parameter a. However, as already noticed from Figure 2.2, the variance of η is very large for small values of a, hence considering the presence of null components and proposing a strategy for selecting only the non null components of u could be one way to increase a and thus to diminish the variance of η. q q q q q q q q q q q q q q q 0.001 0.005 0.01 0.05 0.1 0.5 1 0.0 0.2 0.4 0.6 0.8 q q q q 0.001 0.005 0.01 0.05 0.1 0.5 1 0.2 0.4 0.6 0.8 q q q q q q q q q q q q q q q q q q q 0.001 0.005 0.01 0.05 0.1 0.5 1 0.4 0.5 0.6 0.7 0.8 0.9 1.0 q q q q q q q q q q q q q q q q q q 0.001 In order to illustrate the central limit theorem given in Theorem 2, we first display in Figure 2.7 the histograms of τ -1 n n 1/2 (η -η ) along with the p.d.f of a standard Gaussian random variable for η = 0.7, two values of q: q = 0.01 and q = 0.1 and a = 0.5 (top) and two values of a: a = 0.2 and a = 0.5 with q = 0.5 (bottom). Here, τ n is the empirical version of τ (a, η , q) where γ is replaced by γ n and S(a, η ) is replaced by its empirical version with the eigenvalues of R. When a is large (a = 0.5), we can see that the higher q the better the Gaussian p.d.f fits the histograms.

Density -3 -2 -1 0 1 2 0.0 0.1 0.2 0.3 0.4 0.5 Density -3 -2 -1 0 1 2 3 0.0 0.1 0.2 0.3 0.4 Density -3 -2 -1 0 1 2 0.0 0.1 0.2 0.3 0.4 Density -2 -1 0 1 2 0.0 0.1 0.2 0.3 0.4 0.5 Figure 2.7 -Histograms of τ -1 n n 1/2 (η -η
) for a = 0.5 and q = 0.5 (top left), a = 0.1 and q = 0.1 (top right), and for a = 0.1 and q = 0.01 (bottom left), a = 0.05 and q = 0.1 (bottom right).

We also display in Figure 2.8 the values of n -1/2 τ n and the empirical standard deviation of (η -η ) averaged over all the experiments for η = 0.7 and q = 0.5. As shown in Theorem 2, we observe empirically that both quantities are very close. We also display in this figure the value of n -1/2 τ n with q = 1 which boils down to consider the asymptotic standard deviation found in the non sparse model. We can see from this figure that neglecting the term depending on q leads to underestimate the asymptotic variance of η and that this difference is all the more striking that a is close to 1. τ n with the real value of q (q = 0.5) ("•"), q = 1 (dotted line) and the empirical standard deviation of (η -η ) (plain line) for η = 0.7.

Discussion

In the course of this study, we have proposed a methodology for estimating the heritability in high dimensional linear mixed models. This methodology has two main features. Firstly, the theoretical performances of our estimator are established in a non standard theoretical framework where n and N tend to infinity and where the components of the random effect part can be equal to zero. Secondly, the computational burden of our approach is very low which makes its use possible on real data coming from GWAS experiments. As a byproduct of the central limit theorem that we establish for η we can derive a confidence interval for the heritability. However, the confidence intervals depend on q which is the proportion of non null components in u and which is general unknown. For estimating q, several strategies can be considered. One could, for instance, use a GWAS approach to compute the p-values of the correlation tests of each SNP with the observations Y and then keep only the most significant ones. Such a practical approach can be used for providing a lower bound for q. A refinement of this approach has been proposed by [START_REF] Toro | Genomic architecture of human neuroanatomical diversity[END_REF] who observed, through numerical studies, that for a fixed value of the heritability, the minimal p-value is all the more low that the number of causal SNPs is small. Hence, performing a GWAS approach on a given data set allows them to have an idea of the number of SNPs which are likely to be causal. One could also propose another practical method based on a variable selection technique. Such an approach has already been proposed by [START_REF] Fan | Variable selection in mixed effects models[END_REF] in the context of sparse linear mixed models. However, the framework in which their theoretical results are derived is different from the one that is considered in our paper. We are currently working on a paper [START_REF] Bonnet | Improving heritability by a variable selection approach in sparse high dimensional linear mixed models[END_REF] which presents a variable selection method which is adapted to our framework and which could be used for estimating the proportion q of non null components in the random effects.

Moreover, we did not take into account the linkage disequilibrium issue which would require to extend our results to the case where the columns of the random matrix are correlated. This question will be the subject of a future work.

where U (already introduced in Section 2.1) is a n × n orthonormal matrix, V is a N × N orthonormal matrix and √ D is a n × n diagonal matrix having its diagonal entries equal to √ λ i , the λ i 's being the eigenvalues of R = ZZ /N previously defined. Thus, (2.8) rewrites as

Y = U Y = √ D 0 V t + σ 1 -η ε , (2.11)
where ε = U ε is a n × 1 centered Gaussian vector having a covariance matrix equal to identity. We shall use repeatedly the following lemma which is proved in Section 2.7.4.

Lemma 2. Let Y be defined by (2.11) and H be a n × n diagonal matrix, then

Var Y H Y|Z = 2σ 4 Tr H 2 {(1 -η )Id R n + η D} 2 + 3σ 4 η 2 1 q -1 1≤i≤N M 2 ii ,
where

M = V DH 0 0 0 V , and
Var Y H Y|Z ≤ 2σ 4 Tr H 2 {(1 -η )Id R n + η D} 2 + 3σ 4 η 2 1 q -1 Tr[D 2 H 2 ].
Another useful lemma will be the following.

Lemma 3. Under Assumption 1, let h : R + → R + be such that there exist α > 0 and C such that for all n,

E 1 n n i=1 h(λ i ) 1+α ≤ C. Then 1 n n i=1 h(λ i ) = h(λ)dν a (λ) + o p (1).
The proof of this lemma follows from the application of Lemma 1 to the bounded function h1 h≤M , and the Markov inequality applied to the empirical mean of h1 h>M .

Lemma 4. Under Assumption 1 let n, N → ∞ be such that n/N → a > 0. Then there exists C such that for all n,

E 1 n n i=1 λ 2 i ≤ C.
To prove the lemma, notice that n i=1 λ

2 i = Tr[ZZ /N 2 ]. But E Tr (ZZ ) 2 = k =k i,j E(Z i,k Z j,k )E(Z i,k Z j,k ) + k i E(Z 2 i,k ) = nN (N -1) + N (N -1)n(n -1) 1 n -1 2 + n 2 N
where the values of the involved expectations may be found in the proof of Lemma 1 in Section 2.7.4. We thus have

E 1 n n i=1 λ 2 i ≤ 2 + n N
which ends the proof.

Proof of Theorem 1

The first step is to prove the consistency of η. We shall first prove that L n (η) converges uniformly for η ∈ [0, 1 -δ] in probability to L(η) given by

L(η) = -2 log σ -log η (λ -1) + 1 η(λ -1) + 1 dµ a (λ) -log (η(λ -1) + 1) dµ a (λ).
Using Lemma 2 with H with diagonal entries 1/(η(λ i -1) + 1), we get that

Var 1 n n i=1 Y 2 i η(λ i -1) + 1 |Z ≤ σ 4 n 2 n i=1 2 η (λ i -1) + 1 η(λ i -1) + 1 2 + 3 1 q -1 η λ i η(λ i -1) + 1 2 ≤ σ 4 2 + 3 1 q -1 1 n 2 n i=1 λ i + 1 δ 2 since η ∈ [0, 1 -δ]. Now, using Lemma 4 we get that 1 n 2 n i=1 λ i + 1 δ 2 = o P (1) which leads to 1 n n i=1 Y 2 i η(λ i -1) + 1 = E 1 n n i=1 Y 2 i η(λ i -1) + 1 |Z + o p (1) = σ 2 1 n n i=1 η (λ i -1) + 1 η(λ i -1) + 1 + o P (1).
Now, using Lemma 3 we easily get that

1 n n i=1 η (λ i -1)+1 η(λ i -1)+1 converges in probability to [ η (λ-1)+1 η(λ-1)+1 ]dµ a (λ) and 1 n n i=1 log[(η(λ i -1) + 1)] converges in probability to log(η(λ -1) + 1)dµ a (λ) so that L n (η) = L(η) + o P (1).
In order to prove the uniform convergence of L n to L in probability on [0, 1 -δ], we shall use the following lemma which is proved in section 2.7.4.

Chapter 2 -Heritability estimation in high dimensional sparse linear mixed models Lemma 5. Assume that for any η ∈ [0, 1 -δ], L n (η) converges in probability to L(η) and that sup η∈[0,1-δ] L n (η) = O P (1), as n tends to infinity, (2.12) then sup

η∈[0,1-δ]
|L n (η) -L(η)| = o P (1), as n tends to infinity.

Let us now prove that sup

η∈[0,1-δ] |L n (η)| = O P (1). Note that L n (η) = 1 n n i=1 Y 2 i (λ i -1) {η(λ i -1) + 1} 2 1 n n i=1 Y 2 i η(λ i -1) + 1 -1 - 1 n n i=1 λ i -1 η(λ i -1) + 1 . (2.13) A study of η → 1 n n i=1 Y 2 i (λ i -1) {η(λ i -1)+1} 2 1 n n i=1 Y 2 i η(λ i -1)+1 -1
shows that it is decreasing and that it takes negative values for η ∈ [0, 1-δ], so that its absolute value is maximum for η = 1-δ. Thus

sup η∈[0,1-δ] L n (η) ≤ 1 δ 1 n n i=1 Y 2 i |λ i -1| 1 n n i=1 Y 2 i -1 + 1 nδ n i=1 |λ i -1| ≤ 2 δ + 1 δ 1 n n i=1 Y 2 i λ i 1 n n i=1 Y 2 i -1 + 1 nδ n i=1 λ i .
By Lemma 2 with H = Id, we get

1 n n i=1 Y 2 i = E 1 n n i=1 Y 2 i |Z + o p (1) = σ 2 n n i=1 [η (λ i -1) + 1)] + o p (1) = σ 2 (η(λ -1) + 1)dµ a (λ) + o p (1),
where the last equality comes from Lemma 3. In the same way, we get by using Lemma 2 with H having its diagonal entries equal to λ i and Lemma 3 that

1 n n i=1 Y 2 i λ i = σ 2 λ(η(λ -1) + 1)dµ a (λ) + o p (1) = O P (1).
Finally, we get from Lemma 3 that

1 n n i=1 λ i = λdµ a (λ) + o p (1) = O P (1)
which ends the proof of (2.12). By Lemma 5, we thus have proved that

sup η∈[0,1-δ] |L n (η) -L(η)| = o P (1). (2.14)
Now, using Jensen's inequality, we easily get that for all η ∈ [0, 1], L(η) ≤ L(η ), with equality if and only if η = η . This together with (2.14) gives η = η + o P (1).

(2.15)

The next step is to prove that √ n(η -η ) = O P (1). Let us first note that η satisfies the following equation:

√ n(η -η ) = - √ nL n (η ) L n ( η) , η ∈ (η, η ) . (2.16)
We first prove the asymptotic convergence of L n ( η).

Lemma 6. Let Y = (Y 1 , . . . , Y n ) satisfy Model (2.8) with η > 0 and the entries W i,j of W satisfy Assumption 1. Then, for all q in (0, 1], as n, N → ∞ such that n/N → a ∈ (0, 1], for any random variable η such that η ∈ (η, η ),

L n ( η) = -σ 2 γ 2 (a, η ) + o P (1).
Lemma 6 is proved in Section 2.7.4. Let us now focus on the properties of L n (η ). Using the following notation

U i = Y i η (λ i -1) + 1 , (2.17)
we see that √ nL n (η ) can be rewritten as follows:

   1 √ n n i=1   U i 2 - 1 n n j=1 U 2 j   g(η , λ i )    1 n n i=1 U 2 i -1 =    1 √ n n i=1   U i 2 -1 +   1 - 1 n n j=1 U 2 j     g(η , λ i )    1 n n i=1 U 2 i -1 = 1 √ n n i=1 U i 2 -1 g(η , λ i ) 1 n n i=1 U 2 i -1 -    1 √ n n j=1 U j 2 -1    1 n n i=1 g(η , λ i ) 1 n n i=1 U 2 i -1 , where g(η, λ) = λ -1 η(λ -1) + 1 .
But using Lemma 2 and Lemma 3 we get

Var   n -1/2 n j=1 (U j 2 -1)|Z   = O P (1)
Moreover, by Lemma 3, n -1 n i=1 g(η , λ i ) converges in probability to g(η , λ)dµ a (λ). Thus,

√ nL n (η ) = 1 √ n n i=1 U i 2 -1 g(η , λ i ) -g(η , λ)dµ a (λ) + o P (1), as n → ∞ . (2.18)
Chapter 2 -Heritability estimation in high dimensional sparse linear mixed models

Using again Lemma 2 and Lemma 3 we obtain

√ nL n (η ) = O P (1).
This, together with Lemma 6 and (2.16) ends the proof of Theorem 1.

Proof of Theorem 2

Notice first that all previous results may be used, replacing Assumption 1 by the assumption that the Z i,j are i.i.d. standard Gaussian. Indeed, in this case, Lemma 1 reduces to the original result of [START_REF] Marchenko | Distribution of eigenvalues for some sets of random matrices[END_REF], Lemma 3 only involves Lemma 1 and truncation arguments, and the computations leading to Lemma 4 still hold. Thus, Theorem 1 and Lemma 6 also still hold.

Let us now prove that √ nL n (η ) converges in distribution to a centered Gaussian. Define H the diagonal n × n matrix with diagonal entries

H i = 1 η (λ i -1) + 1 g(η , λ i ) -g(η , λ)dµ a (λ) . Define L n = 1 √ n Y H Y.
Then using (2.18) and Lemma 3 we have

√ nL n (η ) = L n -E[L n |Z] + o P (1).
Now using Lemma 2 we get that setting

γ 2 n = Var [L n |Z], γ 2 n = 2σ 4 1 n Tr H 2 ((1 -η) Id R n + η D) 2 + 3σ 4 η 2 1 q -1 1 n N i=1 M 2 i,i = 2σ 4 1 n n i=1 g(η , λ i ) -g(η , λ)dµ a (λ) 2 + 3σ 4 η 2 1 q -1 1 n n i=1 n k,l=1 λ k λ l H k H l V 2 i,k V 2 i,l .
The first term in this sum converges as n, N → ∞ to 2σ 4 γ 2 (a, η ). Under the assumption that the Z i,j are i.i.d. standard Gaussian, the matrix of eigenvectors V is Haar distributed on the orthonormal matrices, and is independent of (λ i ) 1≤i≤n , see [START_REF] Bai | Spectral analysis of large dimensional random matrices[END_REF] chapter 10. Conditionally to the eigenvalues (λ i ) 1≤i≤n , we thus get that

E   1 n n i=1 n k,l=1 λ k λ l H k H l V 2 i,k V 2 i,l |D   = 1 N n k=1 λ k H k 2 (1 + o(1)) converges to a 2 λ(λ -1) (η (λ -1) + 1) 2 dµ a (λ) - λ (η (λ -1) + 1) dµ a (λ) λ -1 (η (λ -1) + 1) dµ a (λ)
Chapter 2 -Heritability estimation in high dimensional sparse linear mixed models and

Var   1 n n i=1 n k,l=1 λ k λ l H k H l V 2 i,k V 2 i,l |D   = o P (1) so that γ 2 n = 2σ 4 γ 2 (a, η ) + 3σ 4 η 2 1 q -1 S(a, η ) + o P (1). Denote ∆ the diagonal N × N -matrix with diagonal entries ∆ i = σ √ η √ q π i . Let us now write L n -E(L n |Z) = L n -E [L n |∆, Z] + E [L n |∆, Z] -E [L n |Z] .
We first have

E [L n |∆, Z] -E [L n |Z] = σ 2 η 1 √ n N i=1 π 2 i q -1 M i,i
whose variance, conditionally to Z is

s 2 n,1 = σ 4 η 2 1 q -1 1 n N i=1 M 2 i,i .
In the same way as for γ 2 n we get that

s 2 n,1 = σ 4 η 2 1 q -1 S(a, η ) + o P (1).
Let

ξ i = π 2 i q -1 M i,i = π 2 i q -1 n k=1 λ k (λ k -1) (η (λ k -1) + 1) 2 V 2 i,k . Since η > 0, the function λ → λ(λ-1) (η (λ-1)+1) 2 is bounded, and n k=1 V 2 i,k ≤ N k=1 V 2 i,k = 1. Also, the variables π 2 i q -1 are uniformly bounded by 1/q. Thus 1 n n i=1 E ξ 2 i 1 |ξ i |≥cn |Z = 0
for large enough n. Then, by Lindeberg's Theorem, conditionally to Z,

1 s n,1 (E [L n |∆, Z] -E [L n |Z])
converges in distribution to N (0, 1).

Let us now set

s 2 n,2 = γ 2 n -s 2 n,1
and notice that s 2 n,2 converges to

2σ 4 γ 2 (a, η ) + 2σ 4 η 2 1 q -1 S(a, η ).
We shall prove that, conditionally to Z and ∆, (L n -E(L n |∆, Z))/s n,2 converges in distribution to N (0, 1), and thus also unconditionally. Write

L n = 1 √ n (w ε ) B w ε where B is the (N + n) × (N + n)-matrix B = ∆ 0 0 σ (1 -η ) 1 2 Id R n     V DH 0 0 0 V Ṽ√ DH H √ D Ṽ H     ∆ 0 0 σ (1 -η ) 1 2 Id R n .
Here, Ṽ is the N ×n-matrix which consists of the first n columns of V. Let φ be the characteristic function of (L n -E(L n |∆, Z))/s n,2 conditionally to Z and ∆. Notice first that if b j , j = 1, . . . , n + N are the eigenvalues of B, we may write

L n -E [L n |∆, Z] = 1 √ n N +n j=1 b j (e 2 j -1).
for random variables e j i.i.d. standard Gaussian. Thus

φ (t) = N +n j=1 1 -2i tb j s n,2 √ n -1/2 exp -i tb j s n,2
√ n and Taylor expansion leads to

log φ (t) = N +n j=1 - 1 2 log 1 -2i tb j s n,2 √ n -i tb j s n,2 √ n = -t 2 1 ns 2 n,2 N +n j=1 b 2 j + O   1 n √ ns 3 n,2 N +n j=1 b 3 j   .
We shall now prove that 1

ns 2 n,2 N +n j=1 b 2 j converges to 1/2. Tedious computations give N +n j=1 b 2 j = Tr(B 2 ) = Tr(∆M ∆ 2 M ∆) + σ 4 (1 -η 2 Tr(H 2 ) + 2σ 2 (1 -η ) Tr[∆ 2 ṼDH 2 Ṽ ].
Using the distribution of V and its independence on D we get

E   N +n j=1 b 2 j |D   = 2σ 4 Tr H 2 ((1 -η) Id R n + η D) 2 +2σ 4 η 2 1 q -1 1 N n k=1 λ k H k 2 (1 + o(1))
Chapter 2 -Heritability estimation in high dimensional sparse linear mixed models so that

E   1 n N +n j=1 b 2 j |D   = 2σ 4 γ 2 (a, η ) + 2σ 4 η 2 1 q -1 S(a, η ) + o P (1).
Moreover, tedious computations again give

Var   1 n N +n j=1 b 2 j |D   = o P (1),
and we obtain that

1 ns 2 n,2 N +n j=1 b 2 j = 1 2 + o P (1).
We shall now prove that

1 n √ ns 3 n,2 N +n j=1 b 3 j = o P (1). To do so, it is enough to prove that max j |b j | = o P ( √ n). Notice that for any normed vector A = (A 1 , A 2 ) in R N +n where A 1 ∈ R N and A 2 ∈ R n , max j |b j | ≤ A BA. Now, A BA = A 1 (∆M ∆)A 1 + 2σ 1 -η A 1 (∆ Ṽ√ DH)A 2 + σ 2 (1 -η )A 2 HA 2 .
First, since η > 0, all entries of H and D and HD are uniformly bounded and so are all entries of ∆. We thus get A 2 HA 2 = O(1) and A 1 (∆ Ṽ√ DH)A 2 = O(1). Then, using the distribution of V and its independence on D we get

E A 1 (∆M ∆)A 1 |D = O 1 N n i=1 λ i H i and Var A 1 (∆M ∆)A 1 |D = o P (1), so that A BA = O P (1). We have thus proved that max j |b j | = O P (1) = o P ( √ n).
Thus φ(t) converges in probability for all t to exp -t 2 2 and the convergence may be strengthened by contradiction to an a.s. convergence, so that conditionally to Z and ∆, (

L n -E(L n |∆, Z))/s n,2 converges in distribution to N (0, 1). Now, conditionally to Z and ∆, (L n -E(L n |∆, Z))/s n,2 converges in distribution to a Gaussian random variable independent of ∆. Thus conditionally to Z, L n -E [L n |∆, Z] and E [L n |∆, Z] -E [L n |Z]
converge in distribution to independent Gaussian variables, so that their sum converges in distribution to a centered Gaussian with variance the sum of the variances, namely the limit of γ 2 n , and Theorem 2 is proved.

Proof of Theorem 3

Using Lemma 6 and (2.16) , there remains to prove that √ nL n (η ) converges in distribution to N (0, 2σ 4 γ 2 (a, η )) and that γ 2 n converges in probability to γ 2 (a, η ). Notice first that when q = 1, (U 1 , . . . , U n )|Z is a centered Gaussian vector with a covariance matrix equal to σ 2 times the identity matrix. We shall prove that conditionally to Z, √ nL n (η ) converges in distribution to N (0, 2σ 4 γ 2 (a, η )) so that the result still holds unconditionally. Using (2.18), it is only needed to prove it for 1

√ n n i=1 U i 2 -1 g(η , λ i ) -g(η , λ)dµ a (λ) .
Now, conditionally to Z, the variance of

n i=1 U i 2 -1 g(η , λ i ) -g(η , λ)dµ a (λ) is γ 2 n = 2σ 4 n n i=1 g(η , λ i ) -g(η , λ)dµ a (λ) 2 .
Since η > 0, g(η , λ) is a bounded function of λ, and using Lemma 3,

γ 2 n = 2σ 4 γ 2 (a, η )) + o P (1). Also, setting ξ i = U i 2 -1 g(η , λ i ) -g(η , λ)dµ a (λ) and C an upper bound of |g(η , λ)|, we get that for any c > 0, 1 n n i=1 E ξ 2 i 1 |ξ i |≥cn |Z ≤ 4C 2 σ 4 E U 1 2 -1 2 1 2C|U 1 2 -1|≥cn |Z = 4C 2 σ 4 E U 1 2 -1 2 1 2C|U 1 2 -1|≥cn = o(1),
where the first equality comes from the fact that the distribution of (U 1 , . . . , U n )|Z does not depend on Z and is thus also the distribution of (U 1 , . . . , U n ). Then, using Lindeberg's Theorem, conditionally to Z, √ nL n (η ) converges in distribution to N (0, 2σ 4 γ 2 (a, η )) and thus also unconditionally. The fact that γ 2 n converges in probability to γ 2 (a, η ) is a straightforward consequence of Taylor expansion, the fact that g(η , λ) and its derivative with respect to η in the neighborhood of η are bounded functions of λ, and Slutzky's Lemma.

Proofs of technical lemmas

Proof of Lemma 1

As a byproduct of Theorem 1.1, Corollary 1.1 and Remark 1.1 of [START_REF] Bai | Large sample covariance matrices without independence structures in columns[END_REF], we use the following result to prove Lemma 1.

Theorem [START_REF] Bai | Large sample covariance matrices without independence structures in columns[END_REF]). Let Z be a matrix of size n × N which columns, denoted by Z 1 , . . . , Z N , are independent and let us denote Z = 1

N N k=1 Z k . Let us also recall that R = ZZ /N and F R is its empirical spectral distribution defined by F R (x) = 1 n n k=1 1 {λ k ≥x} ,
where λ 1 , . . . , λ n are the eigenvalues of R. As N → ∞, assume the following:

Chapter 2 -Heritability estimation in high dimensional sparse linear mixed models

2. 1 N max i =m E( Zi,j Z m,j ) 2 → 0 uniformly in j ≤ N . 3. 1 N 2 Λ E( Zi,j Z m,j -t m,i )(Z i ,j Zm ,j -t i ,m ) 2 → 0 uniformly in j ≤ N , with Λ = {(i, m, i , m ) : 1 ≤ i, m, i , m ≤ n}\{(i, m, i , m ) : i = i = m = m or i = m = i = m}. 4. n N → a ∈ (0, +∞).
5. The norm of T is uniformly bounded and F T tends to a degenerate distribution with mass at 1/a.

Then, with probability 1, F R converges to the Marchenko-Pastur distribution defined in (2.9).

Observe that for all j = 1, . . . , N ,

n i=1 Z i,j = 0 (2.19) and n i=1 Z 2 i,j = n.
(2.20)

Moreover, for each j, the random variables (Z i,j ) 1≤i≤n are exchangeable. Thus, we deduce from (2.20) that for all i = 1, . . . , n and j = 1, . . . , N , E(Z 2 i,j ) = 1. Hence, by (2.19), we get that

0 = n i=1 Z i,j 2 = n i=1 Z 2 i,j + 1≤i =m≤n Z i,j Z m,j ,
which, by (2.20), implies that for all j = 1, . . . , N and i = m = 1, . . . , n,

E(Z i,j Z m,j ) = - n n(n -1) = - 1 n -1 .
(2.21) Thus, the matrix T = T n defined in Theorem [START_REF] Bai | Large sample covariance matrices without independence structures in columns[END_REF]) is equal to

T = n/(n - 1)Id R n -J n /(n -1)
, where J n is a n × n matrix having all its entries equal to 1. Hence the eigenvalues of T are 0 with multiplicity 1 and n/(n -1) with multiplicity (n -1), which gives Condition 5. of Theorem [START_REF] Bai | Large sample covariance matrices without independence structures in columns[END_REF]).

Let us then check Condition 2. of Theorem [START_REF] Bai | Large sample covariance matrices without independence structures in columns[END_REF]). Observe that,

for i = m, E[(Z i,j Z m,j -t m,i ) 2 ] = E(Z 2 i,j Z 2 m,j ) -t 2 m,i
. By (2.20), for all j = 1, . . . , N ,

n 2 = n i=1 Z 2 i,j 2 = n i=1 Z 4 i,j + 1≤i =m≤n Z 2 i,j Z 2 m,j .
Since the (Z i,j ) 1≤i≤n are exchangeable for each j = 1, . . . , N , we get that for all j = 1, . . . , N ,

n = E[Z 4 1,j ] + (n -1)E[Z 2 1,j Z 2 2,j ] .
Thus, for all j = 1, . . . , N , E[Z 2 1,j Z 2 2,j ] ≤ n/(n -1), which with the definition of the t m,i 's gives the result.

Let us now check Condition 3. of Theorem [START_REF] Bai | Large sample covariance matrices without independence structures in columns[END_REF]). Since the random variables (Z i,j ) 1≤i≤n are exchangeable, it is enough to prove that, uniformly in k, Chapter 2 -Heritability estimation in high dimensional sparse linear mixed models

(i) E[Z 4 1,k ] = o( √ n), (ii) E[Z 2 1,k Z 2 2,k ] -1 = o(1), (iii) E[Z 3 1,k Z 2,k ] = o(1), (iv) √ nE[Z 2 1,k Z 2,k Z 3,k ] = o(1), (v) nE[Z 1,k Z 2,k Z 3,k Z 4,k ] = o(1) , as n → ∞.
Observe that (i) implies (ii). Using (2.19), by expanding 0 = (

n i=1 Z i,k ) 2 n i=1 Z 2 i,k
and taking the expectation, we get that (i) and (iii) imply (iv). By expanding 0 = ( n i=1 Z i,k ) 4 , which comes from (2.19), and by taking the expectation, (i) and (iii) imply (v). Hence, it is enough to prove (i) and (iii) to conclude the proof of Lemma 1.

Let us first prove (i). By the definition of Z 1,k given in (2.4), we get that for all k, Z 2 1,k ≤ n. Hence,

Z 2 1,k ≤ (W 1,k -W k ) 2 2σ 2 k 1 {s 2 k ≥σ 2 k /2} + n1 {s 2 k <σ 2 k /2} ,
and, by the assumptions on the W i,k 's and on the σ k 's,

E(Z 4 1,k ) ≤ W 2 M κ 2 /2 + n 2 P(σ 2 k -s 2 k > σ 2 k /2) .
Theorem A of (Serfling, 1980, p. 201) implies that the second term of the previous inequality tends to zero as n tends to infinity uniformly in k, which concludes the proof of (i).

Let us now prove (iii). Using (2.19), we get Z 3 1,k ( n i=1 Z i,k ) = 0. By expanding this equation and taking the expectation, we obtain that

E(Z 4 1,k )+ n i=2 E(Z 3 1,k Z i,k ) = 0. Since the (Z i,k ) 1≤i≤n are exchangeable: E(Z 3 1,k Z 2,k ) = -E(Z 4 1,k )/(n -1) = o(n -1/2
), where the last equality comes from (i).

Proof of Lemma 2

Using (2.11) and the independence assumptions, we get

Var( Y H Y|Z) = Var v V DH 0 0 0 V v + 2σ 1 -η v V √ D 0 H ε + σ 2 (1 -η ) ε H ε|Z = Var v Mv|Z + 4σ 2 (1 -η ) Var v V √ D 0 H ε|Z + 2σ 4 (1 -η ) 2 Tr(H 2 ) , (2.22)
where M = V DH 0 0 0 V . Using the independence assumptions, we get that

4σ 2 (1 -η ) Var v V √ D 0 H ε|Z = 4σ 4 η (1 -η ) Tr(BB ) = 4σ 4 η (1 -η ) Tr(DH 2 ) , (2.23)
Chapter 2 -Heritability estimation in high dimensional sparse linear mixed models

where B = V √ D 0 H. Moreover, E(v Mv|Z) = σ 2 η Tr(D 2 H 2 ) and E (v Mv) 2 |Z = σ 4 η 2 q 2   2q 2 1≤i =j≤N M 2 ij + q 2 1≤i =i ≤N M ii M i i + 3q 1≤i≤N M 2 ii   = σ 4 η 2   2 Tr(M 2 ) -2 1≤i≤N M 2 ii + Tr(M) 2 - 1≤i≤N M 2 ii + 3 q 1≤i≤N M 2 ii   = σ 4 η 2   2 Tr(D 2 H 2 ) + Tr(M) 2 + 3 1 q -1 1≤i≤N M 2 ii   .
Thus,

Var v Mv Z] = σ 4 η 2   2 Tr(D 2 H 2 ) + 3 1 q -1 1≤i≤N M 2 ii   .
(2.24)

The proof of the equality in Lemma 2 follows from (2.22), (2.23) and (2.24). The proof of the inequality in Lemma 2 follows now from

1≤i≤N M 2 ii ≤ 1≤i,j≤N M 2 ij = Tr[D 2 H 2 ].

Proof of Lemma 5

Let > 0 and let

{η 1 < • • • < η K( ) } be a grid of [0, 1 -δ] such that |η j -η j+1 | < for all j ∈ {0, . . . , K } then sup η∈[0,1-δ] |L n (η) -L(η)| ≤ sup j∈{1,...,K } sup η ∈[η j ,η j+1 ] |L n (η ) -L n (η j )| + |L n (η j ) -L(η j )| + sup η ∈[η j ,η j+1 ] |L(η j ) -L(η )| ≤ sup η∈[0,1-δ] |L n (η)| + sup j∈{1,...,K } |L n (η j ) -L(η j )| + ω( ),
where ω( ) is the modulus of continuity of L, which is continuous on the compact [0, 1 -δ] and thus uniformly continuous on this compact. Since sup

η∈[0,1-δ]
|L n (η)| = O P (1) then, for every β > 0, there exists C such that for all n, P( sup

η∈[0,1-δ] |L n (η)| ≥ C) ≤ β.
Let α > 0 and let us consider the -grid such that ≤ α/3C and ω( ) ≤ α/3, thus we get that

P( sup η∈[0,1-δ] |L n (η) -L(η)| ≥ α) ≤ P( sup η∈[0,1-δ] L n (η) ≥ C) + P( sup j∈{1,...,K } |L n (η j ) -L(η j )| ≥ α -C -ω( )) ≤ P( sup η∈[0,1-δ] L n (η) ≥ C) + P( sup j∈{1,...,K } |L n (η j ) -L(η j )| ≥ α 3 ) ≤ P( sup η∈[0,1-δ] L n (η) ≥ C) + K j=1 P(|L n (η j ) -L(η j )| ≥ α 3 ),
which concludes the proof of Lemma 5 since each term tends to zero as n tends to infinity.

Proof of Lemma 6

The second derivative of L n is given by

L n (η) = - 2 n n i=1 Y 2 i (λ i -1) 2 {η(λ i -1) + 1} 3 1 n n i=1 Y 2 i {η(λ i -1) + 1} -1 + 1 n n i=1 Y 2 i (λ i -1) {η(λ i -1) + 1} 2 2 1 n n i=1 Y 2 i {η(λ i -1) + 1} -2 (2.25) + 1 n n i=1 (λ i -1) 2 {η(λ i -1) + 1} 2 .
In particular for η = η , we have

1 n n i=1 Y 2 i {η (λ i -1) + 1} = 1 + o P (1),
and using as previously Lemma 2, Lemma 3 and the fact that all functions of λ involved in the empirical means are bounded since η > 0, we get

2 n n i=1 Y 2 i (λ i -1) 2 {η(λ i -1) + 1} 3 = 2σ 2 n n i=1 (λ i -1) 2 {η(λ i -1) + 1} 2 + o P (1) = 2σ 2 (λ -1) 2 {η(λ -1) + 1} 2 dµ a (λ) + o P (1) and 1 n n i=1 Y 2 i (λ i -1) {η(λ i -1) + 1} 2 = σ 2 n n i=1 (λ i -1) {η(λ i -1) + 1} + o P (1) = σ 2 (λ -1) {η(λ -1) + 1} dµ a (λ) + o P (1) leading to L n (η) = -σ 2 γ 2 (a, η ) + o P (1).
Using Slutzky's Lemma and η = η + o P (1), there just remains to prove that for small enough α > 0, sup

|η-η |≤α |L n (η) -L n (η)| = O p (α).
But this comes easily from

sup |η-η |≤α |L n (η) -L n (η)| ≤ α sup |η-η | |L (3) n (η)|
where L

(3) n (η) is the third derivative of L n (η), and a similar handling of empirical means as before. Indeed, all functions of λ involved are bounded as soon as α is such that η ≥ 2α.

Chapter 3

Improving heritability estimation by a variable selection approach in high dimensional sparse linear mixed models

The content of this chapter is contained in the article sumbitted for publication: A. Bonnet, C. Lévy-Leduc, E. Gassiat, R. Toro, and T. Bourgeron. Improving heritability by a variable selection approach in sparse high dimensional linear mixed models, 2016, http://arxiv.org/abs/1507.06245v3. The method which is presented is implemented in the EstHer R package, available on the CRAN. Chapter 3 -Improving heritability estimation by a variable selection approach in high dimensional sparse linear mixed models

Abstract

Motivated by applications in neuroanatomy, we propose a novel methodology for estimating the heritability which corresponds to the proportion of phenotypic variance which can be explained by genetic factors. Estimating this quantity for neuroanatomical features is a fundamental challenge in psychiatric disease research. Since the phenotypic variations may only be due to a small fraction of the available genetic information, we propose an estimator of the heritability that can be used in high dimensional sparse linear mixed models. Our method consists of three steps. Firstly, a variable selection stage is performed in order to recover the support of the genetic effects -also called causal variants -that is to find the genetic effects which really explain the phenotypic variations. Secondly, we propose a maximum likelihood strategy for estimating the heritability which only takes into account the causal genetic effects found in the first step. Thirdly, we compute the standard error and the 95% confidence interval associated to our heritability estimator thanks to a nonparametric bootstrap approach. Our contribution consists in providing an estimation of the heritability with standard errors substantially smaller than methods without variable selection when the genetic effects are very sparse. Since the real genetic architecture is in general unknown in practice, we also propose an empirical criterion which allows the user to decide whether it is relevant to apply a variable selection based approach or not. We illustrate the performance of our methodology, implemented in the R package EstHer, on synthetic and real neuroanatomic data coming from the Imagen project. We also show that our approach has a very low computational burden and is very efficient from a statistical point of view.

Introduction

For many complex traits in human population, there exists a huge gap between the genetic variance explained by population studies and the variance explained by specific variants found thanks to genome wide association studies (GWAS). This gap has been called by [START_REF] Lynch | Personal genomes: The case of the missing heritability[END_REF] and [START_REF] Manolio | Finding the missing heritability of complex diseases[END_REF] the "dark matter" of the genome or the "dark matter" of heritability. Various population studies have shown that up to 80% of the variability of neuroanatomical phenotypes such as the brain volume could be explained by genetic factors, see for instance Stein et al. (2012). This result is very important since several psychiatric disorders are shown to be associated to neuroanatomical changes, for instance macrocephaly and autism [START_REF] Steen | Brain volume in first-episode schizophrenia[END_REF] or reduced hippocampus and schizophrenia [START_REF] Amaral | Neuroanatomy of autism[END_REF]. Estimating properly the impact of the genetic background on neuroanatomical changes is a crucial challenge in order to determine afterwards if this background can either be a risk factor or a protective factor from developing psychiatric disorders. The GWAS studies performed for instance by Stein et al. (2012) identified genetic variants involved in the neuroanatomical diversity, which contributes to understand the impact of genetic factors. However, in the course of these studies, it is shown that this approach only explains a small proportion of the phenotypic variance. In order to understand the nature of the genetic factors responsible for major variations of the brain volume, [START_REF] Toro | Genomic architecture of human neuroanatomical diversity[END_REF] used linear mixed models (LMM) to consider the effects of all the Chapter 3 -Improving heritability estimation by a variable selection approach in high dimensional sparse linear mixed models common genetic diversity characterized by the Single Nucleotide Polymorphisms (SNPs). This approach had been suggested by [START_REF] Yang | GCTA: A tool for genome-wide complex trait analysis[END_REF] to study the effects of the SNPs on the height variations. The model they considered is a LMM defined as follows:

Y = Xβ + Zu + e , (3.1)
where Y = (Y 1 , . . . , Y n ) is the vector of observations (phenotypes), X is a n × p matrix of predictors, β is a p × 1 vector containing the unknown linear effects of the predictors, Z is the genetic information matrix, u and e correspond to the random effects. More precisely, Z is a version of W with centered and normalized columns, where W is defined as follows: W i,j = 0 (resp. 1, resp. 2) if the genotype of the ith individual at locus j is qq (resp. Qq, resp. QQ) where p j denotes the frequency of the allele q at locus j. In (3.1), the vector e corresponds to the environment effects and the vector u corresponds to the genetic random effect, that is the j-th component of u is the effect of the j-th SNP on the phenotype. In the modeling of [START_REF] Yang | GCTA: A tool for genome-wide complex trait analysis[END_REF], all the SNPs have an effect on the considered phenotype, that is

u ∼ N 0, σ u 2 Id R n and e ∼ N 0, σ e 2 Id R n . (3.2)
The covariance matrix of Y can thus be written as:

Var(Y) = N σ u 2 R + σ e 2 Id R n , where R = ZZ N ,
and the parameter η defined as

η = N σ u 2 N σ u 2 + σ e 2 (3.3)
is commonly called the heritability [START_REF] Yang | GCTA: A tool for genome-wide complex trait analysis[END_REF], [START_REF] Pirinen | Efficient computation with a linear mixed model on large-scale data sets with applications to genetic studies[END_REF]), and corresponds to the proportion of phenotypic variance which is determined by all the SNPs. Since all SNPs are not necessarily causal, it seems more realistic to extend the previous modeling by assuming that the genetic random effects can be sparse, that is only a proportion q of the components of u are non null:

u i i.i.d. ∼ (1 -q)δ 0 + qN (0, σ u 2 ), for all 1 ≤ i ≤ N, (3.4)
where q is in (0, 1], and δ 0 is the point mass at 0. Then the definition of η has to be adjusted as follows:

η = N qσ u 2 N qσ u 2 + σ e 2 .
(3.5)

It corresponds to the proportion of phenotypic variance which is due to a certain number of causal SNPs which are, obviously, unknown. Let us emphasize that, in most applications, the proportion q of causal SNPs is also unknown, and that it may happen that the scientist has no idea how small q is. When q = 1, that is when considering the modeling (3.2), most proposed approaches to estimate the heritability derive from a likelihood methodology. We can quote for instance the REstricted Maximum Likelihood (REML) strategies, originally proposed by [START_REF] Patterson | Recovery of inter-block information when block sizes are unequal[END_REF] and then developed in [START_REF] Searle | Variance Components. Wiley Series in Probability and Statistics[END_REF]. Several approximations of the REML algorithm have also been proposed, see for instance the software EMMA proposed by [START_REF] Pirinen | Efficient computation with a linear mixed model on large-scale data sets with applications to genetic studies[END_REF] or the software GCTA [START_REF] Yang | GCTA: A tool for genome-wide complex trait analysis[END_REF], [START_REF] Yang | Common snps explain a large proportion of the heritability for human height[END_REF]).

Chapter 3 -Improving heritability estimation by a variable selection approach in high dimensional sparse linear mixed models

We proposed in [START_REF] Bonnet | Heritability estimation in high-dimensional sparse linear mixed models[END_REF] another method based on a maximum likelihood strategy to estimate the heritability and implemented in the R package HiLMM. We proved in [START_REF] Bonnet | Heritability estimation in high-dimensional sparse linear mixed models[END_REF] the following theoretical result: though the computation of the likelihood is based on the modeling assumption (3.2), the estimator is consistent (unbiased) under the less restrictive modeling assumption (3.4). We believe this consistency result remains true for the estimators produced using the algorithms REML, EMMA, GCTA. But we also proved that, when q = 1, the standard error is not the one computed by the softwares when q = 1 and may be very large. We obtained a theoretical formula for the asymptotic variance of the estimator (depending in particular on q) and conducted several numerical experiments to understand how this asymptotic variance gets larger depending on the various quantities, in particular with respect to q and the ratio n/N . We observed that this variance indeed gets larger when q gets smaller, so that the accuracy of the heritability estimator is slightly deteriorated when all SNPs are not causal. Thus, a first problem is to find a method able to produce an estimator with smaller standard error than those obtained using only likelihood strategies. Also, since this standard error depends on q, a second problem is to produce a confidence interval one could trust without knowing q.

The goal of this paper is to address both problems. The results we obtained in [START_REF] Bonnet | Heritability estimation in high-dimensional sparse linear mixed models[END_REF] suggest the following. If we knew the set of causal SNPs, then, considering only this (small) subset in the genetic information matrix, we would obtain with HiLMM an estimator having a smaller standard error than when using all SNPs in the genetic information matrix. Thus, our new practical method contains a variable selection step.

Variable selection and signal detection in high dimensional linear models have been extensively studied in the past decade and there are many papers on this subject. Among them, we can quote [START_REF] Meinshausen | Stability selection[END_REF] and [START_REF] Beinrucker | Extensions of stability selection using subsamples of observations and covariates[END_REF] about variable selection and references therein. The case of high dimensional mixed models has received little attention. As far as variable selection methods in the random effects of LMM are concerned, we are only aware of the work of [START_REF] Fan | Variable selection in mixed effects models[END_REF] and [START_REF] Bondell | Joint variable selection for fixed and random effects in linear mixed-effects models[END_REF]. Let us mention that regarding the estimation of heritability with possible sparse effects, there is also the bayesian approach of [START_REF] Guan | Bayesian variable selection regression for genome-wide association studies and other large-scale problems[END_REF] and [START_REF] Zhou | Polygenic modeling with bayesian sparse linear mixed models[END_REF], which proposes an interesting estimator for the heritability but which is computationally very demanding. Notice that, in our framework, we are not far from the situation for which it is proved in [START_REF] Verzelen | Minimax risks for sparse regressions: Ultra-high dimensional phenomenons[END_REF] that the support cannot be fully recovered, which happens when N q log(1/q) >> n. The variable selection step we propose takes elements from both ultrahigh dimension methods [START_REF] Fan | Sure independence screening for ultrahigh dimensional feature space[END_REF], [START_REF] Ji | UPS delivers optimal phase diagram in high-dimensional variable selection[END_REF], [START_REF] Meinshausen | Stability selection[END_REF]) and classical variable selection techniques [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]).

The second step of our method is to apply HiLMM using the selected subset of causal SNPs produced by the first step. Finally, we propose a non parametric bootstrap procedure to get confidence intervals with prescribed coverage. The whole procedure requires only a few minutes of computation.

To conclude, we propose in this paper a very fast method to estimate the heritability and construct a confidence interval substantially smaller than without variable selection when the genetic effects are very sparse. Since the real genetic architecture is in general unknown in practice, we also propose an empirical criterion which allows the user to decide whether it is relevant to apply a variable selection based approach or not. Our method has also the advantage to return a list of SNPs possibly involved in the variations of a given quantitative feature. This set of SNPs can further be analyzed from a biological point of view.

The paper is organized as follows. Section 3.2 describes the data set which motivated our

Chapter 3 -Improving heritability estimation by a variable selection approach in high dimensional sparse linear mixed models work. Section 3.3 provides the detailed description of the method, and Section 3.4 displays the results of the numerical study. They were obtained by using the R package EstHer that we developed and which is available from the Comprehensive R Archive Network (CRAN). The simulation results illustrate the performance of our method on simulations and show that it is very efficient from a statistical point of view. In Section 3.5, we provide an empirical criterion to help the user to decide whether it is relevant to apply a variable selection based approach or not. In Section 3.6, we propose a thorough comparison of our approach with other methods in terms of statistical and numerical performances. Finally, the results obtained on the brain data described in Section 3.2 can be found in Section 3.7. We also provide a discussion section at the end of the paper.

Description of the data

We worked on data sets provided by the European project Imagen, which is a major study on mental health and risk taking behaviour in teenagers. The research program includes questionnaires, interviews, behaviour tests, neuroimaging of the brain and genetic analyses. We will focus here on the genetic information collected on approximately 2000 teenagers as well as measurements of the volume of several features: the intracranial brain volume (icv), the thalamus (th), the caudate nucleus (ca), the amygdala (amy), the globus pallidus (pa), the putamen (pu), the hippocampus (hip), the nucleus accubens (acc) and the total brain volume (bv). In the following, our goal will thus be to provide a method to estimate the heritability of these neuroanatomical features.

Chapter 3 -Improving heritability estimation by a variable selection approach in high dimensional sparse linear mixed models

Description of the method

The method that we propose can be split into two main parts: the first one consists in a variable selection approach and the second one provides an estimation of the heritability and the associated 95% confidence interval which is computed by using non parametric bootstrap.

At the beginning of this section we shall consider the case where there is no fixed effects, that is

Y = Zu + e (3.6)
but we explain at the end of this section how to deal with fixed effects. Let us first describe our variable selection method which consists of three steps.

Variable selection

Inspired by the ideas of [START_REF] Fan | Sure independence screening for ultrahigh dimensional feature space[END_REF], we do not directly apply a Lasso type approach since we are in an ultra-high dimension framework. Hence, we start our variable selection stage by the SIS (Sure Independence Screening) approach, as suggested by [START_REF] Fan | Sure independence screening for ultrahigh dimensional feature space[END_REF], in order to select the components of u which are the most correlated to the response Y and then we apply a Lasso criterion which depends on a regularization parameter λ. This regularization parameter is usually chosen by cross validation but here we decided to use the stability selection approach devised by [START_REF] Meinshausen | Stability selection[END_REF] which provided better results in our framework.

Step 1: Empirical correlation computation

The first step consists in reducing the number of relevant columns of Z by trying to remove those associated to null components in the vector u. For this, we use the SIS (Sure Independence Screening) approach proposed by [START_REF] Fan | Sure independence screening for ultrahigh dimensional feature space[END_REF] and improved by [START_REF] Ji | UPS delivers optimal phase diagram in high-dimensional variable selection[END_REF] in the ultra-high dimensional framework. More precisely, we compute for each column j of Z:

C j = Y i Z i,j ,
and we only keep the N max columns of Z having the largest C j . In practice, we choose the conservative value N max = n, inspired by the comments of [START_REF] Fan | Sure independence screening for ultrahigh dimensional feature space[END_REF] on the choice of N max .

In the sequel, we denote by Z red the matrix containing these n relevant columns. This first step is essential for our method. Indeed, on the one hand, it substantially decreases the computational burden of our approach and on the other hand, it reduces the size of the data and thus makes classical variable selection tools efficient.

Step 2: LASSO criterion and stability selection

In order to refine the set of columns (or components of u) selected in the first step and to remove the remaining null components in the vector u, we apply a Lasso criterion originally devised by [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] which has been used in many different contexts and has been thouroughly theoretically studied. It consists in minimizing with respect to u the following criterion:

Crit λ (u) = Y -Z red u 2 2 + λ u 1 , (3.7)
Chapter 3 -Improving heritability estimation by a variable selection approach in high dimensional sparse linear mixed models which depends on the parameter λ and where x 2 2 = p i=1 x 2 i and x 1 = p i=1 |x i | for x = (x 1 , . . . , x p ). The choice of the regularization parameter λ is crucial since its value may strongly affect the selected variables set. Different approaches have been proposed for choosing this parameter such as cross-validation which is implemented for instance in the glmnet R package.

Here we shall use the following strategy based on the stability selection proposed by [START_REF] Meinshausen | Stability selection[END_REF].

The vector of observations Y is randomly split into several subsamples of size n/2. For each subsample, we apply the LASSO criterion for a fixed parameter λ and the selected variables are stored. Then, for a given threshold, we keep in the final set of selected variables only the variables appearing a number of times larger than this threshold. In practice, we generated 50 subsamples of Y and we chose the parameter λ as the smallest value of the regularization path. As explained in [START_REF] Meinshausen | Stability selection[END_REF], such a choice of λ ensures that some overfitting occurs and hence that the set of selected variables is large enough to include the true variables with high probability.

The matrix Z containing only the final set of selected columns will be denoted by Z final in the following, where N final denotes its number of columns.

The threshold has to be chosen carefully: keeping too many columns in Z final could indeed lead to overestimating the heritability and, on the contrary, removing too many columns of Z could lead to underestimating the heritability. In the "small q" situations where it is relevant to use a variable selection approach a range of thresholds in which the heritability estimation is stable will appear as suggested by [START_REF] Meinshausen | Stability selection[END_REF]. In practice, we simulate observations Y satisfying (3.6), by using the matrix Z, for different values of q and for different values η and we observe that this stability region for the threshold appear for small values of q. This procedure is further illustrated in Section 3.4.

Heritability estimation and confidence interval

Heritability estimation

For estimating the heritability, we used the approach that we proposed in [START_REF] Bonnet | Heritability estimation in high-dimensional sparse linear mixed models[END_REF]. It is based on a maximum likelihood strategy and was implemented in the R package HiLMM. Let us recall how this method works.

In the case where q = 1, which corresponds to the non sparse case,

Y ∼ N 0, η σ 2 R + (1 -η )σ 2 Id R n , with σ 2 = N σ 2 u + σ 2 e and R = Z final Z final /N final ,
where Z final denotes the matrix Z in which the columns selected in the variable selection step described in Section 3.3.1 are kept.

Let U be defined as follows: U U = UU = Id R n and URU = diag(λ 1 , . . . , λ n ), where the last quantity denotes the diagonal matrix having its diagonal entries equal to λ 1 , . . . , λ n . Hence, in the case where q = 1,

Y = U Y ∼ N (0, Γ) with Γ = diag(η σ 2 λ 1 + (1 -η )σ 2 , . . . , η σ 2 λ n + (1 -η )σ 2 ), (3.8)
where the λ i 's are the eigenvalues of R.

Chapter 3 -Improving heritability estimation by a variable selection approach in high dimensional sparse linear mixed models

We propose to define η as a maximizer of the log-likelihood

L n (η) = -log 1 n n i=1 Y 2 i η(λ i -1) + 1 - 1 n n i=1 log (η(λ i -1) + 1) , (3.9)
where the Y i 's are the components of the vector Y = U Y.

We now explain how to obtain accurate confidence intervals for the heritability by using a non parametric bootstrap approach.

Bootstrap confidence interval

We used the following procedure:

-Step 1: We estimate η and σ 2 by using our approach described in the previous subsection.

The corresponding estimators are denoted η and σ.

-Step 2: We compute Y new = Γ-1/2 Y, where Y is defined in (3.8) and Γ has the same structure as Γ defined in (3.8) except that η and σ are replaced by their estimators η and σ, respectively.

-Step 3: We create K vectors (Y new,i ) 1≤i≤K from Y new by randomly choosing each of its components among those of Y new .

-Step 4: We then build K new vectors ( Y samp,i ) 1≤i≤K as follows: Y samp,i = ΓY new,i . For each of them we estimate the heritability. We thus obtain a vector of heritability estimators (η 1 , ..., ηK ).

-Step 5: For obtaining a 95% bootstrap confidence interval, we order these values of ηk and keep the ones corresponding to the 0.975 × K largest and the 0.025 × K smallest, where x denotes the integer part of x. These values define the upper and lower bounds of the 95% bootstrap confidence interval for the heritability η , respectively.

A bootstrap estimator of the variance can be obtained by computing the empirical variance estimator of the ηk 's. In practice, we chose K = 80 replications.

In Step 2 of the previous algorithm, we should be in the non sparse case q = 1 thanks to the variable selection stage. Hence, the covariance matrix of Y new should be close to identity.

Observe that our resampling technique is close to the one proposed by Abney (2015) for building permutation tests in linear mixed models.

Additional fixed effects

The method described above does not take into account the presence of fixed effects. For dealing with such effects we propose to use the following method, which mainly consists in projecting the observations onto the orthogonal of Im(X), the image of X, to get rid of the fixed effects. In practice, instead of considering Y and Z we consider Ỹ = A Y and Z = A Z, where A is a n × (n -d) matrix (d being the rank of the fixed effects matrix), such that AA = P X , A A = Id R n-d and P X = Id R n -X(X X) -1 X . This procedure was for instance used by [START_REF] Fan | Variable selection in mixed effects models[END_REF].

Chapter 3 -Improving heritability estimation by a variable selection approach in high dimensional sparse linear mixed models

Numerical study

We present in this section the numerical results obtained with our approach which is implemented in the R package EstHer.

Simulation process

Since in genetic applications, the number n of individuals is very small with respect to the number N of SNPs, we chose n = 2000 and N = 100000 in our numerical study. We also set σ 2 u = 1, we shall consider different values for q and we shall change the value of σ e in order to have the following values for η : 0.4, 0.5, 0.6 and 0.7. We generate a matrix W such that its columns W j are independent binomial random variables of parameters n and p j , where p j is randomly chosen in [0.1, 0.5]. We compute Z by centering and empirically normalizing the matrix W. The random effects are generated according to Equation (3.4) and then we compute a vector of observations such that Y = Zu + e.

We can make two important comments about the previous simulation process. Firstly, we generated a matrix W with independent columns, that is we assume that the SNPs are not correlated. Since this assumption may not be very realistic in practice, we provide in Section 3.4.2 some additional simulations where the generated matrix W has been replaced by the real matrix W coming from the IMAGEN project. Secondly, we did not include fixed effects but we show some results in Section 3.4.2 when fixed effects are taken into account.

Results in very sparse scenarios

In this section, we shall focus on the performances of our method in a very sparse scenario, that is 100 causal SNPs out of 100,000. We will describe all the results in terms of heritability estimation, support recovery and computational times in this particular case, then we will study other sparsity scenarios.

Choice of the threshold

In order to determine the threshold, we apply the procedure described in Section 3.3.1 and 3.3.2. Figure 3.2 displays the mean of the absolute value of the difference between η and the estimated value η for different thresholds and for different values of η obtained from 10 replications. We can see from this figure that in the case where the number of causal SNPs is relatively small: 100, that is q = 10 -3 , our estimation procedure provides relevant estimations of the heritability for a range of thresholds around 0.75. Moreover, the optimal threshold leading to the smallest gap between η for different values of η is 0.76. We will use this value in the following numerical study. However, the way of choosing the threshold will be further discussed, especially in the section dedicated to the study of the genetic data.

Confidence intervals

We use the non parametric boostrap approach described in Section 3.3 in order to compute the confidence intervals associated to the estimations of the heritability. Table 3.1 shows that the 95% confidence intervals obtained by bootstrap and the empirical confidence intervals are very similar. The empirical confidence intervals are computed as follows: the different estimations .2 -Absolute difference between η and η for thresholds from 0.6 to 0.9 and for q = 10 -3 (100 causal SNPs). of η obtained along the different replications are ordered, the 0.975 × M largest and the 0.025 × M smallest values correspond to the upper (resp. lower) bound of the 95% empirical confidence interval. Here, x denotes the integer part of x and M is the number of replications. From Table 3.1, we can see that the empirical confidence intervals are included in the bootstrap intervals, which means that our approach provides conservative intervals.

Comparison between the methods with and without selection

Our results are compared to those obtained if we do not perform the selection before the estimation, that is with the method implemented in HiLMM ("without"), but also with an approach which assumes the position of the non null components to be known (oracle). The results are displayed in Figure 3.3 and in Table 3.2. In this table, the confidence intervals displayed for the lines "Oracle" and "without" are obtained by using the asymptotic variance derived in [START_REF] Bonnet | Heritability estimation in high-dimensional sparse linear mixed models[END_REF] which corresponds to the classical inverse of the Fisher information in the case q = 1. We observe that our method without the selection step provides similar results, that is almost no bias but a very large variance due to the framework N n. Our method EstHer considerably reduces the variance compared to this method and exhibits performances close to those of the oracle approach which, contrary to our approach, knows the position of the non null components.

Additional fixed effects

We generated some synthetic data according to the process described in Section 3.4.1 but we added a matrix of fixed effects containing two colums. Figure 3.4 (a) displays the corresponding Chapter 3 -Improving heritability estimation by a variable selection approach in high dimensional sparse linear mixed models Table 3.2 -95 % confidence intervals for η obtained by our approach, GCTA, the oracle approach and the approach without selection ("without"). (c) (d)

η
Figure 3.3 -Estimation of the heritability and the corresponding 95% confidence intervals when q =10 -3 , and for different values of η : (a) η = 0.4, (b) η = 0.5, (c) η = 0.6, (d) η = 0.7. The means of the heritability estimators (displayed with black dots), the means of the lower and upper bounds of the 95% confidence intervals are obtained from 20 replicated data sets for the different methods: without selection ("without"), "oracle" which knows the position of the null components and EstHer. The horizontal gray line corresponds to the value of η .

results which show that the presence of fixed effects does not alter the heritability estimation.

Simulations with the matrix W of the IMAGEN data set

We conducted some additional simulations in order to see the impact of the linkage disequilibrium, that is the possible correlations between the columns of Z. Indeed, in the previous numerical study, we generated a matrix W with independent columns. The matrix W that we use now to generate the observations is the one from our genetic data set, except that we truncated it in order to have n = 2000 and N = 100000. The results of this additional study Chapter 3 -Improving heritability estimation by a variable selection approach in high dimensional sparse linear mixed models are presented in Figure 3.4 (b). We can see that they are similar to those obtained previously in Figure 3.3, which means that our method does not seem to be sensitive to the presence of correlation between the columns of W.

Computational times

The implementation that we propose in the R package EstHer is very efficient since it only takes 45 seconds for estimating the heritability and 300 additional seconds to compute the associated 95% confidence interval. These results have been obtained with a computer having the following configuration: RAM 32 GB, CPU 4 × 2.3 GHz.

Recovering the support

When the number of causal SNPs is reasonably small, our variable selection method is efficient to estimate the heritability and we wonder if it is reliable as well to recover the support of the random effects. In Figure 3.5, we see the proportion of support estimated by our method when there are 100 causal SNPs: our method selects around 130 components. We then focus on the proportion of the real support which has been captured by our method: we see that it may change according to η . Indeed, the higher η , the higher this proportion. Nevertheless, even in the worst case, that is η = 0.5, Figure 3.6 shows that even if we keep only 30% of the real non null components, we select the most active ones. The ability of recovering the support in linear models has been studied by Verzelen (2012) in ultra high dimensional cases. The author shows that with a non null probability, the support cannot be estimated under some numerical conditions on the parameters q, N and n (namely if there are considerably more variables N than observations n, and if the number of non null components qN is relatively high). In this simulation study, even when we consider small Chapter 3 -Improving heritability estimation by a variable selection approach in high dimensional sparse linear mixed models values of q (for instance q = 10 -3 , that is 100 causal SNPs), we are not far from to the ultra high dimensional framework described in [START_REF] Verzelen | Minimax risks for sparse regressions: Ultra-high dimensional phenomenons[END_REF], which can explain the difficulties to recover the full support.

Results when the number of causal SNPs is high

In subsection 3.4.2 we show the performance of our method in the case where the proportion of causal SNPs q is small, that is around 10 -3 . In this subsection, we focus on a more polygenic scenario, that includes the cases where thousands of SNPs or ten of thousands of SNPs are causal.

Results when there are SNPs with moderate and weak effects

We first focus on the statistical performance of EstHer when there are a lot of SNPs (1000 or 10000) with small effects (for example, that explain 5% of the phenotypic variations), and a small number (around 100) with moderate effects. We can see from Figure 3.7 that, in this case, EstHer provides unbiased estimations with a small variance.

Results when all SNPs have moderate effects

If all causal SNPs have moderate effects and if the number of these causal SNPs is high, namely greater than 1000, EstHer underestimates the heritability. These results are displayed in Figure 3.8. Moreover, we can see from Figure 3.9 that there is no threshold choice that can provide accurate estimations of heritability for all values of η .

Chapter 3 -Improving heritability estimation by a variable selection approach in high dimensional sparse linear mixed models Figure 3.6 -Barplots of the proportion of components found by our method as function of the most efficient variables. For example, the first bar is the proportion of the 10 % higher components that we captured with our selection method. The histograms are displayed for several values of η : 0.5 (a), 0.6 (b), 0.7 (c).

A criterion to decide whether we should apply EstHer or HiLMM

On the one hand, we observed that applying HiLMM provides unbiased estimations of the heritability, no matter the number of causal SNPs. However, the main drawback of this estimator is its very large variance. On the other hand, if the number of causal SNPs is not too high, EstHer provides unbiased estimations of the heritability with standard errors substantially smaller than HiLMM. However, if the number of causal SNPs is high, EstHer underestimates the heritability. These observations are similar to those made by [START_REF] Zhou | Polygenic modeling with bayesian sparse linear mixed models[END_REF], who built an hybrid estimator able to deal with both sparse and non sparse scenario, to which we will compare our approach in Section 3.6. Therefore, we propose hereafter a rule to decide whether it is better to apply EstHer or HiLMM. We can see from Figure 3.2 that when there are 100 causal SNPs, there is a large range of threshold values which provide an accurate estimation of η , but when there are 1000 or 10000 causal SNPs, see Figure 3.9), the estimations are very different even for close thresholds. This observation gave us the idea of quantifying the stability of the estimations around the threshold that we determined as the optimal one. More precisely, for each threshold, we have an estimation of heritability with a 95% confidence interval, and we count the number of thresholds for which the confidence intervals overlap. Figure 3.10 confirms the stability around the best threshold for different values of η and Table 3.3 displays the number of ovelapping confidence intervals. We empirically determine the following criterion: if the mean number of thresholds is greater than 10 (over 16 tested thresholds), we apply EstHer, if not, we apply HiLMM. The results obtained by using this criterion are displayed in Figure 3.11. 3.9 -Absolute difference |η -η| for thresholds from 0.6 to 0.9 and for 1000 (left) and 10000 (right) causal SNPs.

Results after applying the decision criterion and comparison to other methods

Statistical performances

In this section we show the results obtained after applying the criterion described in Section 3.5. We compare these results to those obtained using HiLMM, but also with the software GEMMA described in [START_REF] Zhou | Genome-wide efficient mixed model analysis for association studies[END_REF] q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 n=2000, N=10^5, 100 causaux without 0.71 0.73 0.75 0.77 0.79 0.81 0.83 0.85 q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 n=2000, N=10^5, 100 causaux without 0.71 0.73 0.75 0.77 0.79 0.81 0.83 0.85 q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 n=2000, N=10^5, 1000 causaux without 0.71 0.73 0.75 0.77 0.79 0.81 0.83 0.85 q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 n=2000, N=10^5, 1000 causaux without 0.71 0.73 0.75 0.77 0.79 0.81 0.83 0.85 q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 n=2000, N=10^5, 1000 causaux without 0.71 0.73 0.75 0.77 0.79 0.81 0.83 0.85 q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 n=2000, N=10^5, 10000 causaux without 0.71 0.73 0.75 0.77 0.79 0.81 0.83 0.85 q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 n=2000, N=10^5, 10000 causaux without 0.71 0.73 0.75 0.77 0.79 0.81 0.83 0.85 q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 n=2000, N=10^5, 10000 causaux without 0.71 0.73 0.75 0.77 0.79 0.81 0.83 0.85 q q q q q q q q q q q q q q q q q Figure 3.10 -Estimation of the heritability with 95% confidence intervals for η from 0.4 to 0.6 (from left to right), and from 100, 1000 and 10000 causal SNPs from top to bottom. Each graph shows the heritability estimations with 95% confidence intervals computed with HiLMM ("without") and for thresholds between 0.7 and 0.85.

BSLMM in the sequel. As explained in [START_REF] Zhou | Polygenic modeling with bayesian sparse linear mixed models[END_REF], BSLMM can deal with very sparse and also with very polygenic scenarios. We can see from the bottom part of Figure 3.11 that, in very polygenic scenarios (q = 0.1, namely 10,000 causal SNPs), all the methods provide similar results: the four estimators are indeed empirically unbiased, but with a very large variance.

In sparse scenarios (q = 10 -3 , namely 100 causal SNPs), we can see from the top part of Figure 3.11 that EstHer provides better results than HiLMM and GEMMA-LMM which exhibit similar statistical performances. In sparse scenarios, the variance of the BSLMM estimator is larger than the one provided by EstHer and smaller than the one provided by GEMMA-LMM and HiLMM. However, the performances of BSLMM could perhaps be improved by changing Chapter 3 -Improving heritability estimation by a variable selection approach in high dimensional sparse linear mixed models the MCMC parameters. Here, for computational time reasons, we used the default parameters that is 100,000 and 1,000,000 for the number of burn-in steps and the number of sampling, respectively. The results are obtained with 10 replications.

Computational times

The computational times in seconds for one estimation of the heritability with BSLMM and the heritability estimation for 16 thresholds as well as the associated confidence intervals with our method EstHer are displayed in Figure 3.12. We chose this number of thresholds since we applied the criterion defined in Section 3.5. It should be noticed that the computational times for EstHer could be reduced by diminishing the number of thresholds. For BSLMM we used the default parameters for the number of burn-in steps and the number of sampling. We can see from this figure that the gap between EstHer and BSLMM is all the more important that N is large. Contrary to our approach, BSLMM seems to be very sensitive in terms of computational time to the value of N .

Applications to genetic data

In this section, we applied our method to the neuroanatomic data coming from the Imagen project. In this data set, n = 2087 individuals and N = 273926 SNPs. For further details on this data set, we refer the reader to Section 3.2.

Chapter 3 -Improving heritability estimation by a variable selection approach in high dimensional sparse linear mixed models

Calibration of the threshold

We start by finding the threshold which is the most adapted to the Imagen data set. We use the same technique as the one described in Section 3.4.2: for several values of η and several thresholds, we display the absolute value of η -η, see Figure 3.13. The only difference with Section 3.4.2 is that we generated the observations by using the matrix W coming from the IMAGEN data set. According to Figure 3.13, we can find a reliable range of thresholds for estimating the heritability for all η from 0.4 to 0.7 when the number of causal SNPs is smaller than 100. This optimal threshold is equal to 0.79. We shall use this value in the sequel. Chapter 3 -Improving heritability estimation by a variable selection approach in high dimensional sparse linear mixed models

Application of the decision criterion

Since we determined in the previous section that the optimal threshold is 0.79, we apply Es-tHer for thresholds around this value, that is from 0.7 to 0.85. We then count the number of overlapping confidence intervals, as explained in Section 3.5. The results are displayed in Table 3.4. We observe from this table that the sensitivity to the choice of the threshold varies substantially from one phenotype to another. Hence, we choose to apply our EstHer approach to the most stable phenotypes with respect to our criterion, namely pa, amy and acc. For the other phenotypes we recommand to apply HiLMM or another similar approach such as GCTA or GEMMA-LMM.

Results

Figure 3.14 (a) shows the heritability estimation with 95 % confidence intervals for all phenotypes, using either EstHer or HiLMM according to the outcome of our decision criterion. Figure 3.14 (b) shows the results obtained by using HiLMM, namely without any variable selection step. We compare our results with the ones obtained by [START_REF] Toro | Genomic architecture of human neuroanatomical diversity[END_REF] who estimated the heritability of the same phenotypes by using the software GCTA. On the one hand, we can see from Figure 3.14 that in the cases where EstHer is used the confidence intervals given by our methodology are substantially smaller and included in those provided by either HiLMM or [START_REF] Toro | Genomic architecture of human neuroanatomical diversity[END_REF]. On the other hand, when HiLMM is used our results are on a par with those obtained by [START_REF] Toro | Genomic architecture of human neuroanatomical diversity[END_REF]. Moreover, our approach provides a list of SNPs which may contribute to the variations of a given phenotype and which could be further analyzed from a biological point of view in order to identify new biological pathways.

Conclusion

We propose in this paper a practical method to estimate the heritability in sparse linear mixed models using variable selection tools, as well as confidence interval obtained thanks to a non parametric bootstrap approach. Our approach is implemented in the R package EstHer which is available from the Comprehensive R Archive Network (CRAN) and from the web page of the first author. In the course of this study, we showed that our approach has two main features q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 bv icv th pu pa hip amy acc ca q q q q q q q q q (a) (b)

Figure 3.14 -(a) Heritability estimations of bv, icv, th, pu, pa, hip, amy, acc, and ca with 95% confidence intervals obtained using EstHer or HiLMM according to the outcome of our decision criterion. (b) Heritability estimations of bv, icv, th, pu, pa, hip, amy, acc and ca with 95% confidence intervals obtained using HiLMM.

which makes it very attractive. Firstly, it is very efficient from a statistical point of view since it provides confidence intervals considerably smaller than those obtained with methods without variable selection. Secondly, its very low computational burden makes its use feasible on very large data sets coming from quantitative genetics. Moreover, we observed that the statistical performance of the EstHer approach are all the more impressive that the level of sparsity is high that is when q is small. For this reason, we also proposed an empirical criterion which allows the user to decide whether it is better to apply an approach that takes into account the sparsity and starts with a variable selection stage, namely EstHer, or an approach which ignores the potential sparsity in the observations, namely HiLMM.

Chapter 4 Application: heritability estimation of the size of juvenile trouts

The work presented in this chapter comes from a collaboration with Niklas Tysklind, who is a researcher in the UMR EcoFoG, INRA in Kourou. This work will be soon submitted to an international journal of ecology.

Data and motivation

We have a dataset containing the size of 142 juvenile trouts, the genotype of which is described by 3917 SNPs. We also have several possible explaining variables: the river they lived in, the month the trouts have been captured, and several features regarding the environmental conditions: latitude, longitude, alkalinity, hardness, catchment area (it gives an idea of the potential amount of environment available for the sea trout), shreve, strahler (both are complexity values of the river in terms of how many affluents and branches each river has). Finally we have an additional information called the "mean +0 and > 0+" juvenile densities, which are included as a proxy for the effect of competition between the juvenile trouts.

Our goal to estimate the heritability of the trouts size, that is the proportion of variability of the size which can be explained by genetic factors.

Model and method

Model

We propose to use the following sparse linear mixed model:

Y = Xβ + Zu + e
where • Y is a vector of length 142 which contains the observations of the size of the trouts.

• X is a matrix associated to the fixed effects: it contains the explaining variables (river, month...).

• β corresponds to the coefficients associated to X.

• Z is a n × N matrix where Z i,j is the value of the SNP j of the individual i.

• u and e are the random effects. We assume that e ∼ N 0, σ e 2 Id R n , and since we do not know the SNPs which may have an effect on the observations, we shall assume that

u i i.i.d.
∼ (1 -q)δ 0 + qN (0, σ u 2 ) , for all i = 1..n where u = (u 1 , . . . , u n ).

Method

The method that we used is described in Chapter 3. However, the first step of the procedure, called the Sure Independance Screening (SIS), is specific to ultra high dimensional datasets.

According to numerical results obtained with and without the SIS step, which are not presented here, we propose to skip the first step of our method. The results of this approach are shown in the following simulation study.

Simulation Study

We propose a simulation study in the following framework: n = 200, N = 4000, which are close to the size of the trouts dataset that we want to analyze. We assume that the number of causal SNPs varies from 4 (q = 10 -3 ) to 400 (q = 0.1).

Results obtained without selection

The estimations obtained using our method HiLMM, described in Chapter 2, are shown in Figure 4.1. The large empirical variance of this estimator motivates an additional variable selection step before estimating heritability.

Calibration of the threshold in the stability selection

As explained in the detailed procedure in Chapter 3, our method requires to determine the optimal threshold in the stability selection. For this purpose, we apply the procedure described in Section 4.2.2. Figure 4.2 displays the mean of the absolute value of the difference between η and the estimated value η for different thresholds and for different values of η obtained from 10 replications. When the number of causal SNPs is small (namely, smaller than 20), applying our selection method with a threshold around 0.35 seems to provide accurate estimations for different values of the heritability η . Otherwise, there is no common threshold that provides accurate estimations for all values of the heritability η . In this scenario, we shall estimate directly heritability without performing a variable selection approach.

Results regarding heritability estimation

We observe firstly from Figure 4.3 that when the number of causal SNPs is small (namely, smaller than 20), applying our selection method substantially improves the heritability estimations compared to the method without selection. Secondly, when the number of causal SNPs is high, both methods have similar performances.

Support recovery when the selection is applied

When the number of causal SNPs is small, our criterion leads us to apply the variable selection approach. We focus on the support that we captured with our method: Figure 4.4 shows the Chapter 4 -Application: heritability estimation of the size of juvenile trouts 4 causal SNPs 8 causal SNPs q q q q q q q EstHer HiLMM 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q EstHer HiLMM 0.0 0.2 0.4 0.6 0.8 1.0 20 causal SNPs 40 causal SNPs q q q q q q q q q q q q q EstHer HiLMM proportion of the recovered support and the number of variables selected by our method. We notice that choosing a threshold around 0.35 seems to be a relevant trade-off between a large true positive rate and a small false positive rate of non zero recovered components. We recover indeed at least 30% of the support and up to 60% when the number of causal SNPs is very small Chapter 4 -Application: heritability estimation of the size of juvenile trouts (namely, 4 or 8 SNPs).

4 causal SNPs 8 causal SNPs 20 causal SNPs q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 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.2 0.3 0.4 0.5 0.6 0.7 0.8 0.0 0.2 0.4 0.6 0.8 1.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 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.0 0.2 0.4 0.6 0.8 q q q q q q q q q q q q q q 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.0 0.1 0.2 0.3 0.4 0.5 q q q q q 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 10 20 30 40 50 60 q q q q q q q q q q 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.2 0.3 0.4 0.5 0.6 0.7 0.8 0 10 20 30 40 50 60 q q q q q q q q 0.2 0.

Results on the juvenile size dataset 4.4.1 Results obtained without taking fixed effects into account

First we propose to estimate the heritability without taking any fixed effect into account. Figure 4.5 shows firstly that the heritability of the trouts size is very high (80%) despite the large confidence interval (40-100%). Moreover, we can see that EstHer and HiLMM provide the same estimation, which suggests that the number of causal SNPs must be high.

Including fixed effects

We estimated the heritability when considering as a fixed effect the river where the trouts were captured in. The results are shown in Figure 4.6: including the river as a fixed effect provides a higher heritability estimation (0.93 instead of 0.79) and more accurate (the length of the confidence interval is substantially smaller). It means that the river and the SNPs explain a very large proportion of the size variations. Including the month as a fixed effect does not seem to have a relevant effect: a possible explanation is that the month effect was already taken into account in the river effect, since all trouts of the same river have been captured the same month. Then, we considered phylogeographic groups instead of considering the specific river: the 30 rivers are contained in 9 phylogeographic groups with similar environmental features. The Figure 4.5 -Estimated value of the heritability with 95 % confidence intervals, using EstHer and HiLMM.

results using these phylogeographic groups as a fixed effect are displayed in Figure 4.6. It seems that these groups bring less information than the rivers, but if we join the groups with the month effect, we have similar results than those obtained with either the river or the combination river and month. This could mean that the information we lost when considering the groups instead of the rivers was contained in the month effect. Then, we wonder if we can determine the specific features of the river that can explain the increasing proportion of variability explained when knowing the river. Figure 4.7 shows the results of considering different groups of fixed effects. The groups are the following:

-Group 1: longitude, latitude, catchment area, shreve, strahler

Chapter 4 -Application: heritability estimation of the size of juvenile trouts -Group 2: longitude, latitude, catchment area, shreve, strahler, alkalinity, hardness -Group 3 : longitude, latitude, catchment area, shreve, strahler, mean 0+ and mean > 0+ These groups were chosen according to the availability of the data: indeed, in the first group we have all the information for the 142 trouts. Alkalinity and hardness were missing for 12 individuals, and mean 0+ and mean > 0+ for 16 individuals (unfortunately, not the same than alkalinity and hardness).

The results are displayed in Figure 4.7. We see that Group 1 explains a proportion of variability bigger than no fixed effects but smaller than the river. Same result when we add mean 0+ and mean > 0+. The confidence interval increases when considering the additional effect of alkalinity and hardness. However, this could possibly be the effect of removing individuals from the study. This assumption is consistent with Figure 4.8, which shows the results obtained with the variables of Group 1 but without the individuals that we had to remove when we studied Group 2, compared to the results obtained with Group 2. Similarly, Figure 4.9 shows the results obtained with the variables of Group 1 but without the individuals that we had to remove when we studied Group 3, compared to the results obtained with Group 3. We see a substantial difference in the estimations, which suggests that mean 0+ and mean > 0+ could be relevant variables. no cov gp1/month gp2/month gp3/month q q q q q Figure 4.7 -Estimated value of the heritability with different covariates from left to right: river/month, no covariates, group 1/month, group 2/month, group 3/month

Conclusion

• The length variations have an important genetic component (≥40%).

• This genetic component might be due to large number of causal SNPs (at least, more than 20 according to the simulation study), which makes impossible the recovery of these causal SNPs.

• Including fixed effects can refine the heritability estimations: the river seems indeed to be a key factor in the length variation. Once the river has been taken into account, more Figure 4.9 -Estimated value of the heritability with and without mean 0+ and mean > 0+ as fixed effects, but removing 16 individuals from the study.

than 75% of the remaining variation is due to genetic effects. However, it is important to notice that we cannot separate entirely the effects of the month and the river.

• The reporting groups alone do not improve the heritability estimations, but associated to the month effect, they explained as much variability as either the river, or the combined effect river/month.

• In the environmental factors, the combined effect of catching area, longitude, latitude, shreve, strahler explain more than no fixed effects, but less than the river.

Chapter 5

Heritability estimation in case-control studies

The work contained in this chapter will be soon submitted to a journal of statistics.

Abstract

In the genetic field, the concept of heritability refers to the proportion of variations of a biological trait or disease that can be explained by genetic factors. Quantifying the heritability of a disease is a fundamental challenge in human genetics. Although the litterature regarding heritability estimation for binary traits is less rich than for quantitative traits, several methods have been proposed to estimate the heritability of diseases. However, to the best of our knowledge, the existing methods raise at least one of the two concerns (generally both): either they have not been validated theoretically or they are severely affected by the oversampling of the number of patients compared to controls in a medical study. We propose in this paper to investigate the theoretical properties of the method developed by [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF], which is very efficient in practice, despite the oversampling of patients. Our main result is the proof of the consistency of this estimator. We also provide a numerical study to compare two approximations leading to two heritability estimators.

Introduction

In the genetic field, the concept of heritability refers to the proportion of variations of a biological trait or disease that can be explained by genetic factors. Quantifying the heritability is a major information for diseases that are suspected to have a strong genetic component but which causes are often vague and multiple. Indeed, determining a high value of heritability for a disease is a powerful argument in favor of further research for genetic causes, but it also opens the possibility of predicting a risk of illness based on the genetic background.

There exist several methods to estimate the heritability of quantitative traits, which we will describe hereafter, with interesting theoretical and practical properties. Regarding binary traits, such as the presence or absence of a disease, a few methodologies have been proposed, but as far as we know, most of them have not been validated theoretically. [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF] developed a method to estimate heritability of binary traits that they compared to recent methodologies and which was shown to be very efficient in practice. The aim of this paper is to investigate the theoretical properties of [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF]'s method.

Let us first recall the main existing methods to estimate heritability of quantitative traits, which will be strongly linked to the methods used for binary traits. Linear Mixed Models (LMMs) have been widely used for estimating the heritability of quantitative traits. Indeed, [START_REF] Yang | Common snps explain a large proportion of the heritability for human height[END_REF] proposed for instance to estimate the heritability of human height by using a classical LMM defined by

Y = Xβ + Zu + e, (5.1)
where Y = (Y 1 , . . . , Y n ) is the vector of observations of a phenotype of interest, X is a n × p matrix of predictors (or fixed effects), β is a p × 1 vector containing the unknown linear effects of the predictors, and u and e correspond respectively to the genetic and the environmental random effects. We assume that u and e are Gaussian random effects with variances σ 2 u and σ 2 e respectively. Moreover, Z is a n × N matrix which contains the genetic information. They proposed to estimate the parameter

η = N σ 2 u N σ 2 u + σ 2 e , (5.2)
commonly considered as the mathematical definition for heritability since it determines how the variance is shared between u and e. Several methods were developed to estimate the parameter η , see [START_REF] Patterson | Recovery of inter-block information when block sizes are unequal[END_REF], [START_REF] Searle | Variance Components. Wiley Series in Probability and Statistics[END_REF], [START_REF] Yang | GCTA: A tool for genome-wide complex trait analysis[END_REF][START_REF] Pirinen | Efficient computation with a linear mixed model on large-scale data sets with applications to genetic studies[END_REF], [START_REF] Zhou | Genome-wide efficient mixed model analysis for association studies[END_REF]. From a theoretical point of view, [START_REF] Bonnet | Heritability estimation in high-dimensional sparse linear mixed models[END_REF] showed the asymptotic normality of the maximum likelihood estimator of η as well as a central limit theorem leading to confidence population that is affected by the disease (when the disease is rarer, the bias increases). The estimator also underestimates the heritability when the true heritability is high. [START_REF] Weissbrod | Accurate liability estimation improves power in ascertained case-control studies[END_REF] introduced a maximum likelihood based strategy to rebuild the underlying liability before estimating the heritability. However, all these methods do not take into account an essential element of case-control studies: indeed, in a medical study, the number of patients is similar to the number of controls even though the studied disease might be rare, which means that the proportion of cases in the study does not reflect the proportion of cases in the population. This oversampling of the cases has been noticed and handled by the approach of [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF], who proposed a moment based method to estimate the heritability. They computed an approximate quantity of the expectation of W i W j , for two individuals i and j, W i being a centered and normalized version of the binary data Y i , and conditionally to the fact that individuals i and j are in the study.

Since the method of [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF] presented very good numerical results but was not supported by theoretical grounds, we propose in this paper a method which is strongly inspired from the one proposed by [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF]: we obtained two approximations (depending on the order of the approximation) of the expectation E of W i W j conditionally to the fact that both individuals are in the study. We show that the first order approximation provides an estimator with good theoretical properties: indeed, we prove that it is a consistent estimator of η . We also propose a simulation study to compare the numerical performances of the estimators obtained with both approximations. The model we study and the main definitions are given in Section 5.2. Section 5.3 contains the first order approximation of the expectation E with the corresponding estimator of η and Section 5.4 presents our consistency result for this estimator. The second order approximation of E is given in Section 5.5 and the numerical comparison of the two estimators can be found in Section 5.6. In Section 5.7, we discuss the results and potential perspectives. Finally, the proofs are given in Section 5.8.

Model and definitions

Liability model

Let us denote K the prevalence of a disease in a population, that is the proportion of the population affected by the disease. Let Y i be the random variable such that Y i = 1 if the individual i is ill (then, individual i is called a case) and Y i = 0 if the individual i is healthy (then individual i is called a control). We assume that the Y i 's are linked to unobserved variables l i as follows

Y i = 1 {l i >t} , (5.7)
where t is a given threshold, related to the prevalence K by (5.6), and the l i 's are defined as

l = Zu + e, (5.8)
where l = (l 1 , . . . , l n ), u and e are random effects such that u ∼ N (0, σ 2 u Id R N ) and e ∼ N (0, σ 2 e Id R n ). The vector u corresponds to the genetic effects and e to the environmental effects. Moreover, Z is a n × N random matrix which contains the genetic information, and which is such that the Z i,k are normalized random variables in the following sense: they are defined from a matrix A = (A i,k ) 1≤i≤n, 1≤k≤N by

Z i,k = A i,k -A k s k , i = 1, . . . , n, k = 1, . . . , N , (5.9)
where

A k = 1 n n i=1 A i,k , s 2 k = 1 n n i=1 (A i,k -A k ) 2 , k = 1, . . . , N .
(5.10)

In (5.9) and (5.10) the A i,k 's are such that for each k in {1, . . . , N } the (A i,k ) 1≤i≤n are independent and identically distributed random variables and such that the columns of A are independent.

In practice, the matrix A contains the genetic information about all the individuals in the study. More precisely, for each k, A i,k = 0 (resp. 1, resp. 2) if the genotype of the ith individual at locus k is qq (resp. Qq, resp. QQ). In this paper, we consider a more general case with mild assumptions on the distribution of the random variables A i,k , which are described in Section 5.4. With the definition (5.9), the columns of Z are empirically centered and normalized, and one can observe that

Var(l|Z) = N σ u 2 R + σ e 2 Id R n , where R = ZZ N .
The heritability at the liability scale, which is the parameter we want to estimate, is defined as the ratio of variances:

η = N σ 2 u N σ 2 u + σ 2 e .
(5.11)

The variance of l conditionally to Z can then be rewritten with respect to η and σ 2 = N σ 2 u +σ 2 e as:

Var(l|Z) = η σ 2 R + (1 -η )σ 2 Id R n .
We will assume in the sequel without loss of generality that σ 2 = 1. Indeed, if σ 2 = 1, we can consider the variable l i = l i σ and then, instead of estimating t from the prevalence K with the relationship (5.6), we estimate directly t/σ .

Case control study

Since the prevalence P in the study can be very different from the prevalence K in the general population (the cases are substantially oversampled in a case-control study), it is essential to consider that the observations that we have access to depend on the probabilities for both cases and controls to be selected in the study. Indeed, if p control denotes the probability for a control to be selected in the study, we can define the corresponding variable U i ∼ B(0, p control ) which is equal to 1 if individual i is a control who is selected in the study. Similarly we define the probability p case for a case to be selected for the sudy and the corresponding variable V i ∼ B(0, p case ). Then for any individual i, we define the variable i by which is equal to 1 if individual i belongs to the study and 0 if not. We assume that the variables U 1 , . . . , U n , V 1 , . . . , V n are independent and independent of Y 1 , . . . , Y n and Z. Since we do not observe Y i for the whole population but only for the individuals who belong to the study, we will work with the variables W i defined by

W i = Y i -P P (1 -P ) i ,
which are centered versions of Y i in the study and are non-zero only if individual i belongs to the study. The probabilities p case and p control are chosen such that the prevalence in the study is equal to P . Indeed, if we assume that p case = 1, (5.12) it implies that

p control = K(1 -P ) P (1 -K) .
(5.13)

The proof of (5.13) is given in Appendix 5.A. Equation (5.12) means that all cases are accepted in the study and it is usually called a "full ascertainment" assumption (see for instance Golan et al. ( 2014)).

Heritability estimator

l = g + e,
where g is a genetic random effect, which can be correlated across individuals, and e is the environmental random effect, which is assumed to be independent of the genetic effect. Both effects are assumed to be Gaussian: e has a variance equal to (1-η )Id R n and g has a covariance matrix, the diagonal entries of which are equal to η and the non diagonal term (i, j) is equal to η G i,j . The covariance matrix of (l i , l j ) is given by

Σ = 1 η G i,j η G i,j
1 .

The heritability estimator proposed by [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF] is a least square estimator obtained by minimizing i =j

(W i W j -E[W i W j | i = j = 1]) 2 .
(5.14)

Since the expression of E[W i W j | i = j = 1] has no explicit formula as we shall see hereafter, [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF] proposed to take advantage of the fact that the correlations G i,j are small for i = j. The ground of the method is to write

E[W i W j | i = j = 1] = 1-P P P(Y i = Y j = 1) -K(1-P ) P (1-K) P(Y i = Y j ) + K 2 (1-P ) P (1-K) 2 P(Y i = Y j = 0) P(Y i = Y j = 1) + K(1-P ) P (1-K) 2 P(Y i = Y j = 0) + K(1-P ) P (1-K) P(Y i = Y j ) (5.15)
and to propose approximations of P(Y i = Y j ), P(Y i = Y j = 0) and P(Y i = Y j ) thanks to Taylor developments around the quantity G i,j . The computations leading to (5.15) can be found in Appendix 5.B. This approximation, plugged in the least square criterion (5.14), led to the heritability estimator given by

η =    i =j W i W j G i,j c i =j G 2 i,j ∧ 1    ∨ 0, (5.16)
where

c = φ(t) 2 P (1 -P ) K 2 (1 -K) 2 ,
(5.17)

φ being the density of the standard Gaussian distribution.

Our method

In Model (5.7) that we consider, the variance matrix Σ (N ) of (l i , l j ) conditionally to Z can be written as

Σ (N ) = 1 + η (G N (i, i) -1) η G N (i, j) η G N (i, j) 1 + η (G N (j, j) -1) ,
where for all 1 ≤ i, j ≤ n,

G N (i, j) = 1 N N k=1 Z i,k Z j,k .
(5.18)

Note that in the model we consider, G N (i, j) is a random variable, which is not the case of the quantity G i,j in the model studied by [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF]. A key element is to notice that Σ (N ) is close to the n × n identity matrix, more precisely

Σ (N ) = 1 + η A N (i) √ N η B N (i,j) √ N η B N (i,j) √ N 1 + η A N (j) √ N (5.19) where A N (i) = O p (1), A N (j) = O p (1) and B N (i, j) = O p (1)
. The proof of (5.19) can be found in Appendix 5.C. Then, following the idea of [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF], we propose to approximate

E[W i W j |Z, i = j = 1]
defined in Equation (5.15) thanks to Taylor developments around

A N (i) √ N , A N (j)
√ N and B N (i,j) √ N . The detailed computations are devised in Section 5.8.1. A first order approximation of E[W i W j |Z, i = j = 1], plugged in (5.14), leads to the same estimator η(1) as the one proposed by [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF]. Indeed, we obtain (5.20) where c = φ(t) 2 P (1-P ) K 2 (1-K) 2 . In Section 5.5, we consider the second order approximation, which is different from the one devised by [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF].

η(1) =    i =j W i W j G N (i, j) c i =j G N (i, j) 2 ∧ 1    ∨ 0,

Consistency of the heritability estimator η(1)

In this section, we consider the heritability estimator η(1) defined in Equation (5.20).

Assumption 2. There exist d > 0, C > 0 and a neighborhood V 0 of 0 such that for all λ in V 0

1.1 E[exp λ(A i,k -E[A i,k ]) 2 -σ 2 k ] ≤ C exp(dλ 2 ) 1.2 E[exp (λ(A i,k -E[A i,k ]))] ≤ C exp(dλ 2 ) 1.3 E[exp (λ(A i,k -E[A i,k ])(A j,k -E[A i,k ]))] ≤ C exp(dλ 2 )
for all i = j and for all k, where the A i,k 's are defined in (5.9) and σ 2 k is the variance of A i,k . Assumption 3.

inf

k=1..N σ 2 k = δ min > 0 2.2 sup k=1..N σ 2 k = δ max < +∞ Theorem 4. Let Y = (Y 1 , .
. . , Y n ) satisfy Model (5.5) with A satisfying Assumptions 2 and 3, and η(1) the estimator of η defined in Equation (5.20). Then, as n, N → ∞ such that n/N → a ∈ (0, +∞), η(1) = η + o p (1).

Second order approximation of E[W

i W j |Z, i = j = 1]
The purpose of this section is to study the behaviour of the heritability estimator obtained thanks to a second order approximation of

E[W i W j |Z, i = j = 1].
Instead of computing the approximation till order 1/ √ N , we compute the approximation till order 1/N and we obtain:

E[W i W j |Z, i = j = 1] = η √ N P (1 -P ) K 2 (1 -K) 2 φ(t) 2 B N (i, j) + t 2 4 η 2 N A N (i)A N (j) P (1 -P ) K 2 (1 -K) 2 + η 2 N P (1 -P ) K 2 (1 -K) 2 φ(t) 2 B N (i, j) 2 t 2 2 - (P -K) 2 K 2 (1 -K) 2 + η 2 2N P (1 -P ) K 2 (1 -K) 2 φ(t) 2 B N (i, j)(A N (i) + A N (j)) t 2 -1 - P -K K(1 -K) tφ(t) + O p N
The proof of this computation is detailed in Section 5.8.3. Since the minimizer in η of the quantity

g(η) = i =j W i W j - η √ N P (1 -P ) K 2 (1 -K) 2 φ(t) 2 B N (i, j) - t 2 4 η 2 N A N (i)A N (j) P (1 -P ) K 2 (1 -K) 2 - η 2 N P (1 -P ) K 2 (1 -K) 2 φ(t) 2 B N (i, j) 2 t 2 2 - (P -K) 2 K 2 (1 -K) 2 - η 2 2N P (1 -P ) K 2 (1 -K) 2 φ(t) 2 B N (i, j)(A N (i) + A N (j)) t 2 -1 - P -K K(1 -K) tφ(t)
2 has no explicit form, we use a Newton-Raphson approach to obtain the corresponding heritability estimator η(2) of the second order approximation. Note that the second order approximation, which depends on B N (i, j) but also on A N (i) and A N (j), is different from the one found by [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF].

Numerical study

In this section, we propose to study the numerical performance of the estimators η(1) and η(2) devised respectively in Sections 5.3 and 5.5. Since [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF] already compared the estimator η(1) to the one proposed by [START_REF] Lee | Estimating missing heritability for disease from genome-wide association studies[END_REF] and stated several arguments in favor of their estimator, we will focus on comparing our two estimators in terms of statistical and computational efficiency.

Simulation process

In this simulation study, we generated data sets with n 200, N = 10000 in order to respect the classical scenario where N >> n. The value of the prevalence in the population varies from 0.005 to 0.1. The observations were generated as follows.

• We set the parameters η , K, P = 1/2 and the size of the general population, chosen very large. Notice that the number of individuals selected in the study varies from one sample to another. We chosed in practice a population size in order to have around 100 patients in the study.

• We generated the Gaussian random effects e and u with respective variances σ 2 u = η /N and σ 2 e = 1 -η .

• We generated liabilities, from which we generated binary observations in order to have a certain prevalence K in the population, that is certain number of cases.

• For each individual, we determined those who stayed in the study: the cases are automatically selected (full ascertainment assumption) but each control is selected with probability p control computed in Equation (5.13).

Results

Figure 5.1 displays the estimations of η obtained with both estimators η(1) and η(2) . First, we can notice that both estimators seem empirically unbiased. Second, we observe no obvious improvement of the performance of η(2) compared to η(1) in terms of empirical variance. Finally, we can also note that the estimations seem more accurate when the prevalence K is high, namely K = 0.1. Table 5.1 and Figure 5.2 show the computational performance of both estimators. The computation of the estimator η(2) obtained with the more refined approximation is obviously slower, but for small values of n (namely, n = 100), the time required to compute an estimation of η remains quite small (86 seconds, against 40 seconds for the other estimator). However, when n is larger, the computational time increase substantially and the "slower" estimator needs up to 13500 seconds, that is almost 4 hours, to compute an estimation of η .

In conclusion, both estimators are empirically unbiased and since the computation of the estimator η(2) is slower and does not improve the estimations of η , we are satisfied with the first order approximation and the corresponding estimator η(1) .

Discussion

In this paper, we proposed theoretical grounds to support the heritability estimator in casecontrol studies developed by [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF]. We proved indeed its consistency in the framework where both the number of individuals n and the number N of SNPs go to infinity, remains a challenging issue. Considering possible sparsity in the random effects would also be an interesting improvement and will be the subject of a future work.

Proofs

5.8.1 Taylor development of E[W i W j |Z, i = j = 1] in Model (5.7)
According to Equation (5.15), we only need to compute approximations of P(Y i = Y j = 1|Z), P(Y i = Y j = 0|Z) and P(Y i = Y j |Z) to obtain an approximation of E(W i W j |Z, i = j = 1).

P(Y i = Y j = 1|Z) = ∞ t ∞ t f (x, y)dxdy, P(Y i = Y j = 0|Z) = t -∞ t -∞ f (x, y)dxdy and
P(Y i = Y j |Z) = 2 t -∞ ∞ t f (x, y)dxdy, with f (x, y) = 1 2π |Σ (N ) | -1 2 exp - (x, y)Σ (N )-1 (x, y) t 2 .
where the matrix Σ (N ) is the covariance matrix of (l i , l j ).

We will use the result of Equation (5.19), which will be demonstrated in Appendix 5.C, that is

Σ (N ) = 1 + η A N (i) √ N η B N (i,j) √ N η B N (i,j) √ N 1 + η A N (j) √ N , (5.22)
where

A N (i) = O p (1), A N (j) = O p (1) and B N (i, j) = O p (1). We have f (x, y) = 1 2π|Σ (N ) | -1 2 exp - 1 2|Σ (N ) | x 2 (1 + η √ N A N (j)) + y 2 (1 + η √ N A N (i)) -2xy η √ N B N (i, j) = 1 2π|Σ (N ) | -1 2 exp(- x 2 2 ) exp(- y 2 2 ) exp - x 2 2 1 |Σ (N ) | 1 + η √ N A N (j) -1 - y 2 2 1 |Σ (N ) | 1 + η √ N A N (i) -1 + 1 |Σ (N ) | xy η √ N B N (i, j) .
Using a first order Taylor development,

|Σ (N ) | -1 = 1 -(A N (i) + A N (j)) η √ N + α N and |Σ (N ) | -1 2 = 1 - 1 2 (A N (i) + A N (j)) η √ N + β N ,
where

α N = O p ( 1 N ) and β N = O p ( 1 N ). More precisely, α N = -(A N (i)A N (j) -B N (i, j) 2 ) η 2 N + 1 2 -(A N (i) + A N (j)) η √ N -(A N (i)A N (j) -B N (i, j) 2 ) η 2 N 2 1 (1 + α) 3 , with |α| ≤ |(A N (i) + A N (j)) η √ N + (A N (i)A N (j) -B N (i, j) 2 ) η 2 N |. Similarly, β N = - 1 2 (A N (i)A N (j) -B N (i, j) 2 ) η 2 N + 1 2 - 1 2 (A N (i) + A N (j)) η √ N - 1 2 (A N (i)A N (j) -B N (i, j) 2 ) η 2 N 2 3 4 1 (1 + β) 5 2 , with | β| ≤ | 1 2 (A N (i) + A N (j)) η √ N + 1 2 (A N (i)A N (j) -B N (i, j) 2 ) η 2 N |.
Then,

f (x, y) = 1 - 1 2 (A N (i) + A N (j)) η √ N + β N φ(x)φ(y) × exp - x 2 2 (-A N (i) η √ N + γ N ) - y 2 2 (-A N (j) η √ N + γN ) + xy η √ N B N (i, j) + γN where γ N = -A N (j)(A N (i) + A N (j)) η 2 N + α N (1 + A N (j) η √ N ) = O p 1 N , γN = -A N (i)(A N (i) + A N (j)) η 2 N + α N (1 + A N (i) η √ N ) = O p 1 N and γN = η √ N B N (i, j) -(A N (i) + A N (j)) η √ N + α N = O p 1 N
A Taylor development of the exponential function leads to

f (x, y) = 1 - 1 2 (A N (i) + A N (j)) η √ N + β N φ(x)φ(y) × 1 + x 2 2 η √ N A N (i) + y 2 2 η √ N A N (j) + xy η √ N B N (i, j) + ν N (x) with ν N (x) = - x 2 2 γ N - y 2 2 γN + xy γN + 1 2 x 2 2 η √ N A N (i) + y 2 2 η √ N A N (j) + xy η √ N B N (i, j) - x 2 2 γ N - y 2 2 γN + xy γN 2 exp ũ where |ũ| ≤ | x 2 2 η √ N A N (i) + y 2 2 η √ N A N (j) + xy η √ N B N (i, j) -x 2 2 γ N -y 2 2 γN + xy γN |. Then, ∞ t ∞ t f (x, y)dxdy = 1 - 1 2 (A N (i) + A N (j)) η √ N + β N K 2 + 1 2 η √ N (A N (j) + A N (i))K(K + tφ(t)) + B N (i, j) η √ N φ(t) 2 + µ N = K 2 + 1 2 (A N (i) + A N (j)) η √ N Ktφ(t) + B N (i, j) η √ N φ(t) 2 + µ N where µ N = 1 -1 2 (A N (i) + A N (j) + β N ) η √ N ∞ t ∞ t φ(x)φ(y)ν N (x)dxdy and µ N = µ N + β N K 2 + 1 2 η √ N (A N (j) + A N (i))K(K + tφ(t)) + B N (i, j) η √ N φ(t) 2 - 1 2 (A N (i) + A N (j)) η 2 N B N (i, j)φ(t) 2 .
This remainder and its order will be carefully studied in Section 5.8.2.

Chapter 5 -Heritability estimation in case-control studies Similarly, we can compute P(Y i = Y j = 0|Z) and

P(Y i = Y j |Z): t -∞ t -∞ f (x, y)dxdy = (1 -K) 2 - 1 2 (A N (i) + A N (j)) η √ N (1 -K)tφ(t) + B N (i, j) η √ N φ(t) 2 + μN t -∞ ∞ t f (x, y)dxdy + ∞ t t -∞ f (x, y)dxdy = 2K(1 -K) + (A N (i) + A N (j)) η √ N (1 -2K)tφ(t) -2B N (i, j) η √ N φ(t) 2 + μN .
Replacing these terms in the expression of the numerator of E(W i W j |Z, i = j = 1) given in equation ( 5.15) leads to:

η √ N B N (i, j)φ(t) 2 (1 -P ) P (1 -K) 2 + r N , (5.23)
where r N is a linear combination of µ N , μN and μN . Since there is no constant term in this numerator, we only need the development of order 0 of the denominator of E(W i W j |Z, i = j = 1) to obtain the first order approximation of E(W i W j |Z, i = j = 1).

We obtain that the denominator can be written as

K 2 P 2 + rN ,
where rN is the sum of a term of order 1 √ N and a linear combination of µ N , μN and μN . Thus, we obtain that

E(W i W j |Z, i = j = 1) = η √ N B N (i, j)φ(t) 2 (1-P ) P (1-K) 2 + r N K 2 P 2 + rN (5.24) = η G N (i, j)φ(t) 2 P (1 -P ) K 2 (1 -K) 2 + R N (i, j) (5.25)
where

R N (i, j) = η √ N B N (i, j)φ(t) 2 (1 -P ) P (1 -K) 2 + r N rN + K 2 P 2 r N .
(5.26)

Proof of Theorem 4

Properties of Z

In the following proofs, we will use several properties of the matrix Z, which are stated in Proposition 1.

Proposition 1. Uniformly in k,

(1)

E(Z 1,k Z 2,k ) = -1 n-1 .
(2) E(Z p 1,k ) = O(1), for all p.

(3

) E(Z 2 1,k Z 2 2,k ) = 1 + o(1). (4) E[Z 3 1,k Z 2,k ] = O 1 n . (5) E[Z 2 1,k Z 2,k Z 3,k ] = O 1 n . (6) E[Z 1,k Z 2,k Z 3,k Z 4,k ] = O 1 n 2 . (7) E[Z 5 1,k Z 2,k ] = O 1 n . (8) E[Z 3 1,k Z 3 2,k ] = O(1). (9) E[Z 4 1,k Z 2 2,k ] = O(1). (10) E[Z 4 1,k Z 2,k Z 3,k ] = O( 1 n ). (11) E[Z 3 1,k Z 2 2,k Z 3,k ] = O( 1 n ). (12) E[Z 3 1,k Z 2,k Z 3,k Z 4,k ] = O( 1 n 2 ).
The proof of Proposition 1 is given in Appendix 5.D.

Proof of Lemma 7

Let us prove that, when n and N go to infinity and n/N goes to a,

1 n i =j G N (i, j) 2 P → a,
where P → denotes the convergence in probability.

G N (i, j) 2 = 1 N 2 N k=1 Z 2 i,k Z 2 j,k + 1 N 2 k =l Z i,k Z j,k Z i,l Z j,l
Since Z i,k and Z j,l are independent for any i and j when k = l, we will always consider separately the cases where k = l from the cases where k = l. Indeed, let us show that

1 n 1 N 2 i =j N k=1 Z 2 i,k Z 2 j,k P → a (5.27) and 1 n 1 N 2 i =j k =l Z i,k Z j,k Z i,l Z j,l P → 0.
(5.28)

Note that

E( 1 n 1 N 2 i =j N k=1 Z 2 i,k Z 2 j,k ) = 1 n 1 N 2 i =j N k=1 E(Z 2 i,k Z 2 j,k ) = n -1 N (1 + o(1)) by (3) of Proposition 1 = a + o(1) Moreover, Var( 1 n 1 N 2 i =j N k=1 Z 2 i,k Z 2 j,k ) = 1 n 2 1 N 4 N k=1 i 1 =j 1 i 2 =j 2 E(Z 2 i 1 ,k Z 2 j 1 ,k Z 2 i 2 ,k Z 2 j 2 ,k ) - 1 n 2 1 N 4 N k=1   i =j E(Z 2 i,k Z 2 j,k )   2
(5.29)

The second term of (5.29) can be rewritten as:

1 n 2 1 N 4 N k=1   i =j E(Z 2 i,k Z 2 j,k )   2 = N n 2 (n -1) 2 n 2 N 4 (1 + o(1)) by (3) of Proposition 1 = O 1 n i 1 =j 1 i 2 =j 2 E(Z 2 i 1 ,k Z 2 j 1 ,k Z 2 i 2 ,k Z 2 j 2 ,k ) ≤ E( i 1 ,j 1 ,i 2 ,j 2 Z 2 i 1 ,k Z 2 j 1 ,k Z 2 i 2 ,k Z 2 j 2 ,k ) = E n i=1 Z 2 i,k 4 = n 4
This last equality comes from the definition of Z as a centered and normalized variable given in Equation (5.9), which implies that for all k,

n i=1 Z 2 i,k = n.
Then,

1 n 2 1 N 4 N k=1 i 1 =j 1 i 2 =j 2 E(Z 2 i 1 ,k Z 2 j 1 ,k Z 2 i 2 ,k Z 2 j 2 ,k ) ≤ n 4 N n 2 N 4 = O 1 n .
This proves (5.27).

E(

1 n 1 N 2 i =j k =l Z i,k Z j,k Z i,l Z j,l ) = 1 n 1 N 2 i =j k =l E(Z i,k Z j,k )E(Z i,l Z j,l ) = n(n -1)N (N -1) nN 2 (n -1) 2 by (1) of Proposition 1 = O 1 n Var( 1 n 1 N 2 i =j k =l Z i,k Z j,k Z i,l Z j,l ) = 1 n 2 1 N 4 k =l i 1 =j 1 i 2 =j 2 E(Z i 1 ,k Z i 2 ,k Z j 1 ,k Z j 2 ,k )E(Z i 1 ,l Z i 2 ,l Z j 1 ,l Z j 2 ,l ) - 1 n 2 1 N 4 k =l   i =j E(Z i,k Z j,k )E(Z i,l Z j,l )   2 1 n 2 1 N 4 k =l   i =j E(Z i,k Z j,k )E(Z i,l Z j,l )   2 = N (N -1)n 2 (n -1) 2 n 2 N 4 (n -1) 4 by (1) of Proposition 1 = O 1 n 4
In the first term, {i 1 , i 2 , j 1 , j 2 } can be of cardinal 2, 3 or 4 and counting the number of combinations gives the expression:

i 1 =j 1 i 2 =j 2 E(Z i 1 ,k Z i 2 ,k Z j 1 ,k Z j 2 ,k )E(Z i 1 ,l Z i 2 ,l Z j 1 ,l Z j 2 ,l ) = 2 i =j E(Z 2 i,k Z 2 j,k )E(Z 2 i,l Z 2 j,l ) + 4 i =j 1 =j 2 E(Z 2 i,k Z j 1 ,k Z j 2 ,k )E(Z 2 i,l Z j 1 ,l Z j 2 ,l ) + i 1 =i 2 =j 1 =j 2 E(Z i 1 ,k Z i 2 ,k Z j 1 ,k Z j 2 ,k )E(Z i 1 ,l Z i 2 ,l Z j 1 ,l Z j 2 ,l ) = 2n(n -1)(1 + o(1)) + 4 n(n -1)(n -2) n o(1) + n(n -1)(n -2)(n -3) n 2 o(1) = O(n 2 )
This was obtained by using (3),( 5) and ( 6) of Proposition 1. Finally,

V ar   1 n 1 N 2 i =j k =l Z i,k Z j,k Z i,l Z j,l   = O 1 n 2 .
This completes the proof of (5.28).

with

ν N (x, y) = - x 2 2 γ N - y 2 2 γN + xy γN + 1 2 x 2 2 η √ N A N (i) + y 2 2 η √ N A N (j) + xy η √ N B N (i, j) - x 2 2 γ N - y 2 2 γN + xy γN 2 exp ũ
Integrating the first terms of ν N (x, y) gives

∞ t ∞ t φ(x)φ(y) - x 2 2 γ N - y 2 2 γN + xy γN dxdy = - 1 2 K(tφ(t) + K)(γ N + γN ) + φ(t) 2 γN .
. Moreover, we have the upper bound exp ũ ≤ max(exp

x 2 2 η √ N A N (i) + y 2 2 η √ N A N (j) + xy η √ N B N (i, j) + x 2 2 γ N - y 2 2 γN -xy γN , 1). If max(exp x 2 2 η √ N A N (i) + y 2 2 η √ N A N (j) + xy η √ N B N (i, j) -x 2 2 γ N -y 2 2 γN + xy γN , 1) = 1, ∞ t ∞ t φ(x)φ(y)ν N (x, y)dxdy ≤ ∞ t ∞ t φ(x)φ(y) x 2 2 η √ N A N (i) + y 2 2 η √ N A N (j) + xy η √ N B N (i, j) - x 2 2 γ N - y 2 2 γN + xy γN 2 dxdy = 1 N J, (5.39) where J = ∞ t ∞ t φ(x)φ(y) x 2 2 η A N (i) + y 2 2 η A N (j) + xyη B N (i, j) -x 2 2 γ N √ N -y 2 2 γN √ N + xy γN √ N 2 dxdy is finite. Otherwise, exp(ũ) ≤ exp x 2 2 ( N η + P 1 ( N )) + y 2 2 ( N η + P 2 ( N )) + xy( N η + P 3 ( N )) ,
where P 1 , P 2 , P 3 are polynomial functions. This expression comes from upper bounding the terms A N (i)/N , A N (j)/N and B N (i, j)/N by N in γ N , γN and γN . There exists N 0 , such that for all N ≥ N 0 , N η + P 1 ( N ) ≤ 1 4 , N η + P 2 ( N ) ≤ 1 4 and N η + P 3 ( N ) ≤ 1 4 .

. Thus, similarly to the expression 5.39,

∞ t ∞ t φ(x)φ(y) x 2 2 η √ N A N (i) + y 2 2 η √ N A N (j) + xy η √ N B N (i, j) + x 2 2 γ N - y 2 2 γN -xy γN 2 exp(ũ)dxdy ≤ 1 2π ∞ t ∞ t exp(- x 2 4 ) exp(- y 2 4 ) x 2 2 η √ N A N (i) + y 2 2 η √ N A N (j) + xy η √ N B N (i, j) + x 2 2 γ N - y 2 2 γN -xy γN 2 dxdy ≤ 1 N J
where J is finite.

Similarly to the computations made for α N , β N , γ N , µ N , all the remainder terms can be upper bounded by products of A N (i)/ √ N , A N (j)/ √ N and B N (i, j)/ √ N , which proves (5.38).

Proof of Lemma 10

In this section, all the expectations that we consider are conditionally to the presence of the observed individuals in the study, for instance

{ i = j = 1} or { i 1 = i 2 = i 3 = 1} .
However, for the sake of simplicity, we will not always make explicit such conditioning.

Let us show that

Var( 1 n i =j

(W i W j -E[W i W j |Z])G N (i, j)1 E N ) → 0, that is 1 n 2 i 1 =i 2 i 3 =i 4 E [(E[W i 1 W i 2 W i 3 W i4 |Z] -E[W i 1 W i 2 |Z]E[W i 3 W i 4 |Z])G N (i 1 , i 2 )G N (i 3 , i 4 )1 E N ] → 0 (5.40)
For this purpose, we will separate three cases depending on the cardinal of the set {i 1 , i 2 , i 3 , i 4 } in the sum of Equation (5.40).

-If card({i 1 , i 2 , i 3 , i 4 })=2, the corresponding terms in (5.40) are equal to

1 n 2 i =j E E[(W 2 i W 2 j |Z] -E[W i W j |Z] 2 )G N (i, j) 2 1 E N ≤ 1 n 2 i =j E[(α + ρ N (i, j))G N (i, j) 2 1 E N ]
where α is a positive constant and ρ N (i, j) can be upper bounded by a finite product of G N (i, j), G N (i, i) -1 and G N (j, j) -1, according to proof of Lemma 9. This result is obtained by using a similar decomposition of E[W 2 i W 2 j |Z] than the one that we explicited for

E[W i W j |Z]. Since E[G N (i, j) 2 1 E N ] ≤ 2
N and all terms of ρ N (i, j) are upper bounded by a finite sum of k N , with k greater than 1, which all tend to 0, it is clear that

1 n 2 i =j E E[(W 2 i W 2 j |Z] -E[W i W j |Z] 2 )G N (i, j) 2 1 E N → 0.
-If card({i 1 , i 2 , i 3 , i 4 })=3, the corresponding terms in (5.40) are equal to

1 n 2 i 1 =i 2 =i 3 E (E[W 2 i 1 W i 2 W i 3 |Z] -E[W i 1 W i 2 |Z]E[W i 1 W i 3 |Z])G N (i 1 , i 2 )G N (i 1 , i 3 )1 E N .
(5.41)

Since the sum of Equation ( 5.41) has n(n -1)(n -2) terms, we have the refine the upper bound that we used in the case where the cardinal of {i 1 , i 2 , i 3 , i 4 } was equal to 2. Indeed, we will use the following proposition:

Proposition 2. E[W 2 i 1 W i 2 W i 3 |Z] has no term of order less than 1/ √ N , that is no constant term.
Let us explain why Proposition 2 is enough to prove

1 n 2 i 1 =i 2 =i 3 E (E[W 2 i 1 W i 2 W i 3 |Z] -E[W i 1 W i 2 |Z]E[W i 1 W i 3 |Z])G N (i 1 , i 2 )G N (i 1 , i 3 )1 E N → 0.
(5.42)

Let us first recall that, according to Lemma 9,

E[W i 1 W i 2 |Z]E[W i 1 W i 3 |Z] = c 2 η 2 G N (i 1 , i 2 )G N (i 1 , i 3 ) + cη G N (i 1 , i 3 )R N (i 1 , i 2 ) + cη G N (i 1 , i 2 )R N (i 1 , i 3 ) + R N (i 1 , i 2 )R N (i 1 , i 3 ),
where, if E N holds, all these terms are upper bounded by a finite sum of terms of the form k N , with k ≥ 2. Then,

E [E[W i 1 W i 2 |Z]E[W i 1 W i 3 |Z]G N (i 1 , i 2 )G N (i 1 , i 3 )1 E N ]
can be upper bounded by a finite sum of terms of the form k N , with k ≥ 4. Since

N (N -1)(N -2) 4 N n 2 → 0, it shows that 1 n 2 i 1 =i 2 =i 3 E [E[W i 1 W i 2 |Z]E[W i 1 W i 3 |Z])1 E N ] → 0.
Similarly, according to Proposition 2, each term of

E E[W 2 i 1 W i 2 W i 3 |Z]G N (i 1 , i 2 )G N (i 1 , i 3 )1 E N can be upper bounded by a finite sum of k N , with k ≥ 3. Since n(n -1)(n -2) 3 N n 2 = O 1 N 1/2-3γ → 0,
it achieves the proof of (5.41).

-If card({i 1 , i 2 , i 3 , i 4 })=4, let us first observe that

N (N -1)(N -2)(N -3) 5 N n 2 → 0,
which means that we shall only focus on the approximation of

E[W i 1 W i 2 W i 3 W i4 |Z] -E[W i 1 W i 2 |Z]E[W i 3 W i 4 |Z] of order 1/N . Let us recall that E [W i 1 W i 2 |Z]E[W i 3 W i 4 |Z] = c 2 η 2 G N (i 1 , i 2 )G N (i 3 , i 4 ) + R N (i 1 , i 2 , i 3 , i 4 ),
where

R N (i 1 , i 2 , i 3 , i 4 ) = cη G N (i 1 , i 2 )R N (i 3 , i 4 ) + cη G N (i 3 , i 4 )R N (i 1 , i 2 ) + R N (i 1 , i 2 )R N (i 3 , i 4 )
is a remainder, each term of which is upper bounded by a finite sum of terms of the form k N , with k ≥ 2. In particular, it implies that

E[ N (N -1)(N -2)(N -3) n 2 R N (i 1 , i 2 , i 3 , i 4 )G N (i 1 , i 2 )G N (i 3 , i 4 )] → 0.
Thus, we need to prove that (5.43) To do so, we shall prove first the following proposition:

1 n 2 i 1 =i 2 =i 3 E (E[W i 1 W i 2 W i 3 W i 4 |Z] -c 2 η 2 G N (i 1 , i 2 )G N (i 3 , i 4 ))G N (i 1 , i 2 )G N (i 3 , i 4 )1 E N → 0,
Proposition 3. The terms of order less than or equal to 1/

√ N in E[W i 1 W i 2 W i 3 W i4 |Z] are null. Σ =     1 + η (G N (i 1 , i 1 ) -1) G N (i 1 , i 2 ) G N (i 1 , i 3 ) G N (i 1 , i 4 ) G N (i 1 , i 2 ) 1 + η (G N (i 2 , i 2 ) -1) G N (i 2 , i 3 ) G N (i 2 , i 4 ) G N (i 1 , i 3 ) G N (i 2 , i 3 ) 1 + η (G N (i 3 , i 3 ) -1) G N (i 3 , i 4 ) G N (i 1 , i 4 ) G N (i 2 , i 4 ) G N (i 3 , i 4 ) 1 + η (G N (i 4 , i 4 ) -1)     .
For the sake of clarity, let us denote

A 1 = 1 √ N N k=1 (Z 2 i 1 ,k -1) = √ N (G N (i 1 , i 1 ) -1)
, and similarly we define A 2 , A 3 and A 4 .

Let us also denote C 1,2 = √ N G N (i 1 , i 2 ) and similarly, C 1,3 , . . . , C 3,4 . Then, let us rewrite Σ as:

Σ =       1 + η √ N A 1 C 1,2 √ n C 1,3 √ N C 1,4 √ N C 1,2 √ N 1 + η √ N A 2 C 2,3 √ N C 2,4 √ N C 1,3 √ N C 2,3 √ N 1 + η √ N A 3 C 3,4 √ N C 1,4 √ N C 2,4 √ N C 3,4 √ N 1 + η √ N A 4      
.

The approximation of order 1/ √ n of its inverse matrix is given by

Σ -1 |Σ| -1 ×       1 + η √ N (A 2 + A 3 + A 4 ) - C 1,2 √ N - C 1,3 √ N - C 1,4 √ N - C 1,2 √ N 1 + η √ N (A 1 + A 3 + A 4 ) - C 2,3 √ N - C 2,4 √ N - C 1,3 √ N - C 2,3 √ N 1 + η √ N (A 1 + A 2 + A 4 ) - C 3,4 √ N - C 1,4 √ N - C 2,4 √ N - C 3,4 √ N 1 + η √ n (A 1 + A 2 + A 2 )      
where

|Σ| -1 = 1 -η √ N (A 1 + A 2 + A 3 + A 4 ) + O p 1 N . Let us compute P(Y i 1 = Y i 2 = Y i 3 = Y i 4 = 1|Z) = +∞ t +∞ t +∞ t +∞ t f (w, x, y, z)dwdxdydz, where f (w, x, y, z) = 1 (2π) 2 |Σ| 1 2 exp - w 2 2|Σ| (1 + η √ n (A 2 + A 3 + A 4 )) -• • • - z 2 2|Σ| (1 + η √ N (A 1 + A 2 + A 3 )) + wx |Σ| η √ n C 1,2 + wy |Σ| η √ n C 1,3 + • • • + yz |Σ| η √ N C 3,4 + O p 1 N = 1 (2π) 2 |Σ| 1 2 exp - w 2 2 (1 - η √ N (A 1 + A 2 + A 3 + A 4 ))(1 + η √ N (A 2 + A 3 + A 4 )) -. . . - z 2 2 (1 - η √ N (A 1 + A 2 + A 3 + A 4 ))(1 + η √ n (A 1 + A 2 + A 3 )) +wx η √ N (1 - η √ N (A 1 + A 2 + A 3 + A 4 ))C 1,2 + . . . +yz η √ N (1 - η √ N (A 1 + A 2 + A 3 + A 4 ))C 3,4 + O p 1 N = 1 |Σ| 1 2 φ(w)φ(x)φ(y)φ(z) exp w 2 2 η √ N A 1 + • • • + z 2 2 η √ N A 4 -wx η √ N C 1,2 -• • • -yz η √ N C 3,4 + O p 1 N = 1 |Σ| 1 2 φ(w)φ(x)φ(y)φ(z) 1 + w 2 2 η √ N A 1 + • • • + z 2 2 η √ N A 4 -wx η √ N C 1,2 -. . . -yz η √ N C 3,4 + O p 1 N (5.46)
Finally,

• P(Y i 1 = Y i 2 = Y i 3 = Y i 4 = 1|Z) = 1 |Σ| 1 2 K 4 + K 3 2 (tφ(t) + K) η √ n (A 1 + A 2 + A 3 + A 4 ) +K 2 φ(t) 2 η √ n (C 1,2 + • • • + C 3,4 ) + O p 1 N
Similarly, we compute

• P("1 control, 3 cases") = 1 |Σ| 1 2 4K 3 (1 -K) + K 2 2 ((3 -4K)tφ(t) + 4K(1 -K)) η √ N (A 1 + A 2 + A 3 + A 4 ) +2φ(t) 2 K(1 -2K) η √ N (C 1,2 + • • • + C 3,4 ) + O p 1 N 1 |Σ| 1 2 η √ N (A 1 + A 2 + A 3 + A 4 ) (1 -P ) 2 P 2 K 3 2 (tφ(t) + K) - (1 -P ) P K(1 -P ) P (1 -K) K 2 2 ((3 -4K)tφ(t) + 4K(1 -K)) + K(1 -P ) P (1 -K) 2 3K(1 -K) 2 ((1 -2K)tφ(t) + 2K(1 -K)) - P 1 -P K(1 -P ) P (1 -K) 3 (1 -K) 2 2 ((1 -4K)tφ(t) + 4K(1 -K)) + P 2 (1 -P ) 2 K(1 -P ) P (1 -K) 4 (1 -K) 3 2 (-tφ(t) + 1 -K) = 1 (2π) 2 |Σ| 1 2 (1 -P ) 2 K 4 2P 2 [1 -4 + 6 -4 + 1] + 1 (2π) 2 |Σ| 1 2 (1 -P ) 2 K 3 2P 2 (1 -K) [1 -K -3 + 4K + 3(1 -2K) -1 + 4K -K] = 0
Finally, we regroup all the terms in η

√ n (C 1,2 + • • • + C 3,4 ): 1 |Σ| 1 2 (1 -P ) 2 K 2 P 2 (1 -K) 2 φ(t) 2 (1 -K) 2 -2(1 -K)(1 -2K) + 6K 2 -6K + 1 + 2K(1 -2K) + K 2 = 0
This proves Proposition 3.

Let us prove Proposition 5. Let us denote f 2 (w, x, y, z) the density function defined in (5.46) developed till order 1/N .

f 2 (w, x, y, z) = 1 |Σ| 1 2 φ(w)φ(x)φ(y)φ(z) 1 + w 2 2 ( η √ N A 1 - η 2 N (A 2 1 + C 2 1,2 + C 2 1,3 + C 2 1,2 ) + • • • + z 2 2 ( η √ N A 4 - η 2 N (A 2 4 + C 2 1,4 + C 2 2,4 + C 2 3,4 ) + wx(C 1,2 η √ N - η 2 N [(A 1 + A 2 )C 1,2 + C 1,3 C 2,3 + C 1,4 C 2,4 ]) + . . . + yz(C 3,4 η √ N - η 2 N [(A 3 + A 4 )C 3,4 + C 1,3 C 1,4 + C 2,3 C 2,4 ]) + w 4 8 η 2 N A 2 1 + • • • + z 4 8 η 2 N A 2 4 + w 2 x 2 2 η 2 N (C 2 1,2 + A 1 A 2 2 ) + . . . + y 2 z 2 2 η 2 N (C 2 3,4 + A 3 A 4 2 ) + w 3 x 2 η 2 N A 1 C 1,2 + . . . + z 3 y 2 η 2 N A 4 C 3,4 + w 2 xy η 2 N [ A 1 C 2,3 2 + C 1,2 C 1,3 ] + . . . +z 2 xy η 2 N [ A 4 C 2,3 2 + C 2,4 C 3,4 ] + wxyz η 2 N (C 1,2 C 3,4 + C 2,3 C 1,4 + C 1,3 C 2,4 ) (5.47)
In order to prove Proposition 5, we will show that:

1 n 2 i 1 =i 2 =i 3 =i 4 E(A 2 1 C 1,2 C 3,4 ) → 0 (5.48) 1 n 2 i 1 =i 2 =i 3 =i 4 E(A 1 A 2 C 1,2 C 3,4 ) → 0 (5.49) 1 n 2 i 1 =i 2 =i 3 =i 4 E(A 1 C 2 1,2 C 3,4 ) → 0 (5.50) 1 n 2 i 1 =i 2 =i 3 =i 4 E(A 1 C 1,2 C 13 C 3,4 ) → 0 (5.51) 1 n 2 i 1 =i 2 =i 3 =i 4 E(C 2 1,2 C 2,3 C 3,4 ) → 0 (5.52) 1 n 2 i 1 =i 2 =i 3 =i 4 E(C 1,2 C 1,3 C 2,4 C 3,4 ) → 0 (5.53) 1 n 2 i 1 =i 2 =i 3 =i 4 E(C 3 1,2 C 3,4 ) → 0 (5.54)
We will develop the proof of Equation (5.49). By exchangeability of the (Z i,k ) 1≤i≤n , we can write

E[A 1 A 2 C 1,2 C 3,4 ] = k,l,m,r E[(Z 2 1,k -1)(Z 2 2,l -1)Z 1,m Z 2,m Z 3,r Z 4,r ] = k,l,m,r E[Z 2 1,k Z 2 2,l Z 1,m Z 2,m Z 3,r Z 4,r ] -2N k,m,r E[Z 2 1,k Z 1,m Z 2,m Z 3,r Z 4,r ] + N 2 m,r E[Z 1,m Z 2,m Z 3,r Z 4,r ]
(5.55)

We recall that since Z i,k and Z j,l are independent for any i and j when k = l, we will always consider separately the cases where k = l from the cases k = l. Let us first focus on the last term of (5.55).

m,r E[Z 1,m Z 2,m Z 3,r Z 4,r ] = N m=1 E[Z 1,m Z 2,m Z 3,m Z 4,m ] + m =r E[Z 1,m Z 2,m ]E[Z 3,r Z 4,r ] = N × O 1 n 2 + N (N -1) × 1 (n -1) 2
Then,

1 n 2 1 N 4 i 1 =i 2 =i 3 =i 4 (N 2 m,r E[Z 1,m Z 2,m Z 3,r Z 4,r ]) = (N -1)(n-2)(n-3) N n(n-1) + o(1)
Now let us decompose the second term of (5.55) as:

k,m,r E[Z 2 1,k Z 1,m Z 2,m Z 3,r Z 4,r ] = N k=1 E[Z 3 1,k Z 2,k Z 3,k Z 4,k ] + k =l E[Z 3 1,k Z 2,k ]E[Z 3,l Z 4,l ] + k =l E[Z 2 1,k ]E[Z 1,l Z 2,l Z 3,l Z 4,l ] + k =l E[Z 2 1,k Z 3,k Z 3,k ]E[Z 1,l Z 2,l ] + k =l =m E[Z 2 1,k ]E[Z 1,l Z 2,l ]E[Z 3,m Z 4,m ]
Using the results given by Proposition 1, we obtain that

1 n 2 1 N 4   -2N k,m,r E[Z 2 1,k Z 1,m Z 2,m Z 3,r Z 4,r ]   = - 2(N -1)(n -2)(n -3) N n(n -1) + o(1).
Similarly, we can prove that

1 n 2 1 N 4   k,l,m,r E[Z 2 1,k Z 2 2,l Z 1,m Z 2,m Z 3,r Z 4,r ]   = (N -1)(n -2)(n -3) N n(n -1) + o(1),
by using the properties of Proposition 1 or similar relationships coming from other properties of Z that we have not detailed here. Hence we have shown (5.49). The proofs of (5.48), (5.50), (5.51), (5.52), (5.53), (5.54) are very and n i=1

Z 2 i,k = n.

(5.59)

Moreover, for each k, the random variables (Z i,k ) 1≤i≤n are exchangeable. Thus, we deduce from (5.59) that for all i = 1, . . . , n and k = 1, . . . , N , E(Z 2 i,k ) = 1. Hence, by (5.58), we get that 0

= n i=1 Z i,k 2 = n i=1 Z 2 i,k + 1≤i =j≤n Z i,k Z j,k ,
which, by (5.59), implies that for all k = 1, . . . , N and i = j = 1, . . . , n,

E(Z i,k Z j,k ) = - n n(n -1) = - 1 n -1 , (5.60) that is (1).
The proof of (2) comes from the decomposition:

|Z 1,k | p = |Z 1,k | p 1 {s 2 k > δ min 2 } + |Z 1,k | p 1 {s 2 k ≤ δ min 2 } ≤ |A 1,k -Āk | p δ min 2 p + n p 1 {s 2 k ≤ δ min 2 } Assumption 1.2 implies that sup k E |A 1,k -Āk | p <
+∞ and the upper bound for P(s 2 k ≤ δ) of Equation (5.34) prove (2). By (5.59), for all k = 1, . . . , N ,

n 2 = n i=1 Z 2 i,k 2 = n i=1 Z 4 i,k + 1≤i =j≤n Z 2 i,k Z 2 j,k
Since the (Z i,k ) 1≤i≤n are exchangeable for each k = 1, . . . , N , we get that for all k = 1, . . . , N ,

n = E[Z 4 1,k ] + (n -1)E[Z 2 1,k Z 2 2,k ] ,
which gives us (3) by using (2).

If we take the expectation of

Z 3 1,k n i=1 Z i,k = 0, we obtain E[Z 4 1,k ] + (n -1)E[Z 3 1,k Z 2,k ] = 0. Then, (2) implies (4). Similarly, since Z 1,k Z 2,k n i=1 Z 2 i,k = nZ 1,k Z 2,k , we obtain that 2E[Z 3 1,k Z 2,k ] + (n -2)E[Z 2 1,k Z 2,k Z 3,k ] = nE[Z 1,k Z 2,k ].
Then (1) and ( 4) imply (5). Since

Z 1,k Z 2,k Z 3,k n i=1 Z i,k = 0, we obtain that 3E[Z 2 1,k Z 2,k Z 3,k ] + (n -3)E[Z 1,k Z 2,k Z 3,k Z 4,k ] = 0.
Then, ( 5) implies (6). Since

Z 5 1,k n i=1 Z i,k = 0,
we obtain that Then, (2) implies ( 7). The proof of ( 8) is very similar to the proof of (2) but we use Assumption 1.3 which gives us that sup

k E |(A 1,k -Āk )(A 2,k -Āk )| p < +∞. Since Z 4 1,k n i=1 Z 2 i,k = nZ 4 1,k , E[Z 6 1,k ] + (n -1)E[Z 4 1,k Z 2 2,k ] = nE[Z 4 1,k ].
Then (2) implies (9). Similarly, since 7) and ( 9) imply (10). Since

Z 4 1,k Z 2,k n i=1 Z i,k = 0, we obtain that E[Z 5 1,k Z 2,k ] + E[Z 4 1,k Z 2 2,k ] + (n -2)E[Z 4 1,k Z 2,k Z 3,k ] = 0. Then, (
Z 3 1,k Z 2,k n i=1 Z 2 i,k = nZ 3 1,k Z 2,k , we obtain that E[Z 5 1,k Z 2,k ] + E[Z 3 1,k Z 3 2,k ] + (n -2)E[Z 3 1,k Z 2 2,k Z 3,k ] = nE[Z 3 1,k Z 2,k ].
Chapter 5 -Heritability estimation in case-control studies

Then, ( 7), ( 8) and ( 4) imply (11).

Finally, since Z

3 1,k Z 2,k ( n i=1 Z i,k ) 2 = 0, E[Z 5 1,k Z 2,k ] + E[Z 3 1,k Z 3 2,k ] + 2E[Z 4 1,k Z 2 2,k ] + 2(n -2)E[Z 4 1,k Z 2,k Z 3,k ] + 2(n -2)E[Z 3 1,k Z 2 2,k Z 3,k ] + (n -2) 2 E[Z 3 1,k Z 2,k Z 3,k Z 4,k ] = 0
Then, ( 7), ( 8), ( 9), ( 10) and ( 11) imply ( 12).

Chapitre 6

Résumé en français Plusieurs erreurs concernant l'héritabilité sont dues à l'utilisation du terme dans le langage commun, mais avec un sens différent du terme technique utilisé en génétique. Par exemple, une idée reçue serait de penser que l'héritabilité définit la proportion d'un phénotype qui est transmis des parents aux enfants. Tout d'abord, ce ne sont pas les phénotypes qui sont transmis d'une génération à l'autre mais les gènes. De plus, si la moitié des effets génétiques sont effectivement transmis par chaque parent, cette moitié est spécifique à chaque enfant. [START_REF] Visscher | Heritability in the genomics era -concepts and misconceptions[END_REF] a regroupé ces fréquentes erreurs et idées reçues concernant l'héritabilité. Nous allons à présent détailler le concept d'héritabilité tel qu'il est utilisé dans le domaine de la génétique.

Définition de l'héritabilité

En reprenant les explications élégantes de Visscher et al. (2008), nous considérons la modélisation simple où un phénotype d'intérêt est décrit comme le résultat d'effets génétiques et environnementaux, considérés comme indépendants :

Phénotype (P) = Génotype (G) + Environnement (E).
La variance phénotypique (σ 2 P ) peut être décomposée comme la somme des variances non observées (σ 2 G et σ 2 E ) :

σ 2 P = σ 2 G + σ 2 E .
L'héritabilité (H 2 ) est définie mathématiquement comme un rapport de variances et exprime la proportion de la variance phénotypique qui peut être attribuée à des facteurs génétiques :

H 2 = σ 2 G σ 2 P .
La variabilité génétique peut être partitionnée selon différentes sources, en particulier, la variance σ 2 A des effets génétiques additifs. Ces effets additifs sont caractérisés par l'impact des polymorphismes de nucléotide simple (SNPs), qui sont des différences de la séquence d'ADN aux positions du génome auxquelles il existe une variabilité relativement importante dans la population. Ces positions sont en réalité peu fréquentes par rapport à l'intégralité du génome qui contient environ 3 milliards de paires de bases, parmi lesquelles la grande majorité est identique chez tous les êtres humains.

Dans la suite, nous allons considérer l'héritabilité "au sens strict", c'est-à-dire la proportion de variabilité expliquée seulement par les effets génétiques additifs, définie par l'expression suivante :

h 2 = σ 2 A σ 2 P .
Comme l'accès au génotype de milliers d'individus a été rendu possible par la spectaculaire diminution des coûts du séquençage de l'ADN, l'héritabilité de traits quantitatifs ainsi que de pathologies est devenu un sujet très étudié. [START_REF] Yang | Common snps explain a large proportion of the heritability for human height[END_REF] a par exemple estimé que, chez les humains, environ 45 % des variations de la taille sont expliquées par les SNPs les plus fréquents.

L'héritabilité en génétique humaine

La mise en évidence d'une composante génétique importante d'une maladie est également un argument puissant pour réfuter des croyances populaires quant aux causes de certaines maladies. Par exemple, une vague de mouvements anti-vaccins a été déclenchée par un lien présumé entre le vaccin contre l'hépatite B et la sclérose en plaques. De même, le vaccin contre la rougeole a été accusé de provoquer des cas d'autisme [START_REF] Uno | The combined measles, mumps, and rubella vaccines and the total number of vaccines are not associated with development of autism spectrum disorder: the first case-control study in asia[END_REF].

Bien qu'aucun lien n'ait jamais été démontré [START_REF] Poland | Understanding those who do not understand: a brief review of the anti-vaccine movement[END_REF], les conséquences du refus de nombreux parents de vacciner leurs enfants reste un problème de santé publique majeur. En effet, une étude récente [START_REF] Uno | The combined measles, mumps, and rubella vaccines and the total number of vaccines are not associated with development of autism spectrum disorder: the first case-control study in asia[END_REF] a montré que plus de 25 % des parents américains avaient refusé de vacciner leurs enfants contre des maladies mortelles comme la rougeole. Quant à d'autres causes présumées de l'autisme, la théorie de la "mère réfrigérateur" a été formulée par le psychiatre Leo Kanner dans les années 1940 car il avait observé un manque d'affection maternelle chez les mères de ses patients. Bien que cette théorie ait largement été discréditée depuis, les mères de patients autistes ont souffert d'accusations sévères et injustifiées pendant plusieurs décennies.

L'héritabilité en génétique animale et végétale

En génétique animale ou végétale, l'estimation de l'héritabilité est souvent la première étape de la sélection de traits d'intérêt, généralement liés au rendement d'une ressource précieuse. Nous pouvons par exemple citer les exemples d'optimisation de production du lait [START_REF] Visscher | Genetic parameters for milk yield, survival, workability, and type traits for australian dairy cattle[END_REF] ou du blé [START_REF] Eid | Estimation of heritability and genetic advance of yield traits in wheat (Triticum aestivum L.) under drought condition[END_REF]. Le but de Eid ( 2009) est de déterminer des traits héritables liés au rendement du blé et ensuite d'obtenir un génotype optimal. Un tel génotype peut même être choisi pour être le plus résistant possible à des conditions environnementales extrêmes, comme l'absence d'eau, ce qui est actuellement un problème fondamental.

Si ces pratiques sont généralement bien acceptées dans le cadre de la génétique animale, elles ouvrent une controverse quant à de possibles conséquences de l'estimation de l'héritabilité de traits humains. En effet, plusieurs études ont estimé l'héritabilité du QI, par exemple [START_REF] Toro | Genomic architecture of human neuroanatomical diversity[END_REF], et [START_REF] Davies | Genome-wide association studies establish that human intelligence is highly heritable and polygenic[END_REF] a même annoncé que l'intelligence humaine était fortement héritable. La controverse relative à l'héritabilité du QI est discutée dans [START_REF] Visscher | Heritability in the genomics era -concepts and misconceptions[END_REF], qui énumère des raisons à cette polémique. Ces raisons vont de la très controversée définition du QI en tant que mesure de l'intelligence jusqu'à des abus historiques liés à l'eugénisme. Nous ne discuterons pas plus en détail cette polémique ici, mais nous la mentionnant pour illustrer des fréquentes réactions quand on travaille sur l'héritabilité de traits humains.

Après avoir brièvement argumenté en faveur de l'intérêt général d'estimer l'héritabilité, nous allons à présent introduire les modèles statistiques utilisés pour fournir ces estimations.

6.2 Estimation de l'héritabilité dans les modèles linéaires mixtes

Etat de l'art

Les modèles linéaires mixtes (LMMs) sont utilisés dans de nombreux domaines, en particulier en médecine et en génétique. [START_REF] Yang | Common snps explain a large proportion of the heritability for human height[END_REF] ont notamment proposé d'estimer l'héritabilité de la taille humaine en utilisant un modèle linéaire mixte classique défini comme suit :

Y = Xβ + Zu + e (6.1)
où Y = (Y 1 , . . . , Y n ) est le vecteur d'observations d'un phénotype d'intérêt, X est une matrice de taille n × p de prédicteurs, β est un vecteur de taille p × 1 qui contient l'effet inconnu des prédicteurs (aussi appelés effets fixes), et u et e correspondent à des effets aléatoires gaussiens de variances respectives σ 2 u et σ 2 e . De plus, Z est une matrice de taille n × N qui contient l'information génétique. Plus précisément, les Z i,j sont variables aléatoires normalisées au sens où elles sont définies à partir d'une matrice W = (W i,j ) 1≤i≤n, 1≤j≤N par la relation

Z i,j = W i,j -W j s j , i = 1, . . . , n, j = 1, . . . , N , (6.2) où W j = 1 n n i=1 W i,j , s 2 j = 1 n n i=1 (W i,j -W j ) 2 , j = 1, . . . , N . (6.3)
Dans les équations (6.2) et (6.3), les W i,j sont de telle sorte que pour chaque j dans {1, . . . , N } les (W i,j ) 1≤i≤n sont des variables aléatoires indépendantes et identiquement distribuées et les colonnes de W sont indépendantes.

Dans les applications en génétique, la matrice W contient l'information génétique des tous les individus étudiés. Plus précisément, pour chaque j, les (W i,j ) 1≤i≤n sont des variables binomiales i.i.d de paramètres 2 et p j : W i,j = 0 (resp. 1, resp. 2) si le génotype de l'individu i au locus j est qq (resp. Qq, resp. QQ) où p j est la fréquence de l'allèle Q au locus j.

Avec cette définition, les colonnes de Z sont empiriquement centrées et de variance empirique égale à 1.

Le modèle linéaire mixte apparaît comme un choix intuitif de modèle pour décrire le concept d'héritabilité en tant que ratio des variances génétique et phénotypique. [START_REF] Yang | Common snps explain a large proportion of the heritability for human height[END_REF] et [START_REF] Pirinen | Efficient computation with a linear mixed model on large-scale data sets with applications to genetic studies[END_REF] proposent d'estimer le paramètre

η = N σ 2 u N σ 2 u + σ 2 e , (6.4)
communément considéré comme la définition mathématique de l'héritabilité puisqu'il détermine comment les variances sont réparties entre u et e.

Dans le modèle (6.1), la log-vraisemblance conditionnellement à Z est donnée par la formule :

L(β, σ 2 u , σ 2 e ) = - n 2 log(2π) - 1 2 log(|ZZ σ 2 u + σ 2 e Id R n | - 1 2 (Y -Xβ) (ZZ σ 2 u + σ 2 e Id R n ) -1 (Y -Xβ). ( 6
= N σ 2 u + σ 2 e . Plus précisément, Y ∼ N Xβ, η σ 2 R + (1 -η )σ 2 Id R n , où R = ZZ /N .
Si l'on note U la matrice orthogonale (U U = UU = Id R n ) telle que URU = diag(λ 1 , . . . , λ n ) est une matrice diagonale dont les entrées diagonales sont égales à λ 1 , . . . , λ n . Alors Y = U Y est un vecteur gaussien centré de matrice de variance covariance diag(η σ 2 λ 1 + (1η )σ 2 , . . . , η σ 2 λ n + (1 -η )σ 2 ), où les λ i sont les valeurs propres de R. Notons également X = U X. Finalement, [START_REF] Pirinen | Efficient computation with a linear mixed model on large-scale data sets with applications to genetic studies[END_REF] ont maximisé la log-vraisemblance suivante :

L n (β, σ 2 , η) = - n 2 log(σ 2 ) - 1 2 n i=1 log(η(λ i -1) + 1) - 1 2σ 2 n i=1 ( Y i -Xβ) 2 η(λ i -1) + 1 - n 2 log(2π), (6.6) où Y = ( Y 1 , ..., Y n ).
Les méthodes présentées précédemment ont deux possibles faiblesses dans les applications qui nous intéressent. Tout d'abord, elles ont été validées théoriquement uniquement dans le cadre où N est fixé et n tend vers l'infini. En effet, des résultats classiques du modèle linéaire mixte permettent d'obtenir des résultats de consistance et de normalité asymptotiques pour l'estimateur du maximum de vraisemblance. Cependant, puisqu'en pratique le nombre N de SNPs est largement supérieur au nombre n d'individus, il serait plus adapté de valider ces méthodes dans le cadre où n et N tendent vers l'infini, avec n/N qui tend vers une constante a dans (0, +∞).

De plus, toutes ces méthodes ont été développées dans le cas d'effets aléatoires gaussiens, ce qui impliquerait que toute l'information génétique disponible ait un impact sur le phénotype observé. Cette hypothèse, qui semble assez improbable, a été discutée en particulier par [START_REF] Jiang | High-dimensional genome-wide association study and misspecified mixed model analysis[END_REF], qui ont étudié les conséquences d'une modélisation qui ne prend pas en compte le fait que les effets de certains SNPs (possiblement nombreux) soient nuls.

Nous avons pour l'instant mentionné uniquement les méthodes d'estimation de l'héritabilité dans les modèles linéaires mixtes, mais il existe d'autres façons de définir et d'estimer l'héritabilité. En effet, d'importants résultats théoriques ont été prouvés dans le cas du modèle linéaire

Y = Xβ + ε (6.7)
où la composante aléatoire vient d'une part de l'erreur résiduelle ε, supposée gaussienne centrée et de variance σ 2 mais aussi de la matrice de SNPs X dont les colonnes sont supposées gaussiennes centrées réduites indépendantes. Dans ce modèle, l'héritabilité est définie comme le ratio

η = ||β|| 2 2 σ 2 + ||β|| 2 2 . (6.8)
Un avantage de ce modèle est qu'il n'y a pas d'hypothèse sur la distribution de β, en particulier sur sa parcimonie. Néanmoins, des hypothèses fortes sont requises sur la structure de la matrice X. Plusieurs méthodes ont été proposées pour estimer l'héritabilité dans le modèle (6.7). Dicker (2014) a proposé un estimateur de moments qui est asymptotiquement normal quand n, N → a ∈ (0, +∞). [START_REF] Janson | Eigenprism: Inference for high-dimensional signal-tonoise ratios[END_REF] ont développé la procédure Eigenprism pour Chapitre 6 -Résumé en français construire des intervalles de confiance pour η pour des échantillons finis et dont ils ont également étudié les propriétés asymptotiques quand n, N → a ∈ (0, +∞). [START_REF] Dicker | Maximum likelihood for variance estimation in highdimensional linear models[END_REF] ont étudié les propriétés de l'estimateur du maximum de vraisemblance et ont mené une étude de simulation qui compare les méthodes précédemment citées et qui a montré que l'estimateur du maximum de vraisemblance a une variance empirique plus petite que les deux autres approches. [START_REF] Dicker | Maximum likelihood for variance estimation in highdimensional linear models[END_REF] ont également montré la consistance et la normalité asymptotique de l'estimateur du maximum de vraisemblance et ont calculé une formule explicite de la variance asymptotique.

Toujours dans le modèle linéaire, [START_REF] Verzelen | Adaptive estimation of high-dimensional signal-to-noise ratios[END_REF] ont étudié l'optimalité de différentes procédures selon la parcimonie des effets aléatoires. En effet, [START_REF] Verzelen | Adaptive estimation of high-dimensional signal-to-noise ratios[END_REF] ont comparé les performances d'une approche avec sélection (estimateur Gauss-LASSO) et sans sélection (estimateur dense) dans différents régimes de parcimonie. Ils ont distingué des gammes de valeurs pour la parcimonie et pour chacune, ont évalué le risque minimax : ils ont montré que la meilleure procédure était un estimateur adaptatif qui se comporte comme l'estimateur Gauss-LASSO dans les cas très parcimonieux et comme l'estimateur dense dans les cas peu parcimonieux

Contribution

Notre première contribution a été de proposer un estimateur de l'héritabilité dans le contexte où n et N tendent vers l'infini, avec n/N qui tend vers a dans (0, +∞), et d'en établir les propriétés théoriques. Ce travail est développé dans le chapitre 2 de ce manuscrit et a fait l'objet d'un article publié dans Electronic Journal of Statistics. Nous avons étudié le modèle défini par (6.1) sauf que nous avons supposé que les effets aléatoires pouvaient être parcimonieux, c'est-à-dire que seulement une proportion q des composantes de u étaient non nulles :

u i i.i.d. ∼ (1 -q)δ 0 + qN (0, σ u 2 ) , pour tous 1 ≤ i ≤ N et e ∼ N 0, σ e 2 Id R n , (6.9)
où Id R n est la matrice identité de taille n × n, q est dans (0, 1], et δ 0 la masse de Dirac en 0. Quitte à considérer la projection de Y sur l'orthogonal de l'image de X et pour simplifier les calculs, nous avons étudié le modèle suivant Nous avons implémenté un estimateur de η comme le maximiseur de la fonction de vraisemblance dépendant uniquement du paramètre η :

Y = Zu + e . ( 6
L n (η) = -log 1 n n i=1 Y 2 i η(λ i -1) + 1 - 1 n n i=1
log (η(λ i -1) + 1) , (6.11)

Chapitre 6 -Résumé en français

Nous avons obtenu deux résultats principaux dans le cadre où n et N tendent vers l'infini, avec n/N qui tend vers a dans (0, 1) : premièrement, nous avons prouvé que notre estimateur était √ n-consistant malgré la présence de composantes nulles dans les effets aléatoires. Ce résultat a été obtenu sous de faibles hypothèses sur la matrice W et pour toute valeur de parcimonie q.

Nous avons ensuite établi un théorème de la limite centrale sous l'hypothèse supplémentaire que pour tous i et j, les Z i,j étaient des variables gaussiennes centrées réduites et indépendantes. Nous avons calculé une formule explicite de la variance asymptotique, donnée par

τ 2 (a, η , q) = 2 γ 2 (a, η ) + 3 a 2 η 2 γ 4 (a, η ) 1 q -1 S(a, η ) (6.12) où γ 2 (a, η ) = g(η, λ) 2 dµ a (λ) - g(η, λ)dµ a (λ) 2 et S(a, η ) = λ(λ -1) (η (λ -1) + 1) 2 dµ a (λ) - λ (η (λ -1) + 1) dµ a (λ) λ -1 (η (λ -1) + 1) dµ a (λ) 2 .
Dans l'expression précédente, dµ a (λ) est la densité de Marchenko-Pastur, qui est la distribution des valeurs propres de ZZ /N . Cette distribution, obtenue par [START_REF] Marchenko | Distribution of eigenvalues for some sets of random matrices[END_REF] a été un élément clef pour établir les preuves de nos résultats. Nous avons implémenté notre approche dans le package R HiLMM, disponible sur le CRAN.

Nous avons également mené une étude de simulation sur des échantillons finis dont les tailles correspondent à des cas réalistes en pratiques. Nous avons montré que, bien que la variance asymptotique définie par la formule (6.12) dépende en théorie de la parcimonie q, son influence était à peine remarquable en pratique. Par contre, nous avons observé que la variance de notre estimateur était influencée par la paramètre a = n/N . Plus précisément, quand le nombre d'observations est très petit comparé à la taille des effets aléatoires (ce qui est très souvent le cas dans les études en génétique), la variance de l'estimateur était très grande. Ce constat numérique a motivé l'idée de développer une méthode de sélection de variables afin de réduire la taille des effets aléatoires et d'améliorer la précision des estimations de l'héritabilité.

Sélection de variables dans les effets aléatoires des modèles linéaires mixtes

Suite aux performances numériques de notre estimateur de l'héritabilité présentées dans la section précédente, nous avons décidé d'inclure une étape de sélection de variables à notre méthode. Le but de cette étape serait de retrouver le support des effets aléatoires, qui correspond en pratique aux SNPs impliqués dans les variations du phénotype observé. Nous considérerions ensuite la matrice de SNPs réduite à ces SNPs sélectionnés et nous estimerions l'héritabilité avec un erreur standard plus petite que celle nous aurions obtenue en utilisant la matrice de SNPs entière.

Nous allons dans un premier temps présenter les méthodes et les résultats qui existent sur la sélection de variables dans les effets fixes des modèles linéaires mixtes parcimonieux.

Etat de l'art

Bien que la littérature concernant la sélection de variable dans les modèles linéaires mixtes soit moins riche que celle sur les modèles linéaires, plusieurs méthodes ont été proposées pour sélectionner des variables dans la partie des effets fixes et des effets aléatoires.

Pour ce qui est de la sélection dans les effets fixes, nous renvoyons le lecteur à la revue proposée par [START_REF] Müller | Model selection in linear mixed models[END_REF]. Concernant la sélection dans les effets aléatoires, nous avons connaissance du travail de Fan & Li (2012) and [START_REF] Bondell | Joint variable selection for fixed and random effects in linear mixed-effects models[END_REF]. [START_REF] Bondell | Joint variable selection for fixed and random effects in linear mixed-effects models[END_REF] ont proposé une méthode basée sur un algorithme EM pour sélectionner conjointement des effets fixes et aléatoires. [START_REF] Fan | Variable selection in mixed effects models[END_REF] ont utilisé un critère pénalisé avec une pénalité particulière, appelée SCAD (Smoothly Clipped Absolute Deviation), qui combine des pénalités 1 et 2 . Ces deux méthodes peuvent demander des temps de calcul important en grande dimension, pour l'une à cause de l'algorithme EM et pour l'autre à cause de la validation croisée nécessaire au choix des deux paramètre de régularisation.

La sélection de variable dans des cadres de très grande dimension comme ceux auxquels nous nous intéressons est un problème difficile, comme montré par Verzelen (2012) qui a étudié le cas du modèle linéaire défini par l'équation (6.7). Verzelen (2012) a en effet établi que si la condition N q log(1/q) >> n était vérifiée, ce qui est le cas notamment si le nombre de SNPs causaux (nombre de composantes non nulles dans les effets aléatoires) est grand par rapport au nombre d'individus, le support ne peut pas être entièrement retrouvé.

Concernant l'estimation de l'héritabilité, l'idée d'introduire une étape de sélection de variable au préalable a été proposée par [START_REF] Guan | Bayesian variable selection regression for genome-wide association studies and other large-scale problems[END_REF] dans le cadre Bayésien. [START_REF] Guan | Bayesian variable selection regression for genome-wide association studies and other large-scale problems[END_REF] ont en effet proposé une approche, appelée BVSR (Bayesian Variable Selection Regression), qui est très efficace dans des cas où les effets aléatoires sont très parcimonieux mais qui est biaisé quand le nombre de SNPs causaux est élevé. [START_REF] Zhou | Polygenic modeling with bayesian sparse linear mixed models[END_REF] ont ensuite proposé une approche, appelée BSLMM (Bayesian Sparse Linear Mixed Model), définie comme un estimateur hybride qui a un comportement proche de BVSR dans les cas très parcimonieux et comme l'estimateur du maximum de vraisemblance (sans sélection) sinon.

Ces résultats numériques obtenus par [START_REF] Zhou | Polygenic modeling with bayesian sparse linear mixed models[END_REF] rejoignent les résultats théoriques obtenus par [START_REF] Verzelen | Adaptive estimation of high-dimensional signal-to-noise ratios[END_REF] avec leur procédure adaptative optimale dans le cas du modèle linéaire décrit dans la section 6.2.1.

Contribution

Méthodologie Nous avons proposé une méthode de sélection de variables dans le but d'améliorer la précision de l'estimation de l'héritabilité. Ce travail a fait l'objet d'un article soumis à la publication et qui est décrit dans le chapitre 3 de ce manuscrit. Notre méthode est implémentée dans le package R EstHer, disponible sur le CRAN.

Notre approche présente deux caractéristiques principales : premièrement, elle est très efficace d'un point de vue statistique car elle permet de réduire considérablement la taille des intervalles de confiance par rapport à des méthodes sans sélection. Deuxièmement, elle est très rapide, ce qui la rend facilement utilisable sur des jeux de données très grands qui sont fréquents dans les études génétiques. Notre méthode peut gérer des cas de très grande dimension grâce à une étape de Sure Independence Screening développée par [START_REF] Ji | UPS delivers optimal phase diagram in high-dimensional variable selection[END_REF]. Nous appliquons ensuite un critère LASSO [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] combiné avec la méthode de "stability selection" [START_REF] Meinshausen | Stability selection[END_REF]. Nous proposons également une méthode de bootstrap non paramétrique pour calculer des intervalles de confiance, qui ont ensuite été validés sur des données simulées.

Au cours de notre étude numérique, nous avons observé des conclusions similaires à celles de [START_REF] Zhou | Polygenic modeling with bayesian sparse linear mixed models[END_REF] dans le cadre Bayésien : dans des cas très parcimonieux (par exemple, 200 SNPs causaux sur 100000), l'estimateur qui inclut une étape de sélection de variable est non biaisé et sa variance empiriquement est considérablement plus faible que celle de l'estimateur du maximum de vraisemblance. Cependant, quand le nombre de SNPs causaux est élevé, inclure une étape de sélection de variables conduit à sous-estimer fortement l'héritabilité. Nous avons développé un critère à appliquer sur les données pour avoir une idée de la parcimonie et pour décider s'il est judicieux ou non d'appliquer une technique de sélection de variables avant d'estimer l'héritabilité.

Nous avons donc construit un estimateur hybride capable de s'adapter au régime de parcimonie et nous avons montré sur des données simulées que cette procédure permet de réduire considérablement les intervalles de confiance dans les cas très parcimonieux par rapport à l'estimateur du maximum de vraisemblance.

L'avantage de notre méthode par rapport à celle de Zhou et al. (2013) est surtout d'avoir un temps de calcul beaucoup plus faible et aussi de ne pas avoir à régler les différents paramètres nécessaires à la convergence des algorithmes MCMC.

Applications en neuro-anatomie et en génétique animale

Nous avons appliqué notre méthode à deux jeux de données.

Le premier provient du projet européen IMAGEN, qui est une étude sur la santé mentale des adolescents. Nous avons estimé l'héritabilité du volume du cerveau et des volumes des différentes régions sous-corticales. Six des neuf phénotypes n'ont pas passé le critère de sélection que nous avons proposé, nous pouvons donc penser qu'un grand nombre de SNPs sont impliqués dans leurs variations. Nous avons obtenus avec notre méthode des résultats similaires à ceux obtenus avec des résultats obtenus avec des approches classiques de maximum de vraisemblance, comme ceux obtenus par [START_REF] Toro | Genomic architecture of human neuroanatomical diversity[END_REF] qui ont étudié les mêmes données à l'aide du logiciel GCTA développé par [START_REF] Yang | GCTA: A tool for genome-wide complex trait analysis[END_REF]. Cependant, pour les trois autres phénotypes, nous avons obtenu des estimations de l'héritabilité avec de très petites erreurs standard mais également une liste de SNPs potentiellement impliqués dans les variations et dont la pertinence pourrait être étudiée d'un point de vue biologique. Cette application est décrite dans le Chapitre 3 de ce manuscrit, après la présentation de notre méthode et la validation de celle-ci sur des données simulées.

La deuxième application est l'étude d'une espèce de truites appelée Salmo trutta. Cette truite d'eau douce peut ou non, durant sa vie, décider de quitter l'eau douce pour migrer vers la mer. Cette migration a un impact majeur sur la conservation de la truite, et nous voudrions comprendre les raisons de cette décision. Il apparaît que la croissance durant la phase en eau douce pourrait être un facteur déterminant sur la décision des truites : en effet, si un poisson a une bonne croissance en eau douce, il n'a pas d'intérêt à aller en mer où les chances de survie sont plus faibles. Néanmoins, si la truite n'arrive pas à grandir en eau douce, le bénéfice de rejoindre la mer où les chances de grandir sont plus importantes peut compenser les risques associés à ce changement. Nous souhaiterions alors rechercher la proportion de facteurs génétiques et environnementaux impliqués dans les variations de la taille des truites, et si possible déterminer quels SNPs et quelles variables environnementales sont associées à la croissance des truites. Le jeu de données contient la taille de 192 truites dont le génotype est décrit par 4069 SNPs. Cette application est décrite dans le chapitre 4.

Estimation de l'héritabilité pour des traits binaires

Nous nous intéressons à l'extension des méthodes précédentes à l'estimation de l'héritabilité d'une maladie, auquel cas les observations sont catégorielles (patient ou contrôle). Nous avons trouvé dans la littérature différents modèles utilisés pour estimer l'héritabilité de traits binaires.

6.4.1 Modèle linéaire mixte généralisé et "liability model".

Une généralisation intuitive des travaux précédents à l'estimation de l'héritabilité de traits binaires serait de considérer le modèle linéaire généralisé suivant : Y i ∼ B(q i ), (6.13) avec q i = g(l i ) où g est une fonction de lien et l i est défini par l = Zu + e, (6.14) avec u ∼ N (0, σ 2 u ) et e ∼ N (0, σ 2 e ), comme dans le LMM classique défini à la Section 6.2. Un choix classique de fonction de lien dans le cas de données binaires est par exemple g(x) = exp(x) 1 + exp(x) , ce qui garantit que q i ∈ (0, 1). L'héritabilité peut alors être définie pour la variable continue l, et donc la définition est identique à celle que l'on a considérée précédemment pour des phénotypes gaussiens :

η = N σ 2 u N σ 2 u + σ 2 e .
(6.15)

Un autre modèle, plus ancien et plus fréquemment utilisé, a été proposé par [START_REF] Falconer | The inheritance of liability to certain diseases, estimated from the incidence among relatives[END_REF], qui a fait l'hypothèse que les observations binaires pouvaient être vues comme des indicatrices qu'une certaine variable continue et non observée dépasse un seuil t : Y i = 1 {l i >t} (6.16) avec l i définie par la même expression (6.14) que dans le modèle précédent.

Ce modèle est usuellement appelé le "liability model" [START_REF] Falconer | The inheritance of liability to certain diseases, estimated from the incidence among relatives[END_REF], [START_REF] Lee | Estimating missing heritability for disease from genome-wide association studies[END_REF], [START_REF] Tenesa | The heritability of human disease: estimation, uses and abuses[END_REF]). L'héritabilité est alors également définie à l'échelle continue, avec la même définition que celle donnée dans l'équation (6.15).

Méthodes existantes

Dans la littérature spécifique à l'estimation de l'héritabilité de phénotypes binaires, nous avons trouvé des méthodes associées à chacun des modèles décrits précédemment. Pour le premier modèle (6.13), [START_REF] Villemereuil | Comparing parent-offspring regression with frequentist and Bayesian animal models to estimate heritability in wild populations: a simulation study for Gaussian and binary traits[END_REF] ont proposé d'estimer la variance σ 2 u des effets aléatoires en utilisant des méthodes MCMC développées par [START_REF] Hadfield | MCMC methods for multi-response generalized linear mixed models: The MCMCglmm R package[END_REF], puis d'estimer l'héritabilité par le ratio η = σ2 u σ2 u + 1 + 1 , où le premier 1 du dénominateur représente la variance résiduelle et le deuxième 1 la distribution spécifique à la fonction de lien probit [START_REF] Nakagawa | Repeatability for gaussian and non-gaussian data: a practical guide for biologists[END_REF]. La variance résiduelle est en effet fixée à 1 car les données binaires ne permettent pas de fournir suffisamment d'information pour inférer les deux variances σ 2 u et σ 2 e . Puisque l'expression de la vraisemblance n'est pas possible à optimiser directement, [START_REF] Breslow | Approximate inference in generalized linear mixed models[END_REF] ont proposé de maximiser une quasi-vraisemblance pénalisée, en utilisant une approximation de Laplace de la vraisemblance. Il n'existe pas de résultats théoriques concernant ces méthodes mais leurs performances numériques ont été étudiées et comparées par de [START_REF] Villemereuil | Comparing parent-offspring regression with frequentist and Bayesian animal models to estimate heritability in wild populations: a simulation study for Gaussian and binary traits[END_REF].

Concernant les procédures qui s'appuient sur le second modèle défini par les équations (6.14) et (6.16), [START_REF] Lee | Estimating missing heritability for disease from genome-wide association studies[END_REF] ont proposé d'utiliser une approche de maximum de vraisemblance comme si les traits binaires étaient gaussiens, puis d'appliquer un facteur multiplicatif pour corriger cette approximation. [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF] ont montré que cet estimateur était biaisé dans plusieurs scenarios réalistes, en particulier qu'il était très sensible à la prévalence de la maladie (quand la maladie est rare, le biais croît). Cet estimateur sous-estime également l'héritabilité quand celle-ci est élevée. [START_REF] Weissbrod | Accurate liability estimation improves power in ascertained case-control studies[END_REF] ont proposé une approche du maximum de vraisemblance pour reconstruire la variable gaussienne non observée l avant d'estimer l'héritabilité. A notre connaissance, aucune des méthodes précédemment présentées ne possède de validation théorique. De plus, elles ne prennent pas en compte un élément essentiel des données de cascontrôles qui proviennent des études médicales. En effet, bien que la maladie étudiée puisse être rare, dans une étude médicale, le nombre de patients est à peu près égal au nombre de contrôles, ce qui veut dire que la proportion de patients dans l'étude peut être très différente de la proportion de patients dans la population générale. La méthode proposée par [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF] est une méthode de moments qui prend en compte ce sur-échantillonnage des patients et c'est, à notre connaissance, la seule méthode dans ce cas. Plus précisément, [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF] Puisque E(p i p j |S = 1) n'a pas d'expression que l'on peut calculer explicitement, [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF] ont eu l'idée d'utiliser le fait que les corrélations G i,j soient petites et ont proposé une approximation grâce à des développements de Taylor autour de la quantité G i,j .

Contribution

Comme la méthode de [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF] semblait très efficace en pratique et puisque c'est la seule méthode que nous ayons vue qui prenne en compte la spécificité des données de cascontrôles, nous avons cherché à établir les propriétés théoriques de leur estimateur dans le cas où n et N tendent vers l'infini, et n/N tend vers a ∈ (0, +∞). Nous avons considéré le modèle défini par les équations (6.16) et (6.14), et nous avons supposé que Z était une matrice aléatoire dont les colonnes sont centrées et réduites comme nous l'avions défini à la Section 6.2. Dans ce modèle, la matrice de variance covariance de (l i , l j ) s'écrit :

Σ (N ) = 1 + η (G N (i, i) -1) η G N (i, j) η G N (i, j)
1 + η (G N (i, i) -1), où, pour tous 1 ≤ i, j ≤ n,

G N (i, j) = 1 N N k=1
Z i,k Z j,k . (6.17) L'idée principale est de remarquer que les quantités G N (i, j), G N (i, i) -1 et G N (j, j) -1 sont petites, ce qui implique que la matrice Σ (N ) est proche de l'identité. Inspirés par la méthode de [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF], nous avons proposé des approximations au premier et au second ordre de E[p i p j |Z, S = 1] grâce à des développement de Taylor autour de G N (i, j), G N (i, i) -1 et G N (j, j) -1, qui sont ici des variables aléatoires. Abstract : We study statistical methods to estimate the heritability of a biological trait, which is the proportion of variations of this trait that can be explained by genetic factors. First, we propose to study the heritability of quantitative traits using high-dimensional sparse linear mixed models. We investigate the theoretical properties of the maximum likelihood estimator for the heritability and we show that it is a consistent estimator and that it satisfies a central limit theorem with a closedform expression for the asymptotic variance. This result, supported by an extended numerical study, shows that the variance of our estimator is strongly affected by the ratio between the number of observations and the size of the random genetic effects. More precisely, when the number of observations is small compared to the size of the genetic effects (which is often the case in genetic studies), the variance of our estimator is very large. This motivated the development of a variable selection method in order to capture the genetic variants which are involved the most in the phenotypic variations and provide more accurate heritability estimations. We propose then a variable selection method adapted to high dimensional settings and we show that, depending on the number of genetic variants actually involved in the phenotypic variations, called causal variants, it was a good idea to include or not a variable selection step before estimating heritability. The last part of this thesis is dedicated to heritability estimation for binary data, in order to study the proportion of genetic factors involved in complex diseases. We propose to study the theoretical properties of the method developed by [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF] for casecontrol data, which is very efficient in practice. Our main result is the proof of the consistency of their heritability estimator.

  Context . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.1 Definition of heritability . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.1.2 Heritability in human genetics . . . . . . . . . . . . . . . . . . . . . . . 6 1.1.3 Heritability in vegetal and animal genetics . . . . . . . . . . . . . . . . . 7 1.2 Heritability estimations in high dimensional linear mixed models . 8 1.2.1 State of the art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.2.2 Contribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Variable selection in the random effects of a high dimensional sparse linear mixed model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.3.1 State of the art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Contribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Heritability estimation for binary traits . . . . . . . . . . . . . . . . . 14 1.4.1 Generalized Linear Mixed Model and Liability Model . . . . . . . . . . 1.4.2 Existing methods for heritability estimation in the liability model . . . . 1.4.3 Contribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2 Figure 2

 22 Figure 2.2 -Boxplots of η for different values of a, for η = 0.3 (left), η = 0.5 (middle) and η = 0.7 (right). The horizontal line corresponds to the true value of η . The whiskers of each boxplot correspond to the first and third quartiles.

Figure 2 . 6 -

 26 Figure 2.6 -Boxplots of η for different values of q, with η = 0.7 and a = 1 (top left), a = 0.5 (top right), a = 0.1 (bottom left) and a = 0.01 (bottom right). The boxplots are based on 500 replications. The whiskers of each boxplot are the fist and third quartile.

Chapter 2 -Figure 2

 22 Figure2.8 -Values of n -1/2 τ n with the real value of q (q = 0.5) ("•"), q = 1 (dotted line) and the empirical standard deviation of (η -η ) (plain line) for η = 0.7.

  Content 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.2 Description of the data . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.3 Description of the method . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.3.1 Variable selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.3.2 Heritability estimation and confidence interval . . . . . . . . . . . . . . 53 3.3.3 Additional fixed effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.4 Numerical study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.4.1 Simulation process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.4.2 Results in very sparse scenarios . . . . . . . . . . . . . . . . . . . . . . . 55 3.4.3 Results when the number of causal SNPs is high . . . . . . . . . . . . . 59 3.5 A criterion to decide whether we should apply EstHer or HiLMM . 60 3.6 Results after applying the decision criterion and comparison to other methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.6.1 Statistical performances . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.6.2 Computational times . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.7 Applications to genetic data . . . . . . . . . . . . . . . . . . . . . . . . 64 3.7.1 Calibration of the threshold . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.7.2 Application of the decision criterion . . . . . . . . . . . . . . . . . . . . 66

  Figure 3.1, which comes from Toro et al. (2015), is a schematic representation of these different areas of the brain. The data set contains n = 2087 individuals and N = 273926 SNPs, as well as a set of fixed effects, which in our case are the age (between 12 and 17), the gender and the city of residency (London, Nottingham, Dublin, Dresden, Berlin, Hamburg, Mannheim and Paris).

Figure 3 . 1 -

 31 Figure 3.1 -Different regions of the brain (this figure is taken from Toro et al. (2015)).

Chapter 3 -

 3 Figure3.2 -Absolute difference between η and η for thresholds from 0.6 to 0.9 and for q = 10 -3 (100 causal SNPs).

  Figure 3.4 -Estimated value of the heritability with 95 % confidence intervals. The results are displayed for several values of η : 0.5, 0.6 and 0.7. (a) The data sets were generated including fixed effects. (b) The matrix Z used to generate data sets comes from the IMAGEN data. The black dots correspond to the mean of η over 10 replications and the crosses are the real value of η .

  Figure 3.5 -(a) Boxplots of the length of the set of selected variables with EstHer for 40 repetitions. The real number of non null components is 100. (b) Boxplots of the proportion of the real non null components captured in the set of selected variables.

Chapter 3 -Figure 3 . 7 -

 337 Figure3.7 -Results of HiLMM and EstHer when there are a few causal SNPs with moderate effects and a lot of SNPs with small effects. The proportion of each is 100 out of 1000 (up) and 100 out of 10000 (bottom), with η = 0.4 and 0.6.

Figure 3

 3 Figure 3.11 -Estimations of η with 95 % confidence intervals obtained using EstHer, BSLMM, HiLMM and GEMMA-LMM with 100 causal SNPs (top) and 10,000 causal SNPs (bottom). The results are obtained with 10 replications.

Figure 3

 3 Figure3.12 -Times (in seconds) to compute one heritability estimation with BSLMM (crosses) and EstHer (dots) by using 16 thresholds for n = 2000 and different values of N from 50, 000 to 200, 000.

  Figure 3.13 -Absolute value of the difference between η and η for thresholds from 0.6 to 0.9, and for different values of qN : (a) 50 causal SNPs, (b) 100 causal SNPs. Each difference has been computed as the mean of 10 replications.

Content 4 . 1

 41 Data and motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.2 Model and method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.2.1 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.2.2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.3 Simulation Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.3.1 Results obtained without selection . . . . . . . . . . . . . . . . . . . . . 71 4.3.2 Calibration of the threshold in the stability selection . . . . . . . . . . . 71 4.3.3 Results regarding heritability estimation . . . . . . . . . . . . . . . . . . 71 4.3.4 Support recovery when the selection is applied . . . . . . . . . . . . . . 72 4.4 Results on the juvenile size dataset . . . . . . . . . . . . . . . . . . . 74 4.4.1 Results obtained without taking fixed effects into account . . . . . . . . 74 4.4.2 Including fixed effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

Figure 4

 4 Figure 4.1 -Heritability estimations obtained from 100 replications without selection for η = 0.6 and different values of q from 0.001 to 0.1.

Chapter 4 -Figure 4 . 2 -

 442 Figure 4.2 -Absolute difference |η -η| for thresholds from 0.2 to 0.8 and for different numbers of causal SNPs, for different values of η : 0.4 (green), 0.5 (dark blue), 0.6 (blue) and 0.7 (purple).

Figure 4 . 3 -

 43 Figure 4.3 -Estimated values of η obtained with our variable selection method EstHer with a threshold of 0.35 and without selection (HiLMM) for different numbers of causal SNPs.

Figure 4

 4 Figure 4.4 -Proportion of the support recovered by our selection method Esther (top) and number of selected variables (bottom) for different numbers of causal SNPs and for different thresholds.

Figure 4 . 6 -

 46 Figure 4.6 -Estimated value of the heritability using different fixed effects, from left to right: river, river/month, reporting groups, reporting groups and month.

Figure 4

 4 Figure 4.8 -Estimated value of the heritability with and without alkaninity and hardness as covariates, but removing 12 individuals from the study.

Content 5 . 1

 51 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Model and definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.1 Liability model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.2 Case control study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Heritability estimator . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3.1 Method of Golan et al. (2014) . . . . . . . . . . . . . . . . . . . . . . . . 5.3.2 Our method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Consistency of the heritability estimator η(1) . . . . . . . . . . . . . . 5.5 Second order approximation of E[W i W j |Z, i = j = 1] . . . . . . . . . 5.6 Numerical study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6.1 Simulation process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8 Proofs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8.1 Taylor development of E[W i W j |Z, i = j = 1] in Model (5.7) . . . . . . 5.8.2 Proof of Theorem 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8.3 Second order approximation of E[W i W j |Z, i = j = 1] . . . . . . . . .

5. 3 . 1

 31 Method of Golan et al. (2014) Golan et al. (2014) considered a simplified version of Model (5.4), where the liability is given by

Figure 5

 5 Figure 5.2 -Time in seconds to compute an estimation of η obtained with η(1) (dots) and η(2) (triangles) for n = 100 (left) and n = 1000 (right) and for different values of N (from 1000 to 10 5 ).

  .10) De plus, comme dans notre cas nous nous intéressons uniquement à l'estimation de η , nous avons remplacé σ 2 par son

Titre:

  Estimation de l'héritabilité dans les modèles mixtes en grande dimension : théorie et applications Mots clés : Grande dimension, Héritabilité, Modèles mixtes, Sélection de variables. Résumé : Nous nous intéressons à des méthodes statistiques pour estimer l'héritabilité d'un caractère biologique, qui correspond à la part des variations de ce caractère qui peut être attribuée à des facteurs génétiques. Nous proposons dans un premier temps d'étudier l'héritabilité de traits biologiques continus à l'aide de modèles linéaires mixtes parcimonieux en grande dimension. Nous avons recherché les propriétés théoriques de l'estimateur du maximum de vraisemblance de l'héritabilité : nous avons montré que cet estimateur était consistant et vérifiait un théorème central limite avec une variance asymptotique que nous avons calculée explicitement. Ce résultat, appuyé par des simulations numériques sur des échantillons finis, nous a permis de constater que la variance de notre estimateur était très fortement influencée par le ratio entre le nombre d'observations et la taille des effets génétiques. Plus précisément, quand le nombre d'observations est faible comparé à la taille des effets génétiques (ce qui est très souvent le cas dans les études génétiques), la variance de l'estimateur était très grande. Ce constat a motivé le développement d'une méthode de sélection de variables afin de ne garder que les variants génétiques les plus impliqués dans les variations phénotypiques et d'améliorer la précision des estimations de l'héritabilité. La dernière partie de cette thèse est consacrée à l'estimation d'héritabilité de données binaires, dans le but d'étudier la part de facteurs génétiques impliqués dans des maladies complexes. Nous proposons d'étudier les propriétés théoriques de la méthode développée par Golan et al. (2014) pour des données de cas-contrôles et très efficace en pratique. Nous montrons notamment la consistance de l'estimateur de l'héritabilité proposé par Golan et al. (2014). Title: Heritability estimation in high-dimensional mixed models : theory and applications Keywords : High Dimension, Heritability, Mixed Models, Variable Selection.

  Existing methods for heritability estimation . . . . . . . . . . . . . . 23 2.4 Theoretical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.5 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.7 Proofs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.7.1 Proof of Theorem 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.7.2 Proof of Theorem 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.7.3 Proof of Theorem 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.7.4 Proofs of technical lemmas . . . . . . . . . . . . . . . . . . . . . . . . . 41

Bonnet, E. Gassiat, C. Lévy-Leduc, "Heritability estimation in high dimensional sparse linear mixed models", Electronic Journal of Statistics, 9(2):2099-2129, 2015. The method which is presented is implemented in the HiLMM R package, available on the CRAN. Content 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.2 Model and heritability estimator . . . . . . . . . . . . . . . . . . . . . 21 2.2.1 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.2.2 Heritability estimator . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.3 2.5.1 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.5.2 Results in Model (2.2) when q = 1 . . . . . . . . . . . . . . . . . . . . . 27 2.5.3 Results in Model (2.2) when q < 1 . . . . . . . . . . . . . . . . . . . . . 29 2.6

  3.7.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

Table 3 .

 3 1 -95 % confidence intervals for η obtained empirically and by our Bootstrap method.

	η	0.4	0.5	0.6	0.7
	Bootstrap [0.353 ; 0.503] [0.413 ; 0.565] [0.494 ; 0.654] [0.596 ; 0.738]
	Empirical [0.391 ; 0.470] [0.449 ; 0.542] [0.496 ; 0.645] [0.618 ; 0.720]

Table 3 .

 3 4 -Mean value of the number of overlapping confidence intervals for 16 thresholds from 0.7 to 0.85.

	Phenotype Number of thresholds
	Bv	7.19
	Hip	7.5
	Icv	7.37
	Acc	9.94
	Amy	9.88
	Th	7.5
	Ca	7.13
	Pu	7.13
	Pa	10.75

  Chapter 3 -Improving heritability estimation by a variable selection approach in high dimensional sparse linear mixed models

	1.0
	0.8
	0.6
	0.4
	0.2
	0.0
	bv	icv	th	pu	pa	hip	amy	acc	ca

Table 5

 5

	n	N	1000	10000	50000	10 5
	100 η(1)	0.478	2.390	28.595	40.528
		η(2)	3.148	7.127	56.761	86.156
	1000 η(1) 69.047 327.240 2887.518 7845.281
		η(2) 376.363 936.845 6624.186 13500.510

.1 -Times in seconds to compute an estimation of η obtained with η(1) and η(2) for different values of n (100 and 1000) and N (from 1000 to 10 5 ).

  Tous les caractères biologiques, par exemple la taille ou le poids, sont influencés à la fois par des facteurs génétiques et environnementaux. Quantifier ces deux contributions pour un trait donné est une problématique difficile et fondamentale en biologie. Le concept d'héritabilité désigne la part de la variabilité d'un trait observé (ou phénotype) qui peut être attribuée à des causes génétiques.

6.1 Contexte biologique

  [START_REF] Pirinen | Efficient computation with a linear mixed model on large-scale data sets with applications to genetic studies[END_REF] ont remarqué que le modèle défini par l'équation (6.1) peut être reparamétrisé en fonction de β, η et σ 2

				Chapitre 6 -Résumé en français
				.5)
	Searle et al. (1992) ont regroupé de nombreuses techniques d'optimisation pour estimer les
	paramètres β, σ 2 u et σ 2 e , parmi lesquelles nous pouvons citer les équations d'Henderson ou bien
	des méthodes itératives comme Fisher-Scoring et Newton-Raphson.
	Une idée naturelle pour estimer l'héritabilité est d'estimer les paramètres de variances σ 2 u et σ 2 e afin d'obtenir un estimateur de η en tant que ratio :
	N σu	2 /(N σu	2 + σe	2 ).

  ont considéré une version simplifiée du modèle (6.14), où la variable continue l est définie par l = g + e, où g est un effet génétique aléatoire, dont les composantes sont corrélées, et e est un effet aléatoire environnemental, que l'on suppose indépendant des effets génétiques. Les deux effets sont supposés gaussiens : e a une variance égale à (1 -η )Id R n et g a une matrice de variance covariance dont les entrées diagonales sont égales à η et le terme non diagonal (i, j) est égal à 137 Chapitre 6 -Résumé en français η G i,j . Pour 1 ≤ i = j ≤ n, la matrice de variance covariance de (l i , l j ) est donnée par où P est la proportion de cas dans l'étude et l'événement {S = 1} est réalisé si les individus i et j sont sélectionnés dans l'étude. L'estimateur de l'héritabilité proposé par Golan et al. (2014) est un estimateur des moindres carrés obtenu en minimisant le critère i =j (p i p j -E(p i p j |S = 1)) 2 .

	1 G i,j η	G i,j η 1	.
	Ils ont ensuite défini la variable
	p i =	Y i -P P (1 -P )	,

i = V i Y i + U i (1 -Y i ),

Nous présentons à présent l'intérêt d'estimer l'héritabilité de caractères humains. C'est en effet un premier pas vers la compréhension de maladies complexes, qui ont souvent des causes multiples. Nous faisons référence en particulier à des maladies qui ne sont pas causées par un seul gène affecté mais dont on soupçonne néanmoins qu'elles ont une composante génétique importante, probablement répartie sur plusieurs gènes. Par exemple, les causes de certaines maladies psychiatriques comme l'autisme et la schizophrénie sont encore vagues aujourd'hui. Une composante génétique est suggérée par les résultats des études de jumeaux monozygotes et dizygotes (les jumeaux monozygotes ont des génomes identiques alors que les jumeaux dizygotes ont environ 50% de leurs génomes en commun). Ces études montrent que si l'un des jumeaux souffre de troubles autistiques, l'autre également dans 82 à 92 % des cas pour les jumeaux monozygotes[START_REF] Bailey | Autism as a strongly genetic disorder: evidence from a british twin study[END_REF]) contre 20 % des cas chez les jumeaux dizygotes[START_REF] Hallmayer | Genetic heritability and shared environmental factors among twin pairs with autism[END_REF]. Ces études décrivent l'autisme comme la maladie psychiatrique avec la plus importante composante génétique. Cependant, la sévérité des traits autistiques (troubles des interactions et du langage, handicap mental...) peuvent varier fortement chez deux patients qui présentent des causes similaires, par exemple la même mutation. Ces observations laissent penser, comme c'est le cas pour d'autres maladies génétiques, que le fond génétique module les effets d'une mutation et rend certains individus plus ou moins sensibles au développement de traits autistiques.De plus, toutes ces études montrent que, malgré des patrimoines génétiques identiques, la concordance des symptômes chez les jumeaux monozygotes n'est jamais totale, ce qui confirme une composante épigénétique et/ou environnementale. Néanmoins, quantifier ces possibles causes et leurs interactions potentielles reste une question difficile.

Université Paris-Saclay Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France

Chapter 5 -Heritability estimation in case-control studies intervals for η .

All these methods obviously cannot be applied directly for binary traits, but there exist extensions by assuming an underlying Gaussian variable linked to the binary phenotype. There are two different modelings which connect binary phenotypes to a continuous quantity called the liability. The first one consists in assuming that the observations Y 1 , . . . , Y n are distributed according to the following Generalized Linear Mixed Model (GLMM):

(5.3) with p i = g(l i ) where g is a link function and l i is defined as l = Zu + e (5.4) with l = (l 1 , . . . , l n ), u ∼ N (0, σ 2 u ) and e ∼ N (0, σ 2 e ), like in classical LMM defined in Equation (5.1). The heritability is then defined "at the liability scale", which means for the continuous variable l, and is given by the same expression (5.2) as for quantitative traits. Several methods were established to estimate heritability in Model (5.3): among them we can quote the MCMC method of [START_REF] Hadfield | MCMC methods for multi-response generalized linear mixed models: The MCMCglmm R package[END_REF] and the penalized quasi-likelihood approach of [START_REF] Breslow | Approximate inference in generalized linear mixed models[END_REF]. The numerical performance of these methods can be found in the comparative study of [START_REF] Villemereuil | Comparing parent-offspring regression with frequentist and Bayesian animal models to estimate heritability in wild populations: a simulation study for Gaussian and binary traits[END_REF].

Another modeling and definition for the heritability of a binary trait, which is more frequently used than the previous one, was proposed by [START_REF] Falconer | The inheritance of liability to certain diseases, estimated from the incidence among relatives[END_REF], who assumed that the binary observations could be seen as an indicator function of a Gaussian variable exceeding a given threshold t:

(5.5) with l i defined by the same expression (5.4) than in Model (5.3). Observe that the threshold t is directly linked to the prevalence of the disease in the population, that is the proportion K of the population which is affected by the disease. Indeed, K = P(Y i = 1) = P(l i > t).

(5.6)

The unobserved Gaussian variable l = (l 1 , . . . , l n ) is also called the liability in this modeling, which is usually called the "liability model" [START_REF] Falconer | The inheritance of liability to certain diseases, estimated from the incidence among relatives[END_REF], [START_REF] Lee | Estimating missing heritability for disease from genome-wide association studies[END_REF], [START_REF] Tenesa | The heritability of human disease: estimation, uses and abuses[END_REF]) and has been shown to be a reasonable modeling for complex diseases, for instance by [START_REF] Purcell | Common polygenic variation contributes to risk of schizophrenia and bipolar disorder[END_REF]. The heritability is then also defined as the heritability at the liability scale as in Equation (5.2).

Regarding the procedures based on the second modeling defined in Equations (5.5) and (5.4), [START_REF] Lee | Estimating missing heritability for disease from genome-wide association studies[END_REF] proposed to use a maximum likelihood approach as if the binary traits were Gaussian, and then to apply a multiplicative factor to correct this approximation. [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF] showed that this heritability estimator was strongly biased in several realistic scenarios, in particular it was very sensitive to the prevalence of the disease, which is the proportion of the Chapter 5 -Heritability estimation in case-control studies

The proof of Theorem 4 relies on the following lemmas.

Lemma 7. When n and N go to infinity and n/N goes to a,

G N (i, j) 2 converges in probability to a.

We will then have to focus on

(5.21)

Let E N be the following event

N the complement of the event E N . We consider the following decomposition η(1) = η(1) 1 E N + η(1) 1 E c N . Lemma 8. For all values of q, the probability of

Using the result of Lemma 8, η(1) 1 E c N converges in probability to 0 since

Lemma 9. When n and N go to infinity and n/N goes to a ∈ (0, +∞),

converges in probability to acη , where c is defined in Equation (5.17).

Lemma 10. When n and N go to infinity and n/N goes to a ∈ (0, +∞),

The results of Lemmas 9 and 10 achieve the proof of Theorem 4. The proof of Lemmas 7, 8, 9 and 10 are given in Section 5.8.2.

Chapter 5 -Heritability estimation in case-control studies η = 0.5 η = 0.7 q q q q q q q q q q q q q q q q q q q first approx second approx 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q q q q q q q q q q first approx second approx 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q q q q q q q q q q q q q first approx second approx 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q q q q q q q first approx second approx when the ratio n/N goes to a constant a. It would be interesting to complete this work with theoretical results which could allow the user to compute accurate confidence intervals, similarly to existing results for quantitive traits. As it is often the case in genetic applications, the question of removing strong assumptions such as the Gaussianity of the random effects or the independence of the columns of the SNP matrix

Proof of Lemma 8

Note that

Let δ be a positive real number such that √ δ/2c ∈ V 0 and δ ≤ δ min 4 , where V 0 and δ min are defined in Assumptions 2 and 2.1 respectively.

Observe that n i=1

(5.30)

Let us show that

(5.31)

By Chernoff inequality, for all λ ≥ 0,

Then, by Assumption 1.2, for all positive values of λ in V 0 ,

(5.32)

Chapter 5 -Heritability estimation in case-control studies

The right term of (5.32) is maximum when

, which implies that

Similarly, for all negative values of λ in V 0 ,

(5.33)

The right term of (5.33) is maximum when

, which implies that

which proves (5.31).

P(

Since 4δ -δ min < 0 by assumption on δ, we apply again Chernoff inequality, which gives us that:

This result, combined with (5.31), proves that

Chapter 5 -Heritability estimation in case-control studies

Using Chernoff inequality and Assumption 1.1, we can prove that

and

Moreover,

Using Chernoff inequality and Assumption 1.3, we obtain that

and with Assumption 1.1 we have

with n 2 N 2 N = a 2 N 2+2γ and nN N = aN 3 2 +γ where γ > 0, which implies that the main term in the exponential is -

Similarly, we can show that

This concludes the proof that for all values of q,

Chapter 5 -Heritability estimation in case-control studies

We use similar techniques to otain an upper bound for

Since we have already proved (5.34) and (5.35), we will conclude the proof by showing that

and

(5.37)

(5.36) is obtained using Assumption 1.3 and Chernoff inequality.

which proves (5.37) and achieves the proof of Lemma 8.

Proof of Lemma 9

According to the results of Section 5.8.2, we have

Thus, we just need to prove that

We shall see that R N (i, j)1 E N may be upper bounded by a finite sum of terms of the form

Chapter 5 -Heritability estimation in case-control studies (5.38) with k in 2, 22 and

)1 E N is upper bounded by a finite sum of terms of the form

This achieves the proof of Lemma 9.

Let us explain why Equation (5.38) holds. We need to evaluate |R N (i, j)1 E N |. Then, let us look at the previous remainders which compose R N (i, j), and we will provide upper bounds when E N holds.

The remainders γ N , γN and γN are only products of α N , A N (i), A N (j) and B N (i, j).

The next remainder is µ N , which is defined as

Chapter 5 -Heritability estimation in case-control studies

Chapter 5 -Heritability estimation in case-control studies

The term of order exactly 1/N in

We will demonstrate the propositions:

Proposition 5. For all terms T N (i

Propositions 3, 4 and 5 prove (5.43). Let us prove now Propositions 2, 3, 5 and 4.

Since each case has a probability 1 and each control a probability K(1 -P )/P (1 -K) to be in the study (these probabilities are given in Equation (5.12) and (5.13)),

The computations leading to (5.44) are very similar to those leading to (5.15), which are detailed in Appendix 5.B. The development of order 0 of

Similarly,

Replacing all these expressions in (5.44) gives us that the approximation of order 0 is null, which achieves the proof of Proposition 2. Let us prove now Proposition 3. If card({i 1 , i 2 , i 3 , i 4 })=4, let us compute the approximation of order 1/

• -(1-P ) P if one individual is a control and the three others are cases. • 1 if two individuals are controls and two are cases.

• -P

1-P if one individual is a case and the three others are controls.

(1-P ) 2 if all individuals are controls.

The covariance matrix of (l i 1 , l i 2 , l i 3 , l i 4 ) is

Chapter 5 -Heritability estimation in case-control studies

Regrouping all the first terms in the expression of

Similarly we regroup the terms in η

The term in C 1,2 C 3,4 of P("3 cases, 1 control"|Z) is

The term in C 1,2 C 3,4 of P("1 case, 3 controls"|Z) is

It remains to compute the approximation of the denominator of

Chapter 5 -Heritability estimation in case-control studies

This proves Proposition 4.

Second order approximation of E[W

The density function f can still be written as

but with the explicit term of order 1/N in the expressions of

.

Chapter 5 -Heritability estimation in case-control studies Thus,

with the last term obtained by developing the exponential function. Since

Multiplying by

x 2 y 2 dxdy = (-tφ(t) + 1 -K) 2 .

Then we have

Multiplying by

Chapter 5 -Heritability estimation in case-control studies

.

Finally, we compute similarly

We replace the expressions of P(Y i = Y j = 1|Z), P(Y i = Y j = 0|Z) and P(Y i = Y j |Z) in the expression of E(W i W j |Z, i = j = 1). Since we already computed the terms of order 1 √ N for the numerator, it only remains the terms of order 1 N . Eventually, we find that the numerator can be writen as :

Similarly, we compute the expression of the denominator (at order 1 √ N since the main term of the numerator is of order 1 √ N ). We obtain the following expression:

Chapter 5 -Heritability estimation in case-control studies

By definition, the probabilities p case and p control are linked to the variables i as follows:

and

The ratio of the two following equations:

and

, with the full ascertainment assumption given by (5.12) prove equation (5.13).

5.B Proof of Equation (5.15)

This equation was proved in [START_REF] Golan | Measuring missing heritability: Inferring the contribution of common variants[END_REF], we recall the proof here for the sake of completeness. Conditionally to the event { i = j = 1}, the variable W i W j can take the following values:

Let us write the expectaction of W i W j conditionally to Z and conditionally to { i = j = 1}:

(5.56)

Chapter 5 -Heritability estimation in case-control studies

under the full ascertainment assumption given by Equation (5.12).

Similarly, since we have seen in Equation (5.13) that a control has a probability K(1-P ) P (1-K) to be selected in the study and since i and j are assumed to be independent conditionally to Z, Y i and Y j : If we combine all these computations and we plug them in the expression (5.56), we obtain (5.15).

5.C Proof of Equation (5.19)

Notice first that

Moreover, since the variables (Z i,k ) 1≤i≤n are normalized according to Equation (5.9),

Chapter 5 -Heritability estimation in case-control studies

By taking the expectation and since the variables (Z i,k ) 1≤i≤n are exchangeable, we obtain that E[Z 2 i,k ] = 1.

(5.57)

Using (2) of Proposition 1 and Equation (5.57), we obtain that

is bounded and

Similarly, Z i,k Z j,k = O p (1) for all i = j.

5.D Proof of Proposition 1

Observe that for all k = 1, . . . , N , n i=1 Z i,k = 0 (5.58)