Bivariate causal mixture model quantifies polygenic overlap between complex traits beyond genetic correlation

Accumulating evidence from genome wide association studies (GWAS) suggests an abundance of shared genetic influences among complex human traits and disorders, such as mental disorders. Here we introduce a statistical tool, MiXeR, which quantifies polygenic overlap irrespective of genetic correlation, using GWAS summary statistics. MiXeR results are presented as a Venn diagram of unique and shared polygenic components across traits. At 90% of SNP-heritability explained for each phenotype, MiXeR estimates that 8.3 K variants causally influence schizophrenia and 6.4 K influence bipolar disorder. Among these variants, 6.2 K are shared between the disorders, which have a high genetic correlation. Further, MiXeR uncovers polygenic overlap between schizophrenia and educational attainment. Despite a genetic correlation close to zero, the phenotypes share 8.3 K causal variants, while 2.5 K additional variants influence only educational attainment. By considering the polygenicity, discoverability and heritability of complex phenotypes, MiXeR analysis may improve our understanding of cross-trait genetic architectures.

I n recent years, genome-wide association studies (GWASs) have successfully detected genetic variants associated with multiple complex human traits or disorders, providing important insights into human biology 1 . Understanding the degree to which complex human phenotypes share genetic influences is critical for identifying the etiology of phenotypic relationships, which can inform disease nosology, diagnostic practice, and improve drug development. Most human phenotypes are known to be influenced by multiple genetic variants, many of which are expected to influence more than one phenotypes (i.e., exhibit allelic pleiotropy) 2,3 . This has led to cross-trait analyzes, quantifying polygenic overlap, becoming a widespread endeavor in genetic research, made possible by the public availability of most GWAS summary statistics (p-values and z-scores) 4,5 .
Currently, the prevailing measure to quantify polygenic overlap is genetic correlation. For a pair of traits, polygenic overlap refers to the fraction of genetic variants causally associated with both traits over the total number of causal variants across the two traits considered, while genetic correlation quantifies the correlation coefficient of additive genetic effects for the two traits. The sign of the correlation indicates whether the shared genetic effects predominantly have the same or the opposite effect directions. Available methods can quantify genetic correlation using raw genotypes 6,7 or GWAS summary statistics [8][9][10][11] . However, these methods report overall positive, negative, or no genetic correlation, but fail to capture mixtures of effect directions across shared genetic variants. This scenario is exemplified by the genetic relationship between schizophrenia and educational attainment. Despite consistent estimates of a non-significant genetic correlation 12,13 , many genome-wide significant loci are found to be jointly associated with both phenotypes 14 . Among 25 shared loci 15 , 16 had effects in the opposite direction, while 9 had effects in the same direction. Thus, new statistical tools are needed to improve our understanding of the polygenic architecture of complex traits and their intricate relationships.
Here, we introduce a statistical tool (MiXeR), which quantifies polygenic overlap irrespective of genetic correlation between traits, using summary statistics from GWAS. To evaluate polygenic overlap between two traits, MiXeR estimates the total number of shared and trait-specific causal variants (i.e., variants with nonzero additive genetic effect on a trait). MiXeR bypasses the intrinsically difficult problem of detecting the exact location of causal variants, but rather aims at estimating their overall amount. MiXeR builds upon the univariate causal mixture model [16][17][18][19] , which we extend to four bivariate normal distributions as illustrated in Fig. 1, with two causal components for variants specific to each trait, one causal component for variants affecting both traits, and a null component for variants with no effect on either trait. From the prior distribution of genetic effects, we derive the likelihood function of the observed signed test statistics (GWAS z-scores), incorporating effects of linkage disequilibrium (LD) structure, minor allele frequency (MAF), sample size, cryptic relationships, and sample overlap. The parameters of the mixture model are estimated from the summary statistics by direct optimization of the likelihood function.
We show in simulations that MiXeR provides accurate estimates of model parameters in the presence of realistic LD

