Narrowing the gap on heritability of common disease by direct estimation in case-control GWAS

One of the major developments in recent years in the search for missing heritability of human phenotypes is the adoption of linear mixed-effects models (LMMs) to estimate heritability due to genetic variants which are not significantly associated with the phenotype. A variant of the LMM approach has been adapted to case-control studies and applied to many major diseases by Lee et al. (2011), successfully accounting for a considerable portion of the missing heritability. For example, for Crohn's disease their estimated heritability was 22% compared to 50-60% from family studies. In this letter we propose to estimate heritability of disease directly by regression of phenotype similarities on genotype correlations, corrected to account for ascertainment. We refer to this method as genetic correlation regression (GCR). Using GCR we estimate the heritability of Crohn's disease at 34% using the same data. We demonstrate through extensive simulation that our method yields unbiased heritability estimates, which are consistently higher than LMM estimates. Moreover, we develop a heuristic correction to LMM estimates, which can be applied to published LMM results. Applying our heuristic correction increases the estimated heritability of multiple sclerosis from 30% to 52.6%.


Main text
The mystery of the "missing heritability" is a term commonly used to denote the gap between the expected heritability of many common diseases, as estimated by family and twin studies, and the overall additive (narrow-sense) heritability obtained by accumulating the effects of all singlenucleotide polymorphisms (SNPs) that have been found to be significantly associated with these conditions in genome-wide association studies (GWASs) [6][7][8][9] . Many diseases which comprise a considerable portion of the health-care burden display such a gap, including type-1 and type-2 diabetes, bipolar disorder, schizophrenia, Alzheimer's disease, multiple sclerosis and Parkinson's disease.
Researchers have proposed several hypothetical solutions to this mystery. These theories include rare causative variants, which are undetected by the current GWAS methodology, common variants with small effects, which do not pass the significance threshold and are therefore unaccounted for, gene-gene and gene-environment interactions which are overlooked by the additive model assumed by the GWAS scheme, epigenetic effects and more 7,9,10 .
Clearly, different theories have major implications for our understanding of human disease, and also dictate different strategies for discovery of the underlying genetic causes of disease. For example, identifying rare variants requires a focus on deep sequencing 11 , while detecting small effects of common variants requires increasing the sample sizes dramatically, or conducting massive meta-analyses 12 . Effective planning of genetic research can therefore only be guided by a satisfactory and well founded allocation of missing heritability to its various possible sources and causes. Additionally, accurate estimates of heritability can facilitate better personalized genetic risk predictions 13 .
Yang et al. 1 pioneered the use of LMMs to estimate heritability of continuous traits (e.g. height) from GWAS data, while accounting both for significantly and insignificantly associated SNPs, thus providing an estimate of the total heritability explained by common SNPs. Their method was applied to numerous phenotypes including human height, body mass index, von Willebrand factor 14 , gene expression 15 and intelligence 16 . The LMM method was later adapted to dichotomous disease phenotypes by Lee et al. 2 , assuming the well-known liability threshold model 17 .
They showed that for the case of a random (unascertained) population sample the LMM method can be used to estimate heritability by applying it directly to the observed phenotype as if it were continuous, and correcting the resulting "observed scale" estimate as suggested by Dempster and Lerner 17 .
The ascertained case-control scenario, which is relevant for most GWAS studies, proved more challenging, as the enrichment of cases violates the assumption of normality of the genetic component critical for LMM estimation. Lee et al. proposed a complex mathematical solution for this case, and demonstrated that it worked in limited simulations. They applied their method to three phenotypes from the Wellcome Trust Case-Control Consortium (WTCCC) 18 : Crohn's disease (CD), bi-polar disorder (BD) and type-1 diabetes (T1D), and it has since been applied to many major diseases including schizophrenia 4 , Alzheimer's disease, multiple sclerosis 5 and Parkinson's disease 3 . In all these cases, accounting for insignificant SNPs resulted in much higher estimates of heritability than the estimates obtained by accumulating SNPs which were found significant in GWAS.
The basic idea of heritability estimation methods in GWAS is that individuals who are more correlated genetically are more likely to have similar phenotypes. The strength of this connection depends on the heritability -higher heritability implies stronger connection. The genetic correlation of every pair of individuals can be estimated from their genotypes. Since the individuals in the study are unrelated, the correlations between their genotypes are typically small, but non-zero. By accumulating the information across all pairs of individuals in the study, one can leverage these minor differences to separate the phenotypic variance into its genetic and environmental components, resulting in an estimate of heritability.
Ascertained case-control studies pose a considerably harder challenge to deal with theoretically than non-ascertained (prospective or observational) studies, as the fact that cases are overrepresented in the study creates a wide range of artifacts. The normality of the genetic component, assumed by the liability threshold model, is violated by over-sampling of cases, as is the assumption that the genetic and environmental components are independent. Methods that critically rely on these assumptions, like LMM, can not be applied anymore, without requiring complex and potentially questionable adjustments.
We therefore developed our GCR method for estimating heritability of polygenic phenotypes in ascertained case-control studies. The key novelty of our method is that we model the selection process directly and account for ascertainment by conditioning on the selection of observed in-dividuals. This is contrary to the existing methods of Lee et al. 2 and Zhou et al. 19 , which apply methodologies suited for prospective studies of a continuous normal phenotype, and then attempt to correct the estimates so that they account for ascertainment.
Technically, we adopt the commonly used liability threshold model 2,17 , and derive analytically the relationship between the genetic correlation of any two individuals, and the similarity between their phenotypes, conditional on the fact that both individuals were selected for the study.
While this relationship has no closed-form expression, taking the first-order approximation yields a linear relationship between the expected phenotypic similarity and the genetic correlation which depends on the heritability. Our approach then entails performing a regression of phenotype similarities between pairs of individuals on their genetic correlation, as estimated from GWAS data.
The slope of the regression is then transformed to an estimate of the heritability (Online Methods).
We applied GCR to the same three phenotypes from WTCCC: CD, BD and T1D, and four additional phenotypes: type-2 diabetes (T2D), coronary artery disease (CAD), rheumatoid arthritis (RA) and hypertension (HT). Our method yields considerably higher estimates than the LMM method 2 for CD and BD (Table 1).
To explore the nature of this considerable gap between our estimates and previously published estimates of the additive heritability, we conducted a wide range of realistic simulations.
We simulated an end-to-end generative model, starting from minor allele frequencies (MAFs) and individual SNP effects, through genotypes and phenotypes, and finally the entire selection process (Online Methods).
Our simulations differ from all previously reported simulations of case-control studies in the context of heritability estimation in that they yield realistic, rather than degenerate, correlation structure among the individuals in the study. Since we simulated genotypes, the correlation between individuals is not restricted to a small subset of possible values as in the simulations of Lee et al. 2 . Additionally, when both heritability and ascertainment are high, cases tend to be more genetically similar than expected by chance, which was not accounted for in Lee et al.'s simulations.
In particular, their assumed correlation structure is highly degenerate, with as many as 99.98% of correlations being exactly 0. For further discussion of the different simulation setups see Online Methods and Supp. Material.
Our simulations showed that using the LMM method for heritability estimation in casecontrol studies yielded considerably negatively-biased estimates, while our method consistently generated unbiased estimates ( Fig. 1 and Supp. Figs. 1-7).
The simulations also demonstrated a strong relation between the bias of the LMM approach and the increased variance of genetic effects due to ascertainment. This allowed us to derive a heuristic correction for the LMM heritability estimates, which corrected its bias successfully in all our simulation scenarios (Online Methods, Fig. 2 and Supp. Figs. [12][13][14][15]. To validate our heuristic correction, we applied it to LMM heritability estimates for all seven WTCCC phenotypes and compared the corrected estimates to the ones obtained using GCR. The correlation between the corrected LMM estimates and GCR estimates was 0.98, indicating that our heuristic correction generalizes beyond our particular simulations (See Online methods and Supp. material for more details).
We then applied our correction to the results of several published studies which used the LMM approach. As expected, for studies with low ascertainment or phenotypes with low heritability, the corrected estimates were not substantially different from published estimates (Alzheimer's disease 5 , endometriosis 5 , schizophrenia 4 and Parkinson's disease 3 ). However, the corrected heritability estimate of multiple sclerosis -the most ascertained study we inspected -is 52.6%, compared to the uncorrected estimate of 30%. See Supp. material for more details.
An important aspect of heritability estimation is the inclusion of known (fixed) effects that should be accounted for, like known associated SNPs, known environmental effects or sex. However, estimation of fixed effects from ascertained data under the normality assumption might produce biased estimators of both the fixed effects and the heritability [20][21][22] . Conversely, our method can be rigorously extended to allow for fixed effects (Online Methods). Our simulations suggest that GCR produces accurate estimates of heritability in the presence of fixed effects (see Supp. material for more details). Specifically for the WTCCC phenotypes, previous analysis 18 suggested that there's little population structure (after removing individuals of non British descent). We validated this conclusion using the statistical test described in Patterson et al. (2006) 23 (see Supp. material for more details). We therefore included only sex as a fixed effect in our heritability estimation in Table 1. The estimates of both LMM and GCR for most phenotypes without inclusion of a sex fixed effect are not substantially different from the estimates in Table 1, with the exception of the heart-disease related phenotypes (results without fixed effect not shown). For CAD in particular, it is well known that sex is a major risk factor. Accordingly, the GCR estimate without fixed effect is 71.5%, dropping to about 61% once it is added, as seen in Table 1. GCR is computationally very efficient compared to the LMM approach (running time scales quadratically rather than cubically in the number of individuals), and so the running time is seconds rather than hours on WTCCC cohorts. As the typical size of GWAS continues to increase, this efficiency can prove critical in allowing GCR to remain computationally practical compared to other approaches. Fast running time is also useful in allowing us to use resampling approaches like the jackknife 24 for estimating confidence intervals and standard errors, rather than relying on complex and potentially inaccurate parametric approximations. All standard errors and confidence intervals we report are based on the jackknife.
We also experimented with a more involved version of our approach that included the second order term of the Taylor series expansion. This did not yield more accurate results than our first order GCR in our simulations (Supp. material).
Recently, Zuk et al. 10 have shown an intriguing general result on the connection between the derivative of the dependence of phenotypic similarity on the proportion of identity by descent (IBD) at the average population IBD, and the narrow-sense heritability, regardless of population structure and phenotype-genotype architecture. While Zuk et al.'s result uses IBD and our method relies on estimated genetic correlations, the latter is, in fact, an unbiased estimator of the former, which is unknown. Hence, our method can produce unbiased estimators of narrow-sense heritability, even in the presence of population structure, by multiplying the estimate by a factor of (1-average population IBD).
In conclusion, our new proposed GCR method for estimation of heritability in case-control GWAS, which is based on a regression of phenotype similarities on genotype correlations, is shown to be efficient and accurate. It improves on existing methodology and generates substantially higher estimates of heritability for two major diseases inspected: Crohn's disease and bipolar disorder. Moreover, we provide a heuristic correction for published LMM heritability estimates, which suggests that the heritability of multiple sclerosis is also considerably larger than previously thought.

Methods
Note: An implementation of our code for GCR estimation, simulations, and heuristic correction of LMM can be accessed at: https://sites.google.com/site/davidgolanshomepage/gcr Heritability estimation using genetic correlation regression Liability threshold model -notations. Denote K the prevalence of a condition in the population and P the prevalence in the study.
Under the liability threshold model, we assume that each individual i has an unknown liability l i = g i + e i where g i is a genetic random effect, which can be correlated across individuals, and e i is the environmental random effect, which is assumed to be independent of each other and of the genetic effects. Both effects are assumed to follow a Gaussian distribution with variances σ 2 g and 1 − σ 2 g respectively. A person is then assumed to be a case if her liability exceeds a threshold t = Φ −1 (1 − K), i.e. the phenotype y i is given by y i = I{l i > t}. This definition guarantees that the prevalence in the population is indeed K.
Selection probabilities. When the study is observational, the probability of being included in the study is independent of the phenotype. However, in a case-control study, the proportion of cases is usually greatly ascertained. To model this fact, we define a random indicator variable s i indicating whether individual i was selected to the study.
The commonly used "full" or "complete" ascertainment assumption 25 is P(s i = 1 | y i = 1) = 1. While this assumption can be relaxed, as discussed later, it simplifies subsequent analysis.
Suppose the population is of size n and that the expected size of the study is n s < n. The expected number of cases in the study is nK. Additionally, the proportion of cases in the study is P , so: yielding: denote p the probability of a control being included in the study (i.e. p = P(s i = 1 | y i = 0)). The expected number of controls in the study is n(1 − K)p . Additionally, the proportion of controls in the study is (1 − P ) so: Solving for p yields: From here it follows that the probability of being included in the study for a given individual (with unknown phenotype) is: Our results do not depend strictly on the full ascertainment assumption (that is, P(s i = 1 | y i = 1) = 1). The latter assumption can be relaxed such that for any p * , as long as the probability of being selected as a control is multiplied by the same probability. This can model any step prior to the selection procedure, for example the probability that an individual is approached by the health administration to begin with. For example, a nonascertained study (K = P ) might involve only a proportion p * of the population. In this case deriving the selection probability of a given individual yields p * as expected.
Heritability estimation. Next, consider a pair of individuals in the study, whose genetic effects are correlated and denote by ρ the correlation between their genetic effects.
Denote by Z ij the product of the standardized phenotypes: The variable Z ij can take three values: We write down the expected value of Z ij , conditional on the fact that s i = s j = 1 (the individuals are part of the study) and given ρ: We apply Bayes rule to the first of the three expressions on the right: P(y i = y j = 1 | s i = s j = 1; ρ) = P(s i = s j = 1 | y i = y j = 1; ρ)P(y i = y j = 1; ρ) P(s i = s j = 1; ρ) .
Under the full ascertainment assumption P(s i = s j = 1 | y i = y j = 1; ρ) = 1, and so P(y i = y j = 1 | s i = s j = 1; ρ) = P(y i = y j = 1; ρ) P(s i = s j = 1; ρ) . Similarly: and since a control is selected to the study with probability K(1−P ) P (1−K) , this boils down to: For the case of y i = y j , one individual is a case, and is automatically selected, while the other is a control and is selected with probability K(1−P ) P (1−K) . Hence: Using these results we get: Denote the numerator by A(ρ) and the denominator by B(ρ). We wish to approximate the last equation using a Taylor series around ρ = 0. Since the individuals are unrelated, the correlation is expected to be close to 0, and therefore such an approximation is expected to be good. Such an approximation would take the form: See Supp. Mat. for a discussion of a second order approximation.
Note that with ρ = 0, the phenotypes of the two individuals are i.i.d. and so A(0) = 0.
Therefore, the Taylor approximation can be simplified: Similarly, with ρ = 0 the events of being included in the study are i.i.d. for both individuals, so All that remains is to find A (0). Towards that end, we are interested in computing the probabilities of the three possible combinations of phenotypes: −∞ˆ∞ t f ρ,σ 2 g (l 1 , l 2 )dl 1 dl 2 , and P(y i = y j = 0; ρ, σ 2 where f ρ,σ 2 g is the multivariate Gaussian density, namely: with Σ denoting the covariance matrix of the liabilities, given explicitly by: The determinant of Σ is |Σ| = 1 − ρ 2 σ 4 g and its inverse is so the density function f ρ,σ 2 g can be written as: . Differentiating A(ρ) requires differentiating each of the three double integrals w.r.t. ρ: Setting ρ = 0 in the last expression yields: Explanation: we differentiate and set ρ = 0. By the chain rule, the derivative of any expression with ρ 2 is 0 at ρ = 0, and obviously the derivative of any expression which does not depend on ρ is 0. The only expression whose derivative is therefore not 0 at ρ = 0 is −2l 1 l 2 ρσ 2 g in the numerator of the exponent. The denominator of the exponent is 2 at ρ = 0, and so the derivative at ρ = 0 is Using these results we can write down A (0): and so: Hence, when the error of the approximation is small, the slope obtained by regressing Z ij on G ij is an unbiased estimator of P (1−P ) yields an unbiased estimator of σ 2 g -the liability scale heritability.
Extending the liability threshold model to include fixed effects. It is often desired to include fixed effects in the analysis of a complex phenotype. Such fixed effects might include external information such as sex, diet and exposure to environmental risks, but can also be genetic variants with known effects or estimates of population structure such as projections of several top principal components.
Since the liability threshold model is in fact a probit model, fixed effects can be included in the usual manner: where x i is a vector of the values of the relevant covariates and β is a vector of their respective effect sizes.
An individual is a case if l i > t, as before. However, an equivalent formulation would be to subtract the fixed effect from the threshold, rather than adding it to the liability: thus keeping the previous formulation of the liability as a sum of genetic and environmental effects.
Heritability estimation with known fixed effects. Assume first that the fixed effects are known, and so the t i 's are known. The probability of being a case, and of being included in the study, are no longer equal for all the observed individuals. We denote: the probability that the i'th individual is a case, and: the probability that the i'th individual is a case conditional on being selected for the study, where t i is computed using the fixed effects as described above. We now wish to derive the same first order approximation, while accounting for the newly introduced heterogeneity. We redefine: and follow the same steps as before (see supplementary materials for a full derivation). Using the first order Taylor approximation we conclude that the slope obtained by regressing Z ij on is an estimator of heritability on the liability scale.
Estimating heritability with unknown fixed effects. More often than not, the effects of relevant fixed effects are unknown and must be estimated from the data. However, estimating effect sizes under ascertainment in case-control studies is notoriously problematic. Specifically, under the threshold (probit) model, ignoring the ascertainment yields biased estimators.
A special exception is the case of logistic regression. In their seminal paper, Prentice and Pyke (1979) 26 proved that using a logistic regression to estimate fixed effects from ascertained data yields consistent estimators of these effects in the (unascertained) population, and that the ascertainment only biases the intercept.
We therefore suggest a two-step procedure for estimating heritability. First, we estimate the fixed effects using a logistic regression model. We then correct the effect of the ascertainment, and obtain the individual-specific thresholds. Lastly, we plug the thresholds into the estimation scheme described above.
More elaborately, by Bayes' formula: by the complete ascertainment assumption P(s i = 1 | y i = 1, x i ) = 1, and according to the selection scheme: We can thus solve for K i and express it is a function of P i : We then use logistic regression to obtainP i -a consistent estimator of P i , and use this estimate to obtain an estimate of K i , which is in turn used to estimate the threshold: Estimating the added variance due to fixed effects Lastly, the presence of fixed effects increases the variance of the liability, so σ 2 g no longer equals h 2 . The appropriate definition of heritability is now: where σ 2 t is the variance of the thresholds in the population, and so the estimate of σ 2 g can be transformed to an estimate of the heritability simply by dividing it by 1 + σ 2 t . We discuss how σ 2 t can be estimated from the data in the Supplementary material.

Simulations using a generative model
Our goal was to create a realistic simulation, covering a wide range of combinations of disease prevalence K, case sampling probability P and heritability σ 2 g , and for each one creating "natural" genotypes that realistically recreate the complex correlation structure induced by case-control sampling.
Given these parameters, our simulations proceeded as follows: 1. The MAFs of 10,000 SNPs were randomly sampled from U [0.05, 0.5].
2. SNP effect sizes were randomly sampled from N (0, 3. For each individual, we: (a) Randomly generated a genotype using the MAFs, and normalized it (according to Yang et al.'s model 27 ).
(b) Used the genotype and the effect sizes, to compute the genetic effect.
(d) Computed liability and phenotype.
(e) If the phenotype was a case -the individual was automatically included in the study.
Otherwise the individual was included in the study with probability K(1−P ) P (1−K) .

The genotypes of all included individuals were used to compute
This matrix G was then used to estimate heritability for the LMM using the GCT A software 27 and the correction of Lee et al. 2 ; and using GCR as described above.
We note that our choice to work with SNPs that are in linkage equilibrium was motivated by the analysis of Patterson et al. 23 . They show that for the purpose of generating correlation matrices, using SNPs in linkage disequilibrium (LD) is equivalent to using a smaller number of SNPs in linkage equilibrium. They also suggest a method for estimating the effective number of SNPs (i.e. the number of SNPs in linkage equilibrium leading to the same distribution of correlation matrices as a given set of SNPs in LD). Applying their method to WTCCC data suggests that the effective number of SNPs in linkage equilibrium is roughly one tenth of the actual number of SNPs. Hence, using 10,000 SNPs in our simulations is equivalent to simulating roughly 100,000 SNPs with realistic LD structure.

Heuristic correction for the LMM approach
Denote the estimate of the heritability on the liability scale, obtained through the method of Lee et al. 2 byσ 2 g , and the estimates obtained by our method byσ 2 g . Denote the true variance of the genetic effect under ascertainment by: We note a different analytical expression of σ 2 gcc is given in 2 . We validated correctness of our derivation of this expression numerically.
As detailed in the Supplementary Materials, our simulations demonstrated the following properties of these estimators: 1. They are both unbiased when there is no ascertainment.
2. Our estimateσ 2 g remains unbiased in all situations. However, in presence of ascertainment, σ 2 g is biased, and is not linear in the true σ 2 g for fixed K, P .
4. The bias ofσ 2 g worsens as the ascertainment factor K P decreases.
We therefore performed an extensive analysis of the relationship between K P and the bias of the "linearized" estimateσ 2 g σ 2 gcc σ 2 g in our simulations. Our analysis, as detailed in Supplementary Materials, led us to the following relationship betweenσ 2 g and the true underlying heritability σ 2 g : Correcting published results We have derived our heuristic correction using simulations wherein the true underlying heritability is known. However, contrary to previously used corrections, our correction is not a linear transformation of the estimate. When attempting to correct published estimates, we only know K, P andσ 2 g . We define our corrected estimateσ * 2 g to be the value of σ 2 g for which where the expectation is computed using the approximate relationship we derived previously. In other words, the corrected estimate of heritability is the value of heritability for which the expectation of the estimator is the observed estimate, where the expectation is calculated using our heuristic correction.    1 Heritability estimation using genetic correlation regression

Liability threshold model -notations
Denote K the prevalence of a condition in the population and P the prevalence in the study.
Under the liability threshold model, we assume that each individual i has an unknown liability l i = g i + e i where g i is a genetic random effect, which can be correlated across individuals, and e i is the environmental random effect, which is assumed to be independent of each other and of the genetic effects. Both effects are assumed to follow a Gaussian distribution with variances σ 2 g and 1 − σ 2 g respectively. A person is then assumed to be a case if her liability exceeds a threshold t = Φ −1 (1 − K), i.e. the phenotype y i is given by y i = I{l i > t}. This definition guarantees that the prevalence in the population is indeed K.

Selection probabilities
When the study is observational, the probability of being included in the study is independent of the phenotype. However, in a case-control study, the proportion of cases is usually greatly ascertained. To model this fact, we define a random indicator variable s i indicating whether individual i was selected to the study. The commonly used "full ascertainment" [11] assumption is P (s i = 1 | y i = 1) = 1. While this assumption can be relaxed, as discussed later, it simplifies subsequent analysis.
Suppose the population is of size n and that the expected size of the study is n s < n. The expected number of cases in the study is nK. Additionally, the proportion of cases in the study is P , so: nK n s = P, yielding: denote p the probability of a control being included in the study (i.e. p = P (s i = 1 | y i = 0)). The expected number of controls in the study is n(1 − K)p . Additionally, the proportion of controls in the study is (1 − P ) so: Solving for p yields: From here it follows that the probability of being included in the study for a given individual (with unknown phenotype) is:

Relaxing the "full ascertainment" assumption
In fact, the latter results do not depend strictly on the full ascertainment assumption (that is, P (s i = 1 | y i = 1) = 1). The latter assumption can be relaxed such that for any p * , as long as the probability of being selected as a control is multiplied by the same probability. This can model any step prior to the selection procedure, for example the probability that an individual is approached by the health administration to begin with. For example, a non-ascertained study might involve only a proportion p * of the population. In this case deriving the selection probability of a given individual yields p * as expected.

Heritability estimation
Next, consider a pair of individuals in the study, whose genetic effects are correlated and denote by ρ the correlation.
Denote by Z ij the product of the standardized phenotypes: The variable Z ij can obtain three values: We write down the expected value f Z ij , conditional on the fact that s i = s j = 1 (the individuals are part of the study) and given ρ: We apply Bayes rule to the first of the three summands on the right: Under the full ascertainment assumption P(s i = s j = 1 | y i = y j = 1; ρ) = 1, and so P(y i = y j = 1 | s i = s j = 1; ρ) = P(y i = y j = 1; ρ) P(s i = s j = 1; ρ) . Similarly: and since a control is selected to the study with probability K(1−P ) P (1−K) , this boils down to: For the case of y i = y j , one individual is a case, and is automatically selected, while the other is a control and is selected with probability K(1−P ) P (1−K) . Hence: Using these results we get: Denote the numerator by A(ρ) and the denominator by B(ρ). We wish to approximate the latter equation using a Taylor series around ρ = 0. Such an approximation would take the form: We later discuss a second order approximation as well.
Note that with ρ = 0, the phenotypes of the two individuals are i.i.d. and so A(0) = 0. Therefore, the Taylor approximation can be simplified: Similarly, with ρ = 0 the events of being included in the study are i.i.d. for both individuals, so B(0) = K 2 P 2 . All that remains is to find A (0). We are interested in computing the probabilities of the three possible combinations of phenotypes: −∞ˆ∞ t f ρ,σ 2 g (l 1 , l 2 )dl 1 dl 2 , and g is the multivariate Gaussian density, namely: with Σ denoting the covariance matrix of the liabilities, given explicitly by: The determinant of Σ is |Σ| = 1−ρ 2 σ 4 g and its inverse is and so the density function f ρ,σ 2 g can be written as: .
Deriving A(ρ) requires deriving each of the three double integrals w.r.t. ρ: Setting ρ = 0 in the last expression yields: Explanation: we differentiate and set ρ = 0. By the chain rule, the derivative of any expression with ρ 2 is 0 with ρ = 0, and obviously the derivative of any expression which does not depend on ρ is 0. The only expression whose derivative is therefore not 0 at ρ = 0 is −2l 1 l 2 ρσ 2 g in the numerator of the exponent. The denominator of the exponent is 2 at ρ = 0, and so the derivative is l 1 l 2 σ 2 Similarly: Using these results we can write down A (0): and so: Hence, when the error of the approximation is small, the slope obtained by regressing Z ij on G ij is an unbiased estimator of P (1−P ) yields an unbiased estimator of σ 2 g -the liability scale heritability.

Second order approximation
While the first order approximation yields very satisfactory results in our simulations, one can obtain a better estimator using a better approximation. The second term of the Taylor series takes the form: We have already derived: Now, B(ρ) is the probability that two individuals with genetic correlation ρ are included in the study: P (y i = y j = 0; ρ, σ 2 g ).

7
Using the previous derivatives of the double integrals yields: Computing A (0) requires computing the second derivative of the two-dimensional density at ρ = 0 : Hence: and the second order approximation can be written as: Since K and P are assumed to be known, the estimation problem boils down to a single-variable non-linear regression in σ 2 g .
2 Dealing with fixed effects

Extending the liability threshold model
It is often desired to include fixed effects in the analysis of a complex phenotype. Such fixed effects might include external information such as sex, diet and exposure to environmental risks, but can also be genetic variants with known effects or estimates of population structure such as projections of several top principal components.
Since the liability threshold model is in fact a probit model, these effects can be included in the usual manner: where x i is a vector of the values of the relevant covariates and β is a vector of their respective effect sizes.
An individual is a case if l i > t, as before. However, an equivalent formulation would be to substract the fixed effects from the threshold, rather than adding them to the liability: thus keeping the previous formulation of the liability as a sum of genetic and environmental effects.

Heritability estimation with known fixed effects
Assume the fixed effects are known, and so the t i 's are known. We define: and: to be the probability of the i'th individual being a case when accounting for fixed effects through the adjustment of the threshold t i .
We redefine: so now Z ij can obtain four possible values: We now wish to derive the same first order approximation before, while conditioning on the fixed effects.
Repeating the same steps as before while contitioning on t i , t j the expression for A(ρ) is now: Each phenotype was standardized to have mean 0, conditional on the relevant fixed effects. Additionaly, with ρ = 0 the phenotypes are independet and so the expected value of Z at ρ = 0 is 0. An immediate result is that A(0) = 0. To see this consider the Taylor expansion of but on the other hand: Therefore, A(0;ti,tj ) B(0;ti,tj ) = 0 and so A(0; t i , t j ) = 0. We conclude that the first order Taylor approximation is again of the form: To derive an explicit expression for A (0; t i , t j ) we need to differentiate the double integrals as before. When both i and j are cases, with thresholds t i , t j respectively, the double integral takes the form: differentiating it w.r.t. ρ and setting ρ = 0 yields: and differentiating the other double integrals yields similar results: and: We differentiate A(ρ; t i , t j ) w.r.t. ρ, and using the previous results and some algebra we get: As a sanity check, set P i = P j = P and t i = t j = t to get: which is the expression derived for the no fixed effects case.
Moreover, the the probability of inclusion of individual i in the study, conditional on the fixed effects is now: and so: Plugging the derived expressions for A (0; t i , t j ), B(0; t i , t j ), we conclude that regressing Z ij on yields an estimator of heritability on the liability scale.

Dealing with unknown fixed effects
More often than not, the effects of relevant fixed effects are unknown and must be estimated from the data. However, estimating effect sizes under ascertainment in case-control studies is notoriously problematic. Specifically, under the threshold (probit) model, ignoring the ascertainment yields biased estimators.
A special exception is the case of logistic regression. In their seminal paper, Prentice and Pyke (1979) [10] proved that using a logistic regression to estimate fixed effects from ascertained data yields consistent estimators of these effects in the (unascertained) population, and that the ascertainment only biases the intercept.
We therefore suggest a two-step procedure for estimating heritability. First, we estimate the fixed effects using a logistic regression model. We then correct the effect of the ascertainment, and obtain the individual-specific thresholds. Lastly, we plug the thresholds into the estimation scheme described above.
More elaborately, by Bayes' formula: by the complete ascertainment assumption P(s i = 1 | y i = 1, x i ) = 1, and according to the selection scheme: We can thus solve for K i and express it is a function of P i : We then use logistic regression to obtainP i -a consistent estimator of P i , and use this estimate to obtain an estimate of K i , which is in turn used to estimate the threshold:t i = Φ −1 1 −K i , and the estimates of the individual-wise thresholds are used for estimating the liability-scale heritability.

Estimating the added variance due to fixed effects
Lastly, the presence of fixed effects increases the variance of the liability, so σ 2 g no longer equals h 2 . The appropriate definition of heritability is now: where σ 2 t is the variance of the thresholds in the population, and so the estimate of σ 2 g can be transformed to an estimate of the heritability simply by dividing it by 1 + σ 2 t . Therefore, obtaining an estimate of the heritability in the presence of fixed-effects requires an estimate of σ 2 t in the population. To estimate σ 2 t we use the law of total variance: | y)) where y is the phenotype of the individual. Furthermore: and: Once we condition on a specific phenotype, both the expected and variance of T can be estimated from the data, as they are no longer affected by enrichment of cases in the data. Specifically, we use the estimated thresholdst i to estimate the expected value and the variance of the threshold for cases and controls, and plug the estimates into the equations above to obtainσ 2 t . We discuss how σ 2 t can be estimated from the data in the Supplementary material.

Description of the Lee et al. simulation scheme
Lee et al. perform simulations of ascertained case-control studies in the following manner: Phenotypes are simulated in blocks of 100 individuals Given the true genetic variance σ 2 g , liabilities for each are sampled from a multivariate Gaussian distribution with mean 0 and covariance matrix Σ = Gσ 2 g + I(1 − σ 2 g ), where G is given by: that is -individuals within each block have a genetic correlation of 0.05 while individuals in different blocks are perfectly unrelated. Individuals with liabilities higher than the threshold t are considered cases.
Cases are always included in the study. Controls are included in the study with probability K 1−K . Blocks are generated until 100 individuals are accumulated, and this process is repeated 100 times, so in total 10, 000 cases and controls are accumulated.

Discussion of problems with the simulations in Lee et
al.
There are several key aspects in which the simulations of Lee et al. differ from the true generative process of the data.
First, the resulting correlation matrix is highly degenerate, with most of the correlations being 0 -as much as 99.99\%. Since the underlying idea behind heritability estimation from unrelated individuals is leveraging the minor -but non-zero -correlations among unrelated individuals, restricting the simulations to largely degenerate correlation structures is highly unrealistic and counterproductive.
Second, in reality, correlations between individuals span a wide range of values. while in Lee et al.'s simulations the correlations are either 0 or 0.05. Hence negative correlations are impossible, and the expected correlation is larger than 0, contrary to the unrelated individuals assumption.
Thirdly, the selection procedure, in which cases are up-sampled compared to their prevalence in the population, results in cases being more similar genetically than controls. This is a well known phenomenon which has recently attracted considerable attention in the context of GWAS [14]. This is hardly captured by the simulation process of Lee et al. For example, cases from different blocks always have 0 genetic correlation between them in Lee et al.'s simulations, while in a realistic scenario they should have positive genetic correlation, whose magnitude depends on the heritability of the phenotype and the number of causative loci.
Lastly, for small enough prevalence, all blocks in a typical simulation would contain only a single individual. In this case the genetic effects would be completely interchangeable with the environmental effects since G = I and so heritability would be impossible to estimate. In reality, however, a very small prevalence should yield very closely related individuals, thus generating an intuitively easier estimation problem.
We also note that Lee et al.'s simulations use an unrealistic number of individuals -10,000.

Simulations using a generative model
To see if these problems have a major impact, we ran simulations using the full generative model. This was done as follows: 1. The MAFs of 10,000 MAFs were randomly sampled from U [0.05, 0.5].
2. SNP effect sizes were randomly sampled from N (0, 3. For each individual, we: (a) Randomly generated a genotype using the MAFs, and normalized it (according to Yang's model).
(b) Used the genotype and the effect sizes, to compute the genetic effect.
(e) If the phenotype was a case -the individual was automatically included in the study. Otherwise the individual was included in the study with probability K(1−P ) P (1−K) .
5. The genotypes of all included individuals were used to compute G = ZZ m where Z is the matrix of normalized genotypes.
6. We used GCTA [13] to estimate heritability with G as the genetic correlation matrix, and a 0 − 1 vector of the phenotypes as input.
We ran ten repetitions of this simulation for all combinations of σ    We then applied GCR to the same data. The results are portrayed in figure 3. As can clearly be seen from the figure, our estimators are linear in the underlying true heritability. Indeed, as can be seen in ????, applying the correction derived earlier yielded accurate estimators.

Distribution of estimates around the mean
Apart from being unbiased, GCR estimates display an equal or lower variance to the LMM estimates. This can be seen in figure 5.

Weak ascertainment
To test the robustness of our method to the different values of P and K, we ran a wide range of simulations, with less extreme ascertainment. The results are given in figure 6, showing how the downward-bias of the LMM estimates is smaller for less ascertained studies, and disappears completely when the studies are non-ascertained. in both cases the GCR yields unbiased estimates. True h^2 estimated h^2 (corrected) using matrix regression Figure 6: Comparing corrected LMM estimates (top) to GCR estimates (bottom) for studies with low-to-intermediate ascertainment.

non-ascertained simulations
We also applied our method to non-ascertained simulations, where both the LMM and GCR methods produce unbiased estimators, as can be seen in figure 7.

Simulations with population structure
To study the robustness of our method to the presence of subtle population structure, we modified our simulations to introduce population structure similar to the simulations in Lee et al. For every batch of 100 individuals, we randomly selected 5% of the SNPs and fix them for the entire batch. By doing so, the expected correlation within each batch is 0.05, as in Lee  The resulting GCR estimates were unbiased for all these scenarios, as can be seen in figure 8.

Simulations with different numbers of SNPs
The number of SNPs used in the simulations determines the variance of the realized genetic correlations between individuals. Since GCR relies on a first-order apporximation around zero correlation, high variance of genetic correlations might decrease the accuracy of GCR. To study the robustness of GCR to such changes, we reran our initial simulations using 1,000 or 100 SNPs. GCR still produced unbiased estimates. As expected, when the number of SNPs was smallest (100) and the heritability was highest (0.9), the results displayed cosiderably

Second order simulations
To study the accuracy of GCR, we implemented a version of GCR which utilizes the second order approximation as detailed before. We then compared the first-and second-order estimates of heritability on our initial set of simulated data. The resulting second order estimates were almost identical to the results obtained by the first order approximation (figure 10 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q

Simulations with fixed effects
To test that our method still generates accurate estimates of heritability even in the presence of fixed effects, we used the following simulation scenario: In addition to generating the genotypes as described earlier, for each simulated individual we generated a 0-1 "sex" variable with probability 0.5 for each value. We fixed the heritability to 0.5 and the proportion of cases in the study to 0.3, and defined the risk of individuals with sex = 0 to 0.005. The risk of individuals with sex = 1 was either 0.005, 0.01 or 0.02, i.e. the relative risk (RR) was 1, 2 or 4. The thresholds were adjusted to depend on the sex accordingly. Phenotypes were generated after accounting for the fixed effects by way of changing the threshold as described earlier, and individuals were selected for the study based only on their phenotypes (i.e. independently of their sex). We simulated 100 sets of genotypes/phenotypes/sex for each of the three possible RR values. We then estimated the heritability using LMM and GCR, each time with and without including the fixed effect in the analysis. The results are displayed in figure 11. Our simulations demonstrate that GCR correctly accounts for fixed effects while LMM estimates are biased with and without accounting for fixed effects.  Figure 11: Estimation of heritability with and without accounting for fixed effects using LMM and GCR. The true underlying heritability is 0.5 and the only method which is consistently accurate is GCR with fixed effects.

Inference
Standard inference for linear regression is not applicable in our case due to the breaking of several of the basic assumptions used in linear regression inference. For example, the errors are non-normal and not independent. We therefore use the jackknife [5] to estimate the variance of our estimates. Since our method is three to four orders of magnitude faster than LMM estimation, applying such a procedure is not considerably time-demanding. More specifically, we use the jackknife to estimate the standard deviation of the estimate, and use this standard deviation to construct a confidence interval (CI). At each jackknife iteration we remove one individual from the data, and remove all relevant entries to the genetic correlation matrix. We constructed 95% confidence intervals for our simulations and, these CIs covered the true value of the heritability in 98% of the simulations, indicating the jackknife is conservative in this case, and the true standard error of the estimates is usually smaller than the jackknife estimates.

Heuristic correction of LMM estimates
While we have demonstrated that our method can be used to yield accurate estimates of heritability, it might also be of interest to correct existing estimates of heritability which were obtained using LMM and so might be biased.
In their seminal work Dempster and Lerner [3] derive the relationship between additive heritability on the observed scale and heritability on the liability scale under the liability threshold model. However, D&L's work is done under a no-ascertainment assumption. Lee et al. [8] extend their result to ascertained case-control studies, and obtain the following relationship: Note that the relationship depends on the ratio of variances, so that the relationship is not linear. The variance of the genetic effects under ascertainment is given by: We note a different analytical expression of σ 2 gcc is given in Lee et al. We validated our derivation of this expression numerically. A script demonstrating the correctness of our derivation can be found on our webpage.
Despite the fact that this is a relationship between parameters, and we are interested in correcting estimates, we applied this correction to the estimates obtained in our simulations. Our reasoning was twofold -first, we hoped that the biases introduced by the estimation and the biases introduced by using the observed scale would be interchangeable, and thus correcting for the latter would expose the former, making its analysis easier. Second, from ??????????, it seems that LMM works well when heritability is low, regardless of K and P , and also that LMM works well when ascertainment is low. We speculated that this is due to σ 2 gcc being close to σ 2 g in such scenarios. The results of applying this correction to the simulations is shown in figure  12. As can be seen from the figure, the resulting estimates are still biased, except for the case of non-ascertained studies, but the dependency of the corrected estimates on the true underlying heritability seems linear. we speculated that the transformation indeed corrects for the effect of estimating under the observed scale (rather than under the liability scale), but do not correct for the additional effects introduced by the LMM estimation.
The resulting relationships between the corrected estimates and the true underlying heritabilities seemed linear, but with a slope which differs from 1, we speculated that the difference was due to the normality assumptions made by LMM, which are less accurate when ascertainment is strong. We therefore attempted to model the slope as a function of K and P . Since for the noascertained case K = P and the slope is 1, we speculated that the slope depends on K P , and that this size captures the degree of ascertainment. We therefore estimated the slope for each (K, P ) tuple and explored the dependence of the slope on the ratio K P . Using a Box-Cox plot [1] (figure 13) indicated that K P is the best transformation of the K P variable.  Figure 13: A Box-Cox plot using K P as the dependent variable and the slope as the covariate indicates that a square-root transformation is a good fit.
We then used linear regression to estimate the relationship: slope(K, P ) = 1.3 − 0.3 K P We concluded by applying this correction on top of the previous correction (i.e. dividing the estimates by the slope). The results of this correction are given in figure 14. While the results are not as well behaved as when using GCR, it seems that applying these two corrections one on top of the other reduces the bias considerably for all simulated ascertainment scenarios. While after correcting, the results are unbiased, in terms of variance, GCR still outperforms LMM

Applying the correction to published data
So far, we have derived our heuristic correction using simulations wherein the true underlying heritability is known. However, when attempting to correct published estimates, we only know K, P andσ 2 g LM M . We define our corrected estimateσ 2 g * to be the value for which E[σ 2 g LM M ; σ 2 g =σ 2 g * ] =σ 2 g LM M , i.e. the value of underlying heritability for which the expected value of the estimator is the observed estimate, where the expectation is calculated using our heuristic correction.
Confidence intervals are derived by applying the same procedure to the top and bottom limits of a 95% confidence interval based on the published standard deviation of the estimate.

Testing the correction on WTCCC data
To validate our suggested heuristic correction, and to make sure it generalizes beyond our simulations, we computed LMM estimates of heritability for all 7 WTCCC case-control studies, and applied the correction to the resulting estimate. We also applied our method to the data and compared the result to the corrected estimates. The results are highly correlated (ρ = 0.97) and can be seen in figure 16.

Applying our correction to other published results
We applied our correction to several published case-control studies [8,6,7,4]. The corrected heritability estimates are given in  [8,6,7,4] 6 Accounting for imperfect linkage-disequilibrium Yang et al. [12] suggest that the genetic correlation matrix, estimated based on the genotyped data, is in fact not the true underlying correlation matrix, as the latter should be computed only using causal SNPs. Since causal SNPs are unknown, and often not genotyped, the estimated genetic correlation is a noisy and biased estimate of the true genetic correlation. To correct for this effect, Yang et al. suggest using a different correlation matrix where G is the estimated genetic correlation matrix, and β is a constant correcting the bias between the true and estimated genetic correlations. Using extensive simulations using real data and realistic assumptions they estimate β at roughly 0.875. It is easy to show that using G * in LMM estimation yields an estimate of heritability which is 1 βĥ 2 , whereĥ 2 is the heritability estimate. In the case of GCR, we regress the product of normalized phenotypes Z ij on G ij . Hence, plugging in G * ij = βG ij as the covariate has the exact same effect of increasing the estimated heritability by a factor of 1 β .

Quality Control
Following Lee et al. (2011) we applied a stringent quality control (QC) process on the WTCCC data to avoid detecting spurious heritability due to genotyping differences between cases and controls or between the different control groups. We removed SNPs with MAF>5%, with missing rate >1% and SNPs which displayed a significantly different missing rate between cases and controls (p-value<0.05). We also removed SNPs which deviated from Hardy-Weinberg (HW) equilibrium in the control groups (p-value <0.05). Additionally we removed SNPs which displayed a significant difference in frequency between the two control groups. Sex chromosomes were excluded from the analysis. We removed all the individuals appearing in the WTCCC exclusion lists. These include duplicate samples, first or second degree relatives, individuals which are not of European descent and other reasons. In addition we removed individuals with missing rate >1% and all individual pairs with an estimated relationship >0.05 based on the correlation matrix. The last step is done to ensure individuals in the study are not closely related.

No evidence of population structure
Lee et al. [8] include several top principal component vectors as fixed effects in an attempt make sure heritability estimates are not inflated due to population structure in the WTCCC samples. We suggest that this procedure is not required for two reasons. First, we note that [2] perform a thorough investigation of population structure in the WTCCC samples, and conclude that once individuals of different ancestry are removed from the samples, the effect of population structure is very minute. Second, [9] present a straightforward statistical test to decide how many principal components should be included as covariates in an association study. However, their method assumes that the individuals are randomly sampled from the population. In our case, cases are over-sampled in the study, and so LD regions hosting a causative variant might be mistaken for population structure. One simple way to overcome this problem is to dilute the number of SNPs, making sure each LD block is represented by a small number of SNPs. It is then unlikely that the differences between the case and control groups due to caustive variants would be captured by a top ranking principal component.
We applied the method of [9] do the GWAS data using every k'th SNP for increasing values of k. If there's a real population structure, we expect the number of significant principal components to be non-zero asymptotically. However, this was not the case, in agreement with [2] (data not shown).
We note that including principal components which tag such regions hosting causative variants would result in an underestimation of the heritability, as such a vector would be highly correlated with the phenotype. It is therefore highly undesired to include principal components which do not tag actual population structure.