Results
Simulations studies. In our first set of simulations, we generated synthetic GWAS data that follow model assumptions and assessed the validity of MiXeR estimates (polygenic overlap, π 12 ; correlation of effect sizes within the shared polygenic component, ρ 12 , and genetic correlation, r g ) in the presence of a realistic LD structure (Fig. 2). We observed no bias in the estimates across a wide range of simulation scenarios (Supplementary Figs. 1-3), except for the scenario with low heritability (h 2 = 0.1) and high polygenicity (π u 1 ¼ 3 10 À3 ) in both traits, which represents an insufficiently powered GWAS study. In this scenario, MiXeR shows large variation among the estimates, and reports a nonzero polygenic overlap, when no such overlap exists. Standard errors estimated by the model are shown in Supplementary  Table 1. Similarly, we show that univariate estimates of polygenicity and heritability are correct in all scenarios (Supplementary Fig. 4, Supplementary Table 2), except for the same scenario with low heritability and high polygenicity. In this scenario, a large variation among the polygenicity estimates suggests that the biases can be explained by truncated distribution of the errors, as the polygenicity parameters π 1 , π 2 , and π 12 are bound to be nonnegative. Traits with low heritability should obtain correct MiXeR estimates when the GWAS are sufficiently powered.
In addition, we validated that the model accurately predicts GWAS quantile-quantile (Q-Q) plots ( Supplementary Fig. 5) and detailed Q-Q plots with SNPs partitioned into disjoint groups according to MAF and LD score ( Supplementary Fig. 6a, b). Detailed Q-Q plots show a stronger GWAS signal for SNPs with higher MAF and higher LD score. The model's prediction follows the same pattern, indicating that it correctly captures dependency of GWAS association statistics on MAF and LD score. To validate the accuracy of the predicted bivariate density, we generated conditional Q-Q plots ( Supplementary Fig. 7), showing observed versus expected −log 10 p-values in the primary trait as a function of significance of association with a secondary trait at the level of p ≤ 0.1, p ≤ 0.01, and p ≤ 0.001. We note that data Q-Q plots are closely reproduced by the model predictions across all p-value strata. Interestingly, scenarios without polygenic overlap show an enrichment, arising because GWAS p-values depend on allele frequency and LD structure, though this effect is generally smaller than enrichment arising due to shared causal variants.
Sensitivity analysis. For sensitivity analysis, we conducted simulations with traits that have a shared pattern of differential enrichment of heritability across genomic categories 20 , which are not accounted for by the MiXeR model. Simulations were informed by the enrichment pattern of schizophrenia 21 , as estimated by stratified LD score regression 22 (Supplementary Fig. 8a). In the univariate analysis, polygenicity was underestimated by about 20% (Supplementary Fig. 9), suggesting that the model likely groups adjacent causal variants together and interpret their combined effect as arising from a single variant. In the bivariate analysis, we observed a small upwards bias in the estimate of polygenic overlap (Supplementary Table 3), but it did not exceed 10% of the polygenicity across all sufficiently powered scenarios.
Another assumption in the MiXeR model is that effect sizes are independent of allele frequencies. To asses this assumption, we ran simulations with effect sizes drawn from the BayesS 19 model (see Online Methods). It characterizes MAF-dependent architecture in terms of a single parameter S, ranging from S = 0 (equivalent to the MiXeR assumptions) to S = −1 (equivalent to the LDSR assumption). We simulated three intermediate values, S = −0.25, S = −0.5, and S = −0.75, which cover the range of BayesS estimates observed for real GWAS data. The results ( Supplementary Fig. 10, Supplementary Tables 2, 4) highlight certain biases in MiXeR parameter estimates: for the extreme case of S = −0.75, heritability (h 2 ), univariate polygenicity (π u 1 ), and polygenicity of the shared genetic component (π 12 ) are underestimated by up to 25%; correlation of effect sizes (ρ 12 ) is overestimated by 25%; however, the genome-wide genetic correlation (r g ) appears to have no bias.
Finally, we ran simulations with an incomplete reference, and simulated phenotypes where causal variants were spread across our entire reference panel of N = 11,015,833 variants, but only a fraction (50, 25, or 12.5%) of the variants enter LD structure estimation and fit procedure. The results (Supplementary Table 5) show that the total number of causal SNPs, as well as the heritability, are estimated correctly, while the polygenicity parameter is different from the simulated value, because it reflects the fraction of all tagged causal variants with respect to the reference that went into LD structure estimation. GWAS summary statistics. We applied MiXeR to summary statistics from GWAS of 14 phenotypes, including five psychiatric 21,23-26 and four autoimmune 27,28 diseases, four anthropomorphic traits [29][30][31][32] , and educational attainment 12 (see Supplementary Table 6 for metadata about the studies). MiXeR estimates of genetic correlation (Table 1; Supplementary Table 7) were generally consistent with those of cross-trait LD Score Regression 8 , with the highest genetic correlation observed between schizophrenia and bipolar disorder. As expected, these disorders also exhibit substantial polygenic overlap, sharing 6.2 K out of 8.5 K causal variants involved. Here and below, the numbers of causal variants are reported as 22.6% of their total estimate, which jointly accounts for 90% of SNP heritability in each phenotype, to avoid extrapolating model parameters into the area of infinitesimally small effects (Supplementary Fig. 11).
Furthermore, MiXeR reveals important differences among traits with low genetic correlation, represented as Venn diagrams of shared and unique polygenic components ( Fig. 3; Supplementary Figs. 12, 13a-g). For example, schizophrenia and educational attainment exhibit substantial polygenic overlap, sharing 8.3 K out of 10.8 K of causal variants involved. On the contrary, schizophrenia and height share only about 0.8 K out of 10.6 K causal variants. Educational attainment and height also show low polygenic overlap, sharing 1.8 K out of 12.3 K causal variants. Nevertheless, these traits have a high correlation of effect sizes within the shared component, ρ 12 = 0.52 (0.04), which at genome-wide level is observed as genetic correlation of rg = 0.16 (0.01) according to MiXeR, or rg = 0.14 (0.01) according to LDSR.
MiXeR estimates of the unique polygenic components provide insights into trait-specific genetic architectures. For example, schizophrenia has 2.1 K causal variants not shared with bipolar disorder, less than 0.1 K variants not shared with educational attainment, but as many as 7.5 K variants not shared with height. Also, for the other phenotypes, the number of trait-specific causal variants varies across different pairs of traits (Fig. 3). Figure 4 and Supplementary Fig. 14a-g visualize the observed bivariate density of the GWAS-signed test statistics (z 1j , z 2j ), the predicted densityẑ 1j ;ẑ 2j from the MiXeR model, and the estimated bivariate density of the additive causal effects (β 1j , β 2j ) that underlie model predictions. Figure 4 gives real examples for the three different scenarios of polygenic overlap (genetically independent traits, polygenic overlap with and without genetic correlation, as shown in Fig. 1). Finally, we use conditional Q-Q plots 33,34 to compare the observed and predicted distributions of z-scores, and show that MiXeR-based prediction provides accurate estimates of the data Q-Q plots (Fig. 5), both in univariate and bivariate contexts.

Discussion
MiXeR is a statistical method for cross-trait analysis of GWAS summary statistics, which enables a more complete quantification of polygenic overlap than provided by the LD score regression method 8 . In addition to genetic correlation, MiXeR estimates the total number of shared and trait-specific causal variants, providing new information about the genetic relationships between complex traits and disorders.
MiXeR extends cross-trait LD score regression 8 by incorporating a causal mixture model [16][17][18][19] , thus relying on a biologically more plausible prior distribution of genetic effect sizes compared with the infinitesimal model 35,36 . We show that polygenicity, measured as a total number of causal variants, and discoverability 37 , measured as variance of additive genetic effects across causal variants, have major implications for the future of GWAS discoveries ( Supplementary Fig. 15).
Applying MiXeR to real phenotype data, we provide new insights into the genetic relationships between schizophrenia, bipolar disorder, educational attainment, and height. In line with the strong clinical relationship 38 between schizophrenia and bipolar disorder, and prior genetic studies 25,34,39,40 , we find substantial polygenic overlap between the two disorders. Other studies have reported more substantial genetic differences between the disorders 40 (albeit with strong correlation 40,41 ), likely because they did not specifically model the polygenic overlap. For Table 1 The results of cross-trait analysis with the MiXeR model for schizophrenia (SCZ), bipolar disorder (BIP), educational attainment (EDU) and height Trait  Columns: n 12estimated number of shared causal variants, reported in thousands; n 1 (n 2 )-estimated number of causal variants, unique to trait 1 (trait 2), reported in thousands; ρ 12correlation of effect sizes in shared polygenic component; rggenetic correlation (r g ¼ ρ 12 π 12 = ffiffiffiffiffiffiffiffiffiffi π u 1 π u 2 p , see Online Methods); rg LDSR -estimate of genetic correlation from LD Score Regression. The number of variants (n 12 , n 1 , and n 2 ) are adjusted to explain 90% of heritability in the corresponding component. Parameters are fitted using approximately 1.1 M HapMap3 SNPs  40 performed a combined GWAS of schizophrenia and bipolar disorder, and a differential GWAS of schizophrenia cases versus bipolar disorder cases. Our results indicate a higher polygenicity of schizophrenia than bipolar disorder, which is in line with the recent study by Bansal et al. 14 , who highlighted two schizophrenia subtypes. Our results also indicate that both schizophrenia and bipolar disorder have a small fraction of causal variants conferring disorder-specific risk (Fig. 3). Identifying shared and disorder-specific genetic variants is a subject of our future research, as it could provide critical knowledge about the distinct genetic architectures underlying these psychiatric disorders. Moreover, we find that nearly all causal variants influencing schizophrenia risk also appear to influence educational attainment, despite a genetic correlation close to zero (Table 1). This is in line with recent studies demonstrating shared genetic loci between schizophrenia and educational attainment 15 and a strong genetic dependence between the phenotypes possibly related to different subtypes of schizophrenia 14 . In contrast, while 89% of genetic variants influencing bipolar disorder also appear to influence educational attainment, there is in this case a significant positive genome-wide correlation of 0.191 (0.036), in agreement with the crosstrait LD score regression estimate of 0.188 (0.023) ( Table 1;  Supplementary Table 7). We show that polygenicity is best expressed as a total number of causal variants (Supplementary Table 5). Previous studies presented it as a fraction, which is highly dependent on the reference panel used (1.1 M hapmap in ref. 18 , or 484 K Affymetrix SNPs in ref. 19 ). When expressed as a total number, our estimates of polygenicity are consistent with previously reported results 19 (Supplementary Table 8). In addition, we estimate that under the assumptions of the MiXeR model, only 5% of causal variants are needed to explain 50% of heritability, and 22.6% of causal variants are needed to explain 90% of heritability (Supplementary Fig. 11). These numbers are expected to be less dependent on modeling assumptions, because with finite GWAS samples it is not possible to distinguish small effects from truly null effects. The actual number of causal variants is, potentially, even higher, as our model tends to clump together variants if they are in high LD with each other (Supplementary Tables 3, 4).
Some existing methods can already uncover polygenic overlap in the absence of genetic correlation. For example, conjFDR  33,34 is a non-parametric model-free approach, which detects shared genetic loci regardless of their allelic effect directions, by prioritizing variants with strong associations across more than one GWAS 42 . Other methods, including gwas-pw 43 and HESS 44 , also aim at detecting genomic loci jointly associated with two traits. MiXeR complements these methods by providing an easily interpretable high-level overview of the shared and unique genetic architectures underlying complex phenotypes. Among other notable methods for cross-trait analysis, the Gen-omicSEM 9 and MTAG 11 represent a multi-trait extension of the LD score regression. They can handle two or more traits at a time, but are based on the infinitesimal assumption, and quantify polygenic overlap using genetic correlation. For two-trait analysis, these methods are equivalent to LD score regression, thus we did not perform a formal comparison between GenomicSEM, MTAG, and MiXeR.
MiXeR has some notable advantages compared with the existing methods that implement causal mixture. First, our mathematical model for the likelihood term p(z j |β j ) is conceptually simpler and more flexible, resulting in unbiased estimates of model parameters across a wide range of simulation scenarios (Supplementary Figs. 1-3) and providing accurate prediction of GWAS z-scores across varying ranges of MAF and LD ( Supplementary Fig. 6a, b). Second, MiXeR implementation works well with a large reference of 10 M variants, while other methods have reduced the reference to 1.1 M HapMap SNPs (ref. 18 ) or 484 K Affymetrix SNPs (ref. 19 ). Finally, our model individually processes all SNPs, without grouping them into bins (ref. 16 ).
MiXeR models causal effects as a single Gaussian component, while recent work 18,45 suggests that certain phenotypes, including height, require at least two causal components of small and large effects. We note that the MiXeR model still provides a good fit for SNPs not reaching the GWAS significance threshold (Supplementary Fig. 16) and shows deviations only toward the tail of the distribution. To further investigate the effects of model misspecification, we implemented right-censoring of genome-wide significant SNPs (see Online methods). The results (Supplementary Tables 7, 9) are consistent with our main analysis, except for height which received a lower estimate of heritability (65% instead of 70%), a slight increase in polygenicity, and increased polygenic overlap with other traits. We propose that for a better estimate of height's polygenicity, it would be beneficial to run MiXeR on a   residualized GWAS, after covarying association statistics for genotypes of all genome-wide significant SNPs. Recent work suggests the importance of MAF-and LDdependent genetic architectures 19,46 , which are not directly modeled by MiXeR. Our simulations show certain biases in the estimates of polygenicity parameters (π u 1 and π 12 ), which underestimate the true value by up to 25% in the case of extreme MAFdependent architecture with S = −0.75 (Supplementary Fig. 10). However, these biases tend to cancel one another out when considering the relative size of the polygenic overlap (π 12 =π u 1 ratio). Also, on real data, most BayesS 19 estimates lay between S = −0.25 and S = −0.5, where the bias of the MiXeR estimates is in the order of 10% rather than 25%. On real data, we observe effects of MAF-dependent architectures by drawing Q-Q plots for subsets of SNPs ( Supplementary Fig. 17a-g) partitioned into nine groups according to MAF and LD score, where the model tends to underestimate z-scores in low MAF bins. This effect, however, is quite subtle, and does not manifest itself on the overall Q-Q plots ( Supplementary Fig. 16).
The MiXeR method requires large GWAS studies. Our recommendation is to apply MiXeR to studies with at least N = 50,000 participants, and inspect standard errors reported by MiXeR. In addition, MiXeR applies Bayesian information criterion (BIC) to compare causal mixture model versus the infinitesimal model, as shown in Supplementary Table 9. The cases where BIC selects the infinitesimal model indicate that the GWAS sample size is insufficient to reliably fit the polygenicity parameter. Generally, polygenicity estimation requires more GWAS power than heritability estimation, which can be visually explained by GWAS Q-Q plots ( Supplementary Fig. 16): heritability is determined by the overall departure of the GWAS curve from the null line, while polygenicity is determined by its curvature, i.e., the point where the GWAS curve begins to bend upward from the null line, which is harder to estimate when GWAS signal is weak. This is captured by MiXeR standard errors, which show that individual parameters of the mixture model have lower estimation accuracy than their combinations-for example, relative errors for π 1 and σ 2 β are larger than for the heritability estimate h 2 / π 1 σ 2 β , due to inversely-correlated errors (Supplementary Table 2 In our future work, we plan to incorporate an additional Gaussian component to model small and large effects 18 , and to explicitly account for MAF-dependent architectures 46 . Further extensions may account for differential enrichment for true associations across genomic annotations 20 . Another limitation to address is that the MiXeR model assumes similar LD structure among studies, and is not currently applicable for analysis across different ethnicities. We also aim to extend the MiXeR modeling framework to be used to improve power for discovery of shared and trait-specific SNPs by estimating the posterior effect size of SNPs associated with one trait given the test statistics in another trait, as well as for improving predictive power of polygenic risk scores.
In conclusion, MiXeR represents a useful addition to the toolbox for cross-trait GWAS analysis. By considering the intricate polygenic architectures of complex phenotypes, MiXeR allows for measures of polygenic overlap beyond genetic correlation. We expect this to lead to new insights into the pleiotropic nature of human genetic etiology.

Methods
Bivariate causal mixture model. Consider a simple additive model of genetic effects, ignoring gene-environment interactions, epistasis and dominance effects.
Under these assumptions, the contribution of the genotype to the phenotype is modeled as a sum of individual contributions from genetic variants: y k ¼ P j g jk β j , where y k is a quantitative phenotype or disease liability of k-th individual, g jk is 0, 1, 2-coded number of reference alleles for j-th variant, and β j is the additive genetic effect of allele substitution. We say that a genetic variant is causal for a trait if it has a non-zero effect on that trait (β j ≠ 0).
MiXeR builds upon the univariate causal mixture model 16 , β j $ π 0 N 0; 0 ð Þþπ 1 N 0; σ 2 β , which assumes that only a small fraction (π 1 ) of variants have an effect on the trait, while the effect of the remaining variants is zero. For the mathematical convenience, we chose a Gaussian distribution for the nonnull arm of the causal mixture. A drawback with the gaussian prior is that a large fraction of causal variants will have effect sizes close to zero. We would prefer to count a variant as causal only if it has a sufficiently large effect size, using for example a bi-modal prior distribution with probability mass separated from zero, but for such prior, it was not feasible to accurately model the effects of the LD structure.
In a joint analysis of two traits, we expect some variants to affect both traits; some variants to affect one trait but not the other; and most variants to have no effect on either trait. We assumed that for a given trait, all causal variants follow the same distribution of effect sizes, regardless of what effect a variant has on the other trait. Within the shared component, we model correlation of effect sizes, to account for genetically correlated traits. Based on these assumptions, MiXeR models additive genetic effects β 1j , β 2j of variant j on the two traits as a mixture of four bivariate Gaussian components (Fig. 1): where π 1 and π 2 are weights of the unique components (variants with an effect on the first trait only, and on the second trait only); π 12 is a weighting of the component affecting both traits; and π 0 is a fraction of variants that are noncausal for both traits, π 0 + π 1 + π 2 + π 12 = 1; σ 2 1 and σ 2 2 control expected magnitudes of per-variant effect sizes; and ρ 12 is the correlation coefficient of the effect sizes in the shared component. All parameters are assumed to be the same for all genetic variants.
The effectsβ 1j ;β 2j estimated by a GWAS, represent only proxies of the true causal effects (β 1j , β 2j ), which are distorted by limited sample size (poor statistical power), cryptic relatedness within a GWAS sample, as well as LD between variants.
To disentangle these effects, we derive the likelihood term for observed GWAS signed test statistics (z 1j , z 2j ), incorporating effects of the LD structure (allelic correlation r ij between variants i and j); heterozygosity H j = 2p j (1 − p j ), where p j is the minor allele frequency of the j-th variant; the number of subjects genotyped per variant (N 1j and N 2j ); and variance distortion parameters σ 2 01 , σ 2 02 , and ρ 0 . Specifically (see Supplementary Note 1), The nine parameters of the model (π 1 ; π 2 ; π 12 ; σ 2 1 ; σ 2 2 ; ρ 12 ; σ 2 01 ; σ 2 02 ; ρ 0 ) are fit by direct optimization of the weighted log likelihood, with standard errors estimated from the Observed Fisher's Information matrix.
Forcing π 12 = 1 (so that π 0 = π 1 = π 2 = 0) reduces our model to an infinitesimal assumption that underlies cross-trait LD score regression 8 . Under this constraint, our model predicts that GWAS-signed test statistics follow a bivariate Gaussian distribution with zero mean and variance-covariance matrix i.e., (z 1j , z 2j ) ∼ N(0, Σ j ), where ' j ¼ P i H i r 2 ij is the LD score (adjusted for heterozygosity). This model is consistent with cross-trait LD score regression, with expected chi-square statistics Eðz 2 1j Þ, Eðz 2 2j Þ, and cross-trait covariance E (z 1j z 2j ) being proportional to the LD score of j-th SNP, and parameters ρ 0 , σ 01 , σ 02 playing the role of LD score regression intercepts 47 . The only distinction here is that we choose to model effect sizes that are independent of allele frequency, leading to the incorporation of H i into our model; this factor is absent from the LD score regression model due to the assumption there of effect sizes that are inversely proportional to H i . Thus, MiXeR is a direct extension of cross-trait LD score regression, which relaxes the infinitesimal assumption.
Model for bivariate distribution of GWAS z-scores. We derive two models for GWAS z-scores, which we call "fast model" and "full model". The "fast model" is quicker to run, and we use it to perform an initial search in the space of the model's NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-10310-0 ARTICLE parameters. The "full model" is slower but more accurate, and we use it for a final tuning of model estimates.
The "full model" for GWAS z-scores approximates (z 1j , z 2j ) distribution of a given GWAS SNP as a mixture of K = 20,000 bivariate normal distributions, all having equal weight in the mixture. For each k = 1, …, K, we randomly draw the location of causal variants (π 1 N causal variants specific to the first trait, π 2 N specific to the second trait, and π 12 N shared causal variants, where N denotes the total number of variants in the reference panel), and calculate the variance-covariance matrix Σ ′ kj from equation (3), using estimated LD r 2 correlations between the assumed causal variants and the GWAS SNP. Then The "fast model" is derived from the method of moments (see Supplementary Note 1): where ⊕ denotes convolution of probabilistic distribution functions (so that right-hand side evaluates to a mixture of eight components), π ′ cj ¼ ' j π c =η j is adjusted weight of mixture component (c ∈ 1, 2, 12); ij is the LD score, adjusted for can be interpreted as shape parameter that affects fourth and higher moments of the distribution. This model h i of z-score distribution, and forms a theoretical basis for the mixture model of sparse and ubiquitous effects 49,50 . Of interest is that the "fast model" involves the forth power of allelic correlation r 4 ij , which is directly proportional to kurtosis (measure of heavy tails) of z-score distribution.
LD structure estimation. To estimate the LD structure, we use 489 individuals from the 1000 Genome project 51 (phase 3 data), obtained from the LD score regression website 8,22,52 . In total, 14 individuals were excluded due to relatedness 53 . For simulations, LD scores were estimated from the actual genotypes that we use to produce synthetic GWAS summary statistics. LD r 2 coefficients were calculated using PLINK 54 with LD r 2 cutoff of 0.05 and fixed window size of 50,000 SNPs, corresponding on average to a window of 16 centimorgans. We deliberately chose a larger LD window compared with the LDSR-recommended window of 1 centimorgan, because the later appears to truncate a noticeable part of LD structure. At the same time, we did not observe an effect of using an unbiased estimate 55 of r 2 , thus fall back to the standard Pearson correlation coefficient. We employ small integer compression 56 for efficient storage of the LD matrix.
Fit procedure. We fit the model by direct optimization of weighted log likelihood where θ ¼ π 1 ; π 2 ; π 12 ; σ 2 1 ; σ 2 2 ; ρ 12 ; σ 2 01 ; σ 2 02 ; ρ 0 À Á is a vector of all parameters being optimized, and weights w j chosen by random pruning (64 iterations at LD r 2 0.1). Optimization is done by the Nelder-Mead Simplex Method 57 as implemented in MATLAB's fminsearch. First, we fit univariate parameters separately for each trait, i.e., π u 1 ; σ 2 1 , σ 2 01 for the first trait, and similarly for the second trait. Univariate fit employs a sequence of optimizations to ensure robust convergence: first, we use the "fast model" under constraint π u 1 ¼ 1 to find σ 2 1;inf and to initialize σ 2 01 ; second, we use constraint π u 1 σ 2 1 ¼ σ 2 1;inf to find initial values of π u 1 and σ 2 1 , again with the "fast model". Finally, we use the "full model" and unconstrained optimization to jointly fit π u 1 ; σ 2 1 , σ 2 01 parameters. The same procedure is repeated for the second trait, to find π u 2 ; σ 2 2 , σ 2 02 . To improve convergence, we parametrize univariate log-likelihood as a function of log π u 1 σ 2 1 À Á and log π u 1 =σ 2 1 À Á , which represent almost independent dimensions of the energy landscape. In bivariate optimization, we use the "fast model" and constraint π 12 = 1 to estimate r g and ρ 0 . Then, we proceed with the "full model" optimization of the parameters specific to the bivariate model (π 12 , ρ 12 ), constraining all other parameters to their univariate estimates, and also constraining r g and ρ 0 to the estimates from the infinitesimal model. The additional analysis (Supplementary Tables 7, 9) uses right-censoring 58 of zscores exceeding z t = 5.45, by using cumulative distribution function 59 in the log likelihood: Standard error estimation. We estimate standard errors of all parameters from the observed Fisher's information, based on the "fast model". It is known from the likelihood optimization theory that the observed Fisher's information may not be suitable for a parameter near its boundary, which is applicable to the mixture weights π 1 , π 2 , π 12 , and the correlation of effect sizes ρ 12 . To mitigate this problem, we apply transformations-MATLAB's logit () for π 1 , π 2 , π 12 , exp() for σ 2 1 ; σ 2 2 ; σ 2 01 ; σ 2 02 , and erf() for ρ 0 , ρ 12 , and estimated a variance-covariance matrix of errors in the transformed parameter space. We validated that our estimates based on the observed Fisher's information are in good agreement with block jack-knife estimates. To estimate standard errors for a function of the parameters, such as r g or h 2 , we incorporate linear correlation among parameter errors in the transformed space. We sample N = 1000 realizations of the parameter vector, calculating the function (e.g., r g or h 2 ) on each of them, and report the standard deviations. In cases when joint hessian was not positive definite, we estimate marginal errors of fitted parameters.
Akaike and Bayesian information criteria. We apply standard formulas and estimate AIC = 2k − 2F and BIC ¼ k ln n À 2F, where F is the log-likelihood from Eqs. (6) or (7), k is the number of parameters (k = 2 for an infinitesimal model, and k = 3 for causal mixture model), and n ¼ P j w j is the effective number of SNPs (the sum of weights across all SNPs used to fit the model).
Large LD blocks. The log-likelihood cost function and the Q-Q plots apply a weighting scheme to SNPs to 0avoid overcounting evidence from large LD blocks.
As an alternative to weighting by inverse LD score, we chose to infer the weights by random pruning. This technique is a stochastic procedure which averages loglikelihood function across repeatedly selected subsets of variants, such that for each pair of variants i, j in a subset J the squared allelic correlation r 2 ij falls below a certain threshold. Given T iterations of random pruning the log-likelihood function can be calculated as follows: which is equivalent to weighted log-likelihood F θ ð Þ ¼ P j w j log pdf z j jθ with weights w j = |{t:j ∈ J t }|/T, |S| denotes cardinality of set S. Random pruning with stringent threshold r 2 = 0.1 justify independent modeling of the residuals in Eq. (3) across SNPs, which otherwise would be correlated.
Heritability estimates. In an additive model, SNP heritability is defined as a sum across all causal variants: σ 2 β P j:β j ≠0 2p j 1 À p j , which we approximate from an average heterozygosity of all variants in the reference: π 1 H total σ 2 β , where H total ¼ P j 2p j 1 À p j . To estimate the proportion of causal variants that explain a certain fraction of heritability ( Supplementary Fig. 11), we randomly sample N = 10,000 causal effects from the reference, draw their effects β j from normal distribution, sort according to β 2 j p j 1 À p j , and report the fraction of variants that cumulatively account for 90% of heritability.
Genetic correlation. Parameter ρ 12 in MiXeR defines the correlation of effect sizes within the shared polygenic component. Genome-wide genetic correlation, calculated across all SNPs, is related to ρ 12 by the following formula that involves polygenicity π u 1 ¼ π 1 þ π 12 and π u 2 ¼ π 2 þ π 12 of the traits, and polygenic overlap π 12 : For traits with K-fold difference in polygenicity (π u 1 ¼ Kπ u 2 ), the formula predicts an upper bound on genome-wide genetic correlation: r g ρ 12 = ffiffiffi ffi K p , where equality holds if causal variants of the less polygenic trait form a subset of the higher-polygenic trait.
Quantile-quantile plots. Univariate Q-Q plots and stratified Q-Q plots for the model were constructed from pdf j (z) density as defined by Eq. (3), given fitted parameters of the model and LD structure of j-th SNP, calculated across a fine grid of z-scores ranging from 0 to 38 with 0.05 step. We average pdf j (z) across 1% of randomly sampled SNPs, and numerically integrate the resulting probability density function to convert it into a cumulated distribution function. Error bars on data Q-Q plots represent the 95% binomial confidence interval q ± 1:96 , where q is the probability of observing a p-value as extreme as, or more extreme then the chosen p-value, and n total is the effective number of SNPs after controlling for LD structure, which in our case was calculated as a sum of random pruning weights across all SNPs.
GWAS power curves. Causal mixture model can project the future of GWAS discoveries, by estimating proportion S(N) of narrow-sense heritability captured by genome-wide significant SNPs at a given sample size N. The S(N) is defined as follows: where z t = 5.45 gives z-score corresponding to the standard genome-wide significance threshold 5 • 10 −8 , and C(z, N, j) ≡ P(z, N, j) ⋅ E(δ 2 |z, N, j) denotes a posterior effect size E(δ 2 |z, N, j) of the non-centrality parameter δ 2 for a GWAS SNP j, given certain z-score, multiplied by a prior probability of observing that zscore. Probability density function P(z, N, j) is given by Eq. (4), and E(δ 2 |z, N, j) can be calculated from the Bayesian rule. Thus, C z; N; j ð Þ¼ R δ 2 P zjδ ð ÞP δ; j ð Þdδ, where zjδ $ N δ; σ 2 0 À Á . Analytical expressions for C(z, N, j) and R z:jzj!z t Cðz; N; jÞdz are given in the Supplementary Note 1.
SNPs in the analysis. To enable a direct comparison of our model with LD score regression, we use the same set of SNPs in our log-likelihood optimization, which consists of approx. 1.1 million variants, subset of 1000 Genomes and HapMap3 60 , with MAF above 0.05, ambiguous SNPs excluded, imputation INFO above 0.9, MHC and other long-range LD regions excluded. Calculation of the LD structure, LD scores ' j and shape parameter η j are based on 9,997,231 SNPs from 1000 Genomes Phase 3 data, downloaded from LD score regression website. In simulations we generate GWAS and estimate LD structure on a subset of 11,015,833 SNPs from 1000 Genomes Phase 3, with MAF above 0.002, call rate above 90%, excluding duplicated RS numbers; the fit procedure was constrained to~130 K GWAS SNPs, keeping only HapMap3 SNPs, and pruning SNPs at LD r 2 threshold of 0.1.
LD score regression estimates. For dichotomous phenotypes, we used an effective sample size of N eff = 4/(1/N case + 1/N cont ) to account for imbalanced numbers of cases and controls, both in MiXeR and in LD score regression. In addition, we ran LDSR using MiXeR MAF model (using --per-allele flags in LD score estimation), and show the results alongside with original LDSR estimates (Supplementary Tables 7, 9). For case/control phenotypes heritability is reported on the observed scale.
Simulations. In our simulations, we use a panel of N = 100,000 samples and 11,015,833 SNPs, generated by HapGen2 61 using 1000 Genomes 51 data to approximate the LD structure for European ancestry. To avoid relatedness across individuals, we run HapGen2 for small disjoint chunks of about 2900 SNPs at a time, 3920 chunks in total. The chunks were acting as additional recombination hotspots, causing certain changes in the distribution of the LD scores (Supplementary Fig. 20). However, the total amount of allelic correlation in the HapGen2 panel was still substantial, for example the median LD scores in the HapGen2 panel was 66.4, versus 63.5 in the 1000 Genomes panel, which makes the HapGen2 panel appropriate for our simulations.
We also validated that the HapGen2 panel shows no signatures of cryptic relatedness and sample stratification. The "plink --pca" analysis of the genotype matrix shows no signatures of sample stratification, as shown by the scatter plot of the first and second principal components (Supplementary Fig. 20). The "plink --genome" test found no related individuals (PI_HAT measure was below 0.1 for all pairs of individuals). We use a subset of 115,267 SNPs in the analysis, selected according to steps described in the PCA module of the Ricopili GWAS pipeline.
For each simulation run, we use PLINK to obtain GWAS summary statistics, including Wald's test z-score and p-value, of two synthesized quantitative phenotypes, with complete sample overlap between GWAS samples. Quantitative phenotype y k of k-th sample is calculated via a simple additive genetic model, y k ¼ P j g kj β j þ ϵ k , where g kj is the number of reference alleles for j-th SNP on k-th sample, β j is causal effect size, and ϵ is the residual vector drawn from normal distribution with zero mean and variance chosen in a way that sets heritability h 2 = var(Gβ)/var(y) to a predefined level.
For the simulations with differential enrichment, we simulate three levels of polygenicity (3 K, 30 K, and 300 K causal variants), three levels of heritability (0.1, 0.4, and 0.7), and for each combination, generate 20 pairs of genetically independent traits (except for having shared pattern of enrichment). To simulate the enrichment, we keep a constant variance of effect sizes across all SNPs, but modulate the probability of having causal variant proportionally to LDSR regression coefficient. We use --per-allele flags in LD score estimation to run simulations with the MiXeR MAF model.
For simulations with MAF-dependent architectures, we simulate effect sizes as follows: Parameter S = 0 corresponds to the MiXeR MAF model, S = −1 corresponds to the LDSR MAF model. The same model is implemented in BayesS software 19 , thus we choose parameter S from −0.25, −0.50, and −0.75, which corresponds to the range of BayesS estimates observed on real GWAS data.