{'fraction_non_alphanumeric': 0.03925586756748103, 'fraction_numerical': 0.01552944190067561, 'mean_word_length': 4.591031149301826, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 26, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Contrasting genetic architectures of schizophrenia and other complex diseases using fast variance components analysis

Heritability analyses of GWAS cohorts have yielded important insights into complex disease architecture, and increasing sample sizes hold the promise of further discoveries. Here, we analyze the genetic architecture of schizophrenia in 49,806 samples from the PGC, and nine complex diseases in 54,734 samples from the GERA cohort. For schizophrenia, we infer an overwhelmingly polygenic disease architecture in which ≥71% of 1Mb genomic regions harbor ≥1 variant influencing schizophrenia risk. We also observe significant enrichment of heritability in GC-rich regions and in higher-frequency SNPs for both schizophrenia and GERA diseases. In bivariate analyses, we observe significant genetic correlations (ranging from 0.18 to 0.85) among several pairs of GERA diseases; genetic correlations were on average 1.3x stronger than correlations of overall disease liabilities. To accomplish these analyses, we developed a fast algorithm for multi-component, multi-trait variance components analysis that overcomes prior computational barriers that made such analyses intractable at this scale.


Exact methods
By "exact methods" we mean methods that locally converge to the REML optimum (vs. methods that approximate the REML optimum, which we discuss next). All REML procedures are iterative and compute analytically exact solutions only in the limit of infinite iterations; additionally, all iterative optimizers are at least somewhat sensitive to their starting points and are not guaranteed to find the global optimum from an arbitrary starting point.
First derivative methods. These methods perform an analog of gradient descent and thus are most robust-as they monotonically increase the likelihood and never leave the valid parameter domain-but converge slowly.
• Expectation maximization (EM). Dempster et al. [66] is the canonical reference; see Sec. 8. 3 of Searle's text [67] for an in-depth discussion and derivation of update equations. The popular GCTA software [2] offers EM as a non-default option.
• EM with parameter expansion (PX-EM). Liu et al. [68] first proposed PX-EM as a general approach to speed up EM convergence; Foulley and van Dyk [69] provide a comprehensive treatment of PX-EM for LMMs.
Second derivative methods. These methods use exact or approximate second derivative information about the log likelihood surface and therefore converge much faster than EM when near the optimum but are less robust when far from the optimum. Consequently, many implementations start by taking a few EM steps (GCTA) or PX-EM steps (GEMMA [70], PX-AI [71]). Second derivative methods can be classified by the way in which the Hessian is computed or approximated.
• Newton-Raphson (NR) (exact Hessian). Newton's method is the standard derivative-based optimization procedure; in the context of REML, Zhou and Stephens [70] provided an important computational advance for NR REML by developing an efficient eigendecompositionbased algorithm (implemented in the GEMMA software) to compute the gradient and Hessian of the log likelihood.
• Method of scoring, a.k.a. Fisher scoring (approximate Hessian). This approach replaces the negative Hessian, which is the observed Fisher information, with its expected value, which is the information matrix (independent of the observation y) and is sometimes easier to compute than the Hessian [67]. GCTA offers Fisher scoring as a non-default option.
• Average information (AI) (approximate Hessian). Gilmour et al. [26] proposed using a matrix equal to the average (appropriately defined) of the negative Hessian (used in Newton-Raphson) and the Fisher information (used in the method of scoring), terming this matrix the "average information" as it is an average of the observed and expected information. This matrix is computationally convenient because it eliminates a trace term involving an inverse matrix; however, AI REML still requires O(N 3 )-time matrix operations. GCTA uses AI REML as its default option.
• Broyden's method (approximate Hessian, quasi-Newton). The basic idea is to avoid direct computation of the Hessian by approximating it using a multi-dimensional analog of the secant method. Groeneveld [72] applies Broyden's method to REML in genetics.

Monte Carlo methods
Unlike exact methods, Monte Carlo REML methods use random sampling to approximate derivatives of the log likelihood, with the goal of trading off a small amount of accuracy in favor of greater computational speed.
Existing aproaches. • MC "NR"/AI/Broyden. To our knowledge, Matilainen et al. [23] were the first to use Monte Carlo to directly approximate second derivative-based methods-which offer far faster convergence than EM-in the context of LMM REML, where the computational challenge is problem size rather than model complexity. They investigated three algorithms that all use Monte Carlo sampling to approximate the gradient of the log likelihood and differ in their Hessian approximation methods: • Trust region optimization using approximate change in log likelihood. We achieve robust convergence by applying a trust region method to detect and prevent large, detrimental steps when far from the optimum. (Trust region methods are among the most commonly used numerical methods for robust multi-parameter optimization [77].) The basic idea is to identify where the local quadratic model of the log likelihood (based on its approximate gradient and Hessian) starts breaking down by comparing model-predicted vs. computed values of the log likelihood. Here, computing the log likelihood is infeasible, so instead we apply a gradient-based approximation of the change in log likelihood.
• Simple convergence criterion via reuse of Monte Carlo samples. By reusing pseudorandom Monte Carlo sample components across iterations, the approximate, Monte Carlo gradient we compute varies smoothly with variance parameter estimates. (The AI matrix is exactly computed and does not use Monte Carlo, hence is always a smooth function of the variance parameters.) Consequently, our MC AI REML iteration truly converges rather than jumping around near the optimum (as is typically the case for Monte Carlo parameter optimization methods), allowing easy assessment of convergence. Garcia-Cortes et al.
[25] did take such an approach in their MC EM REML algorithm, but Matilainen et al.
[23] did not, instead suggesting more complex convergence criteria that have been developed for general MC EM methods.
• Efficient software implementation. We implement our MC AI REML algorithm in a welloptimized, open source software package. In particular, BOLT-REML uses a computational framework requiring no O(M N ) or O(N 2 ) memory allocations beyond genotype data storage and performs all rate-limiting computations using multithreaded BLAS Level 3 matrix operations.
When specialized to the single-trait, single-variance component case, BOLT-REML is very similar to the heritability estimation procedure implemented in BOLT-LMM [20]; we expect that our use of secant iteration in BOLT-LMM vs. the trust region approach in BOLT-REML is unlikely to reduce robustness in the univariate case.

Single-trait case
We begin by developing the mathematics of Monte Carlo REML for one phenotype. We will later generalize to multiple phenotypes, but the multi-trait notation is more cumbersome, so we start with the single-phenotype case to convey the intuition behind the approach.

Model and log likelihood
We consider the model [21] y = σ 0 u 0 + σ 1 Z 1 u 1 + · · · + σ K Z K u K , in which y denotes phenotype, σ 1 Z 1 u 1 , . . . , σ K Z K u K are K random genetic effects, and σ 0 u 0 is a random environmental (noise) effect, where σ 0 , . . . , σ K are variance parameters to be estimated, the Z k contain normalized SNPs corresponding to the k-th non-identity variance component, and all entries of all the u k are iid standard normal. (Note that the u k are vectors of different lengths: u k has length equal to the number of SNPs in the k-th variance component. Also note that for the case of a single non-identity variance component containing M SNPs with normalized genotypes X, we may write Z = X/ √ M and u = √ M β, so that β ∼ N (0, σ 2 g /M ) ⇐⇒ u ∼ N (0, σ 2 g ), Xβ = Zu, and XX /M = ZZ is the GRM.) We have not included fixed effect covariates in the model for notational convenience, but covariates can be incorporated by projecting them out of all components and adjusting dimensionality accordingly (the REML approach [24]).
We wish to find the values of the variance parameters σ 2 0 , σ 2 1 , . . . , σ 2 K that maximize the likelihood of observing the phenotype y (a vector of observations over N individuals) under the above model. We can write the model equivalently as y ∼ N (0, V ), where V = σ 2 0 I N + σ 2 1 Z 1 Z 1 + · · · + σ 2 K Z K Z K . 6 The pdf of this multivariate normal distribution is given by where N is the number of samples (which changes to N − C in REML with C independent covariates projected out). Thus, up to a constant that we ignore, the log likelihood of observing y assuming parameters σ 2 0 , σ 2 1 , . . . , σ 2 K is given by = (σ 2 0 , σ 2 1 , . . . , σ 2 K ) = − 1 2 (log det V + y V −1 y).

Monte Carlo approximation of the gradient of the log likelihood; MC REML
We may compute the gradient and Hessian of the log likelihood wrt σ 2 k by using the matrix calculus formulas For large sample sizes, building and inverting V is computationally challenging, leading us to rewrite the trace term in terms of an expectation that does not involve matrix inverses, following equation (14.12) of ref. [78]: where y V is a random vector such that For fixed σ 2 k , sampling such a vector y V is easily done: we build y V = σ 2 0 u 0 + σ 2 1 Z 1 u 1 + · · · + σ 2 where the components of all the u k vectors are iid standard normal. We note in passing that were V a linear combination of variance components with some negative coefficients, we could still apply the above trick by writing V = E[y V + y V + ] − E[y V − y V − ] with y V + and y V − built from terms corresponding to the positive and negative variance components, respectively. The first-order conditions for optimality of the variance parameters are the vanishing of the gradient of the log likelihood, i.e., ∂ ∂σ 2 k = 0. As noted in ref.
[78], σ 2 k Z k V −1 y is the BLUP estimate u k of u k , so the terms appearing in equation (1) are expected and observed squared norms of BLUP estimates, and REML optimization amounts to equating these quantities. We previously presented the rudimentary heritability estimation algorithm used within the BOLT-LMM mixed model association software [20] in this way. Similarly, the parameter updates of EM REML are scaled differences between observed and expected squared norms of BLUP estimates (Sec. 4); comparing equations (1) and (13)(14) show that EM REML is very closely related to gradient descent.
In practice, we estimate the gradient of the log likelihood by replacing the expectation in equation (1) with an empirical average over a finite number S of Monte Carlo samples of y V : where the overline denotes a Monte Carlo average over S samples. For large N , typically only a few samples are required to obtain precise enough estimates of these derivatives to use within an optimization routine. Computing the quantities on the right only requires solving the mixed model equations, which can be done efficiently with conjugate gradient iteration [18,19]. Moreover, by reusing the same S sets of pseudorandom components u 0 , . . . , u K used in equation (2)  becomes a smooth function of the variance parameters σ 2 k , which is convenient for numerical optimization. We refer to the approach of approximately maximizing the log likelihood by solving ∇ = 0 as MC REML.
We note that while we have an efficient Monte Carlo method for estimating the gradient of the log likelihood, this method does not estimate the log likelihood itself, which involves a log determinant. We are not aware of an efficient way to estimate the log determinant that applies here; ref. [79] proposes an approximation based on the expansion where x ∼ N (0, I) are Monte Carlo samples, but this method requires the eigenvalues of A to be between −1 and 1 and appears unsuitable here.

Average information (AI) approximation of the Hessian
We may similarly express the second-order partials (i.e., the entries of the Hessian) as Instead of approximating this matrix using Monte Carlo sampling, we use the average information (AI) approximation of ref.
[26]. This approximation amounts to replacing the random y V in the expectation term of equation (4) with the observed phenotype y, giving

Multiple traits
We now generalize to multiple traits, following ref. [22,70,80]. To avoid notational clutter, we assume a model with only one genetic effect term, but all of the following easily generalizes to multiple random genetic effects.

Model
Consider D correlated traits, each measured at the same N individuals. We denote the "stacked" DN -vector of trait values by y and denote the N -element subvectors corresponding to individual traits by y d , d = 1, . . . , D, so that We assume the linear random effects model where Z is an N × M coefficient matrix (e.g., Z = X/ √ M under our usual SNP normalization), u g is a stacked DM -vector of effect sizes, and u e is a stacked DN -vector of noise residuals such that .
The distributional assumption on u g says that for every effect index m, the d-vector of effect sizes for the d traits is distributed as N (0, V g ), and across different m, effect sizes are independent of each other. An analogous statement holds for the environmental terms u e . Marginalizing across individual traits, each trait satisfies the usual LMM: We treat the genetic and environmental D × D covariance matrices V g and V e , each of which have D(D + 1)/2 independent coordinates (due to symmetry), as unknowns to be estimated. In the general case of K random genetic effects, each random effect as well as the environmental term has its own set of D(D + 1)/2 covariance parameters, so that the model has a total of P = (1 + K)D(D + 1)/2 parameters.
We note that fixed effect covariates need to be included independently for every trait. In particular, for D traits, D "all-1s" covariate vectors need to included, namely, the indicator vectors corresponding to the coordinates of each of the D traits. In REML analysis, doing so is equivalent to projecting covariates out from the genotypes and each phenotype independently (and adjusting dimensionality accordingly).

Monte Carlo gradient approximation
We denote the log likelihood of y assuming variance parameters V g , V e by = (V g , V e ) as before, and we combine the distinct entries of the covariance matrices V g , V e (D(D + 1)/2 entries each) into the parameter vector θ = (θ p ), where p = 1, . . . , P indexes parameters. Thus, To compute a Monte Carlo estimate ∇ of the gradient ∇ = ∂ ∂θp , we replace the expectation with an empirical average over S Monte Carlo samples of y V ∼ N (0, V). Explicitly, we generate y V via y V = (I D ⊗ Z)u g,rand + u e,rand where u g,rand = vec(randn(M, D) · chol(V g ) ), u e,rand = vec(randn(N, D) · chol(V e ) ), so that To make the Monte Carlo gradient ∇ a smooth function of θ, we generate S fixed sets of matrices randn(M, D) and randn(N, D), one per Monte Carlo trial, and reuse these matrices for all parameter values θ, as described in Sec. 2.1.2.

Average information (AI) approximation of the Hessian
As in equation (4), we write to approximate the Hessian using the negative AI matrix. (Note that the second-order partials of V do not appear in equation (6) because under our parameterization, ∂V/∂θ p is a constant matrix for each θ p , so that second-order partials conveniently vanish.) The AI approximation has two significant advantages over Monte Carlo estimation of the Hessian. First, it is guaranteed to be positive semidefinite, leading to more robust numerical optimization of the log likelihood, especially when some covariances are close to singular. Positive semidefiniteness can be seen by horizontally concatenating the P vectors ∂V ∂θp V −1 y into a DN × P matrix in which case we may write Because V is psd (i.e., has a Cholesky factorization), it follows that the same holds for I A . Second, computing I A is cheaper than using Monte Carlo: evaluating the expressions in equation (7) only requires solving the mixed model equations once per variance parameter θ p .

Standard errors
In general, standard errors of maximum likelihood parameter estimates can be obtained by approximating the covariance matrix of the estimates with the inverse of the negative Hessian. (This approximation amounts to fitting a quadratic form to the log likelihood, noting the vanishing of its gradient at the MLE. The negative Hessian is the observed Fisher information and measures the curvature of the log likelihood.) Approximating the negative Hessian with the AI matrix, we have The usual (non-Monte Carlo) AI REML estimated standard error of θ p is the square root of the p-th diagonal entry of the above covariance. For MC AI REML, we multiply the standard error by a factor of 1 + 1/S to take into account the additional variance introduced because MC AI REML essentially fits the observed data to a reference set of S simulated data sets.

Standard errors under transformed coordinates
It is most convenient (for exposition and implementation) to parameterize variance components using entries of the covariance matrices as we have done above: as noted earlier, the (θ p ) parameterization results in simple formulas for ∂V ∂θp and the vanishing of second-order partials of V. However, the quantities we typically wish to report are heritabilities (i.e., fractions of total variance contributed by random effects) or genetic correlations (i.e., ratios of off-diagonal covariance parameters to square roots of products of diagonal entries). Computing point estimates of these quantities is trivial once maximum likelihood variance parameter estimates have been found, but estimating their standard errors requires transforming the Hessian of the log likelihood according to the desired reparameterization, which requires a little more work. We derive the transformation formulas below; once we have the Hessian of the log likelihood with respect to transformed coordinates, we can then estimate error covariances by the usual approach of inverting the negative Hessian.
Consider a scalar function f (x) (where x is vector-valued); denote its Hessian by H x . Suppose we wish to reparameterize in terms of coordinates X via a transformation τ : In general, applying the product rule to ∇ x ((∇ x f ) J) gives both a Hessian term ∇ x (∇ x f ) J = H x J and a term involving ∇ x f and derivatives of J. However, at a critical point of f (such as the MLE if f is a log likelihood), the gradient vanishes, leaving the simple formula We also note that the AI matrix defined in equation (7) obeys the same transformation formula under change of coordinates. We can derive the transformation by observing that for a reparameterization φ → θ, the partial derivatives of V transform according to The partials of θ p with respect to the new parameters φ q are entries of the Jacobian, and expanding gives the transformation I A → J I A J.
In BOLT-REML, we implement the following two sets of transformations to reparameterize from covariance parameters θ p (i.e., entries of V g and V e ) to genetic correlations and heritabilities. We carry out the first set of transformations independently across all variance components (environmental and genetic), after which we carry out the second set independently across all traits.
Transformation from covariances to correlations. For each variance component, we have D(D + 1)/2 independent variables that parameterize its covariance matrix. We wish to transform the D(D − 1)/2 off-diagonal covariance parameters into correlations: that is, we start with D(D + 1)/2 parameters σ ij , 1 ≤ i ≤ j ≤ D (ignoring the index of the variance component under consideration for ease of notation), and we wish to reparameterize using D(D − 1)/2 parameters ρ ij = σ ij / √ σ ii σ jj , 1 ≤ i < j ≤ D, along with the D unchanged diagonal parameters σ ii , the nonzero entries of the Jacobian are given by where J r c denotes the r, c entry of the Jacobian.
Transformation from variances to heritabilities. For each trait, we have 1 + K independent variables σ 2 0 , σ 2 1 , . . . , σ 2 K that parameterize the variances of the environmental and genetic variance components for that trait (ignoring the index of the trait under consideration for ease of notation). We wish to reparameterize with 1 overall variance parameter σ 2 and K heritability fractions h 2 k , 1 ≤ k ≤ K, such that The Jacobian is given by 2.3.2 Breakdown of SE estimates at parameter space bounds, e.g., zero heritability In the above discussion of standard errors and transformed SEs, we need to be careful that if the MLE variance parameter vector is located on the boundary of the parameter space due to domain constraints, the gradient of the log likelihood may not equal zero. In the single-trait case, this situation typically arises when a variance parameter would have been estimated to be negative if negative variances were allowed, but is instead estimated to be zero because of the nonnegativity constraint. An analogous statement holds in the multi-trait case with covariances required to be positive semidefinite.
If the gradient does not vanish at the MLE due to bound constraints, then the quadratic form approximation of the log likelihood, which underlies the analytic standard error from equation (8), breaks down. (Similarly, the transformation formula (9) for change of coordinates no longer holds.) Therefore, analytic standard errors should only be trusted for variance components not at boundary constraints, and we further recommend that in situations with zero-heritability variance components, REML should be rerun with those variance components removed.

Distinction between SE and likelihood ratio test for nonzero variance
We also note that the analytic standard error of a heritability parameter provides only an approximate test of whether the corresponding variance component explains nonzero heritability (by computing z = σ 2 /s.e.( σ 2 )): such a test assumes a quadratic approximation to the log likelihood, i.e., a normal distribution of the posterior density of the heritability, which in reality is only locally valid. The proper way to test for nonzero heritability is to perform a likelihood ratio test; however, a limitation of MC REML is that it does not efficiently compute likelihood ratio test statistics (as it achieves speedup by circumventing likelihood computations).
In particular, for variance components of very low rank (e.g., modeling a small subset of SNPs), we have observed that computing a naive z-score based on the analytic standard error produces an extremely underpowered test for nonzero contribution to variance: even if the component explains a very large fraction of variance (reflected in an appropriately large REML parameter estimate), the analytic standard error is also very large, so that the z-score achieved is only ≈2. The reason is that the posterior distribution on the variance parameter is highly non-normal and asymmetric. For example, consider the case of a variance component containing just one SNP estimated to explain 1% of variance. Despite this huge effect, REML cannot be sure whether the true variance parameter is actually 2% (and the SNP just happened to get a smaller effect size under the random effects model). Thus, REML will report a large standard error (reflecting the uncertainty in the upper tail)-even though a closer examination of the likelihood actually would show that at the lower tail, the variance explained is convincingly >0% (as the likelihood ratio test would show).
For the goal of computing efficient likelihood ratio tests for whether low-rank variance components contribute nonzero variance-which is distinct from our primary goal of estimating variance parameters-we recommend the spectral decomposition approach of FaST-LMM-Set [81], which performs fast set association tests using two random-effects components, one for the SNP set being tested and the other to model potential genome-wide confounders. (In order to achieve computational efficiency, the latter component can only use a small subset of SNPs, so precautions such as including principal components as fixed effect covariates may be necessary [82]. We further note that for case-control ascertained data sets-in which ascertainment induces correlations among causal SNPs-modeling only a small subset of SNPs should help to avoid loss of power when testing SNP sets for associations [27]; in contrast, when estimating heritability contributed by a SNP set, all SNPs need to be modeled in order to avoid inflation of estimates due to the same phenomenon.)

BOLT-REML algorithm
We now describe our algorithm to compute approximate REML estimates of variance parameters. The overall structure of the algorithm is derivative-based iterative optimization using a trust region method for robust convergence. We begin with an initial guess θ = (θ 1 , . . . , θ P ). In each iteration, we propose an update step p = (∆θ 1 , . . . , ∆θ P ) such that θ + p achieves the largest predicted increase in log likelihood based on a local quadratic model of the log likelihood (built from the Monte Carlo approximate gradient and AI approximate negative Hessian), subject to parameter domain and trust region constraints. We evaluate the approximate gradient at the tentative new parameter vector θ + p and use its value to approximate the actual change in log likelihood. We compare the approximate actual change in log likelihood to the predicted change based the local quadratic model, and based on this information, we decide whether to accept the step and/or modify the trust region radius. We declare convergence of the algorithm when the approximate change in log likelihood falls below a specified threshold (corresponding to the distance of the current iterate from the optimum in units of standard errors).

Single trait
Because of the possibility of local optima, a good initial guess is very important to the robustness of an optimization routine. In the single-trait case with K non-identity variance components, we can use the sample variance of the phenotype as a rough guess of the overall variance parameter σ 2 , leaving K heritability fractions to create guesses for as follows: 2. For each variance component k, estimate the total heritability h 2 −k explained by all variance components except k using one step of secant iteration (as in the original BOLT-LMM reestimation of heritability for LOCO reps [20]; the only difference here is that we leave each variance component-rather than each chromosome-out in turn).

Multiple traits
To obtain guesses of the parameters along the diagonal terms of the covariances V g , V e , we apply the single-trait initial guess procedure above to each trait in turn. To obtain guesses of the offdiagonal covariance terms, we compute for each pair of traits the phenotypic correlation between the two traits. We then assume that the environmental correlation and genetic correlation(s) (plural in the case of multiple genetic variance components) all equal the estimated phenotypic correlation.

Trust region iterative optimization
After making an initial parameter vector guess, we proceed to take steps in parameter space designed to iteratively optimize the log likelihood. At each iteration, we need to take into account two important considerations: (1) optimizing the step to try to increase the log likelihood as much as possible (according to the local quadratic model), to achieve fast convergence; but also (2) being careful not to take too ambitious of a step, to ensure robustness. We balance these considerations by using a trust region approach that enforces an adaptive bound ∆ on the maximum allowed size of the step (with coordinates scaled to improve conditioning).

Step computation
In each iteration, we begin with a current parameter vector estimate θ and wish to compute a step p = (∆θ p ) to an updated estimate θ + p. We choose the step to maximize the predicted change in log likelihood according to the local quadratic model (where ∇ is the Monte Carlo gradient and I A is the AI negative Hessian) over steps p satisfying the following constraints: 1. REML parameter domain: Positive semidefiniteness of the covariance matrices (one per variance component) corresponding to the new parameter vector θ + p.
The value of ∆ bounds the region in which the local quadratic approximation of is trusted. Note that the objective function in equation (10) is a simple quadratic, so that if there were no constraints, the optimal p would be the Newton step. On the opposite extreme, if ∆ is very small, the optimal p is approximately parallel to ∇ as in gradient descent. Finding the step p that maximizes the predicted change in log likelihood given in equation (10) subject to the above bounds is itself an optimization problem, but because the objective function is quadratic, we can apply standard numerical optimization methods to solve it at trivial computational cost. To increase robustness, we run three optimization routines implemented in NLopt

Step acceptance criterion and trust region radius update
After computing a potential step p using the above procedure, we need to decide whether to accept the step and also whether and how to update the trust region radius ∆. Usually, trust region methods do so by comparing the predicted change in the objective function (here, the log likelihood) to the actual change, but a challenge that our method faces is that it can only estimate derivatives of the log likelihood: we are unable to compute or estimate the log likelihood itself. Instead, we use the following derivative-based proxy for the change in the objective function: (Note that this quantity is different from the predicted change in log likelihood given in equation (10), which was what we optimized to compute the step.) The rationale is that the change in log likelihood equals the dot product of the step vector p with the mean value of the gradient along the step from θ to θ + p. We can approximate the mean gradient with the average of the (approximate) gradients at the beginning and end of the step, i.e., at θ and θ+p, and this approximation is precisely equation (11). The approximation is an accurate estimate of the change in the log likelihood as long as the approximate gradient is reasonably close to linear over the step from θ to θ + p, which in particular holds if the step size is not too big. We update the trust region radius ∆ according to how well the change in log likelihood predicted by the quadratic model (10) matches the approximate actual change (11). We use the procedure and parameters recommended by ref.
[77] with the additional safeguard that we reject steps that more than double the norm of the gradient (suggesting that the local quadratic model may be invalid). Details are provided in the pseudocode (Sec. 3.4).

Convergence criterion
We declare convergence of the iterative optimization when the predicted change in log likelihood drops below a set threshold (10 −4 by default). Note that suboptimality of the achieved vs. optimal log likelihood has a simple interpretation in terms of the accuracy of the final parameter estimates. If our converged parameter estimates θ final give log likelihood final whereas optimal REML estimates θ opt give log likelihood opt , then if we assume that the log likelihood is locally well-approximated by the quadratic −I A (which is the assumption underlying the analytic standard error estimates discussed in Section 2.3). If we further assume near-independence of coordinates, we have

Number of Monte Carlo samples and speedup via two-phase optimization
In choosing the number of Monte Carlo samples S to use in estimating the gradient of the log likelihood, we need to balance two considerations: (1) accuracy, which increases with S (because the variance of our estimates scales roughly with 1 + 1/S, as discussed in Sec. 2.3); and (2) computational cost, which increases linearly with S plus a constant (Sec. 3.6). In our BOLT-REML implementation, we finesse this trade-off by applying a two-phase optimization using different values of S. First, we run a "coarse" optimization using S = 15 samples and stopping the iteration once the predicted change in log likelihood drops below 0.01. Then, starting from the final iterate of the coarse optimization, which should be close to the true optimum, we run a "fine" optimization using S = 100 samples and iterating to a log likelihood tolerance of 10 −4 .

Computational cost
The cost of evaluating the approximate gradient in equation (3) is dominated by computation of the terms V −1 y V for each of S Monte Carlo simulated phenotype vectors y V and the term V −1 y (for the real phenotype vector y). This computation amounts to solving a batch of S + 1 mixed model equation systems, which we do using conjugate gradient iteration [18,19]. To compute the AI matrix entries given in equation (7), we need to further compute V −1 ∂V ∂θp V −1 y for each parameter θ p , requiring P = (1 + K)D(D + 1)/2 additional solutions of the mixed model equations.
Thus, the total cost of the algorithm per AI iteration is O(S + KD 2 ) conjugate gradient runs, each of which has cost scaling roughly with M N 1.5 D 1.5 according to ref.
[20]. (As defined above, S is the number of Monte Carlo samples, K is the number of variance components, D is the number of traits, M is the total number of SNPs across all variance components, and N is the number of individuals.) In our experience, fewer than 10 AI iterations are typically needed for convergence.

Comparison to other approaches
For context, we now provide a detailed discussion of the computational requirements of BOLT-REML vs. existing approaches. To streamline the exposition, we restrict to the single-trait case D = 1, which is most frequently used in practice, and we focus on the scaling of computational costs with N and M . We primarily compare BOLT-REML to GCTA [2], which is the most widely used and most versatile among publicly available REML software for genetics; we point out other methods that contribute algorithms of particular interest.
Computation of genetic relationship matrices (GRMs). For variance component models built from genetic random-effect terms, most REML algorithms (e.g., GCTA) begin by computing the genetic relationship matrices corresponding to the random effects. (BOLT-REML and the FaST-LMM [86] family of spectral decomposition-based methods are exceptions.) Computing the GRMs takes O(M N 2 ) time, which may appear to be the dominant term in typical settings with M > N ; however, this computation can easily be parallelized across a computing cluster (e.g., with PLINK2 [60]), so in practice, it is not typically a computational barrier.
Computation per REML iteration. For large N , the running time of standard Newton-type REML methods is dominated by the costs (per REML iteration) of inverting the estimated phenotypic covariance (i.e., computing V −1 ) and performing a few N × N matrix multiplications per entry of the Hessian. Both inversion and matrix multiplication are O(N 3 ) operations, and the inversion is particularly difficult to parallelize across multiple machines. BOLT-REML iterations have cost scaling roughly with M N 1.5 (which is also difficult to parallelize across machines), so for typical M < 10 6 in genetics, BOLT-REML achieves a speedup once N exceeds several thousand individuals, with the break-even point increasing with M .
We note that in the special case of only one non-identity variance component, spectral decomposition techniques drastically reduce computational cost per-REML iteration [70,87]. We also note that for situations in which all or all but one of the variance components have low rank, two other approaches exist that can (for certain ranges of M and N ) achieve speedup beyond BOLT-REML. Listgarten et al. [81] observed that if the total number of SNPs in the model is smaller than the sample size (i.e., M < N ), then spectral decomposition can be used to reduce computational cost to O(M 2 N ) (for one SVD) per REML iteration. FaST-LMM-Set [81] implements this approach for mixed model set association tests. Speed et al. [88] observed that if only one variance component has full rank (as in our estimation of per-Mb schizophrenia SNP-heritability), then its GRM only needs to be computed and inverted once at O(M N 2 + N 3 ) cost, after which the Woodbury matrix identity can be used to efficiently perform REML iterations at O(M other N 2 ) cost, where M other is the total number of SNPs in the low-rank components. MultiBLUP [88] implements this approach for mixed model prediction.
Memory use. Standard REML methods compute and work with one N × N GRM per variance component; GCTA currently stores all of these matrices in memory, so that in practice, memory constraints become a computational barrier before running time. For models with many variance components, the memory footprint could be reduced by rereading each GRM from disk when necessary; however, O(N 2 ) memory would still be needed to work with at least a few N ×N matrices. In contrast, BOLT-REML (like BOLT-LMM [20]) stores only raw genotypes in memory, requiring ≈M N/4 bytes of memory (regardless of the number of variance components). We have found this memory requirement to be mild in practice, but it could potentially be reduced by storing genotypes in compressed form and decompressing on-the-fly in each REML iteration (which would especially help for rare SNPs) or by rereading genotypes from disk in each REML iteration; the latter approach would incur a substantial performance penalty, however.

Appendix: EM REML
Here we derive the BLUP estimates (assuming known variance parameters) and EM REML update equations (for variance parameter estimation) for the multi-trait mixed model given in equation (5)

Mixed model equations: BLUP estimates
To solve the mixed model equations and obtain BLUP estimates of random effects, we write out the joint distribution of y, u g , and u e . From the model equations (5) Conditioning on y, we have (by Appendix S.3 of ref. [67]) Therefore, the BLUP random effect estimates given observed y arê

EM REML update equations
Following the derivation of ref.
[67] Sec. 8.3 (noting in particular equation (14) on p. 300), the EM equations for estimating V g and V e amount to the successive updates for p, q = 1, . . . , D, where the expectations are taken over random effect estimates based on randomly generated y according to the assumed model parameters V g and V e . Note that at convergence, the pairwise dot products of estimated random effects (resp. residuals) across all pairs of traits equal their expectations.
Detailed derivation of EM update equations. For p, q = 1, . . . , D, we wish to update V g,pq and V e,pq according to the sufficient statistics We will carry out the derivation for the genetic terms; the residual terms behave analogously. Write using equation (12) for the covariance conditional on y. Now      Figure 4. Probability of a zero heritability estimate as a function of actual heritability explained per region: Robustness to genetic architecture. Plotted data points (×) correspond to empirical data from simulations; solid curves are best-fit sums of two exponentials. The red curve is the same as the red curve in Supplementary Fig. 3a, which we used in our main per-1Mb heritability analyses of schizophrenia ( Figure 2). The simulations used to generate this curve used genetic architectures in which 2K-470K causal SNPs were assigned effect sizes drawn from a normal distribution. The black curve corresponds to an alternative simulated genetic architecture in which effect sizes were drawn from a Laplace distribution ( Supplementary Fig. 6). The green and blue curves combine the simulation data from the Gaussian and Laplace effect size distributions but restrict to simulations with ≤20K causal SNPs (green) or ≥30K causal SNPs (blue). All four curves are nearly identical, indicating that the probability of a zero heritability estimate depends almost solely on the actual heritability explained per region (independent of how the heritability is distributed across SNPs in the region). We simulated three case-control data sets with case over-ascertainment according to the following procedure. We simulated genotypes by generating individuals as mosaics of up to 100 random "ancestors" from the PGC2 controls, resampling ancestors every 250 SNPs. We simulated case-control phenotypes using a liability threshold model in which we first generated continuous phenotypes with h g 2 =0.24 explained by either 4,000, 20,000, or 100,000 markers and then defined cases as individuals with phenotypes exceeding a threshold corresponding to 0.4% population risk of disease. In each simulation, we ascertained 22,537 cases and 28,109 controls, matching the case-control ratio of the PGC2 schizophrenia data we analyzed. (The exact numbers of cases and controls differ slightly from our PGC2 data because they match an earlier round of QC.) (a) For each simulated case-control data set, we estimated h g,1Mb 2 per 1Mb region and computed the fraction P (0) of regions with zero h g,1Mb 2 estimates (solid blue horizontal lines), which we compared to the calibration data we used in our schizophrenia analysis (dashed red curve, as in Fig. 2b) based on quantitative traits simulated for real PGC2 individuals. We observed that the calibration curve closely matches the number of simulated causal SNPs in each simulation. (b) We then applied the same procedure we used in our main analyses ( Fig. 2c and Online Methods) to bound the fraction of SNP-heritability explained by the 1Mb regions explaining the most heritability (dashed red curves), which we compared to the true fractions of SNP-heritability explained by the top 1Mb regions (solid blue curves), known from the simulation procedure. We observed that the our bounding procedure produced a conservative upper bound in each of the three simulations, as expected for a (conservative) 95% confidence bound. 2 equal to its lower bound constraint of zero in schizophrenia (solid) and simulated phenotypes with varying degrees of polygenicity and with h g 2 matching the h g-cc 2 of schizophrenia (dashed). Simulation data plotted are means over 5 simulations; error bars, 95% prediction intervals assuming Bernoulli sampling variance and taking into account s.e.m. The 95% error bar overlaps the observed fraction of zero h g,1Mb 2 estimates for genetic architectures involving >30,000 SNPs. (b) Conservative 95% confidence intervals for the cumulative fraction of SNP-heritability explained by the 1Mb regions that contain the most SNP-heritability. Lower bounds are from a cross-validation procedure involving only the schizophrenia phenotype while upper bounds are inferred from the empirical sampling variance of h g,1Mb 2 estimates (Online Methods).    Figure 9. SNP-heritability of schizophrenia liability partitioned by per-1Mb GC content after stratifying by meQTLs. We partitioned 1Mb regions into 3 × 5 = 15 categories as follows. First, we partitioned the regions into approximate meQTL tertiles (sets containing 1Mb regions with 0-1 meQTLs, 2-4 meQTLs, or 5+ meQTLs). We then partitioned the 1Mb regions in each approximate tertile into GC content quintiles, producing a partition of SNPs into 15 meQTL×GC variance components, on which we performed heritability partitioning using BOLT-REML. Panels (a-c) display the results subdivided by approximate meQTL tertile. The enrichment of heritability with higher GC content appears to be driven by regions in the lowest meQTL tertile.  We performed weighted least-squares regression of estimated h g,MAF 2 for schizophrenia in each of six MAF bins (Fig. 4b) using h g,MAF 2 from simulations (based on a given value of the parameter α) as a lone regressor (with no intercept term) and using regression weights equal to the inverses of estimated h g, MAF 2 variances. For α = 0, -0.25, -0.5, -1, we averaged h g,MAF 2 over 4,000 simulation replicates; for other values of α, we interpolated h g,MAF 2 from the closest two simulated values of α. We plot the mean-squared residual error as a function of α; the value of α minimizing residual error is -0.28. We applied a leave-one-chromosome-out jackknife to estimate the standard error of this estimate and obtained an s.e. of 0.09.  Fig. 4b,c for the three GERA traits with largest h g-cc 2 ; see the caption of Figure 4 for details.) We ran GERA analyses using the 558,778 GERA SNPs with MAF≥2% and ran UK10K simulations (to infer h MAF 2 from h g,MAF 2 ) using the 535,626 SNPs with IDs also present in the UK10K data. We averaged simulation results over the 91% of GCTA runs that converged. Estimates of h MAF 2 are based on best-fit α = -0.47 (0.25), 0, -0 Supplementary Figure 14. Impact of including covariates on SNP-heritability of GERA traits. The heat map depicts the ratio of h g 2 estimated using additional covariates to h g 2 estimated using only our standard set of covariates (age, sex, 10 principal components, and Affymetrix kit type). The first nine rows correspond to augmenting the covariate set with each survey-derived covariate in the GERA data individually; the last row corresponds to including all nine additional covariates. The "general health" covariate denotes self-reported health status at date of survey completion.  Figure 5, (a) before adjusting for BMI as a covariate, and (a) after adjusting for BMI as a covariate. For the regression fits, we used r l as the lone regressor (with no intercept term); we note that standard errors on values of r l are negligible (Supplementary Table 15). Horizontal and vertical bars indicate estimates of r l and r g plus or minus 1 s.e.  Explicitly, for each individual, we simulated genotypes as mosaics of chromosomal segments sampled from 100 random "ancestors" chosen from its ancestry group, resampling every 500 SNPs. We then generated genetic components of disease liabilities using a polygenic model (20,000 causal SNPs with total h g 2 =0.27; details described below). The simulations in Supplementary Fig. 16 did not include ancestry effects (i.e., these simulations contain population structure but not stratification), whereas the simulations in Supplementary Figures 17-19 include population effects explaining 5% of variance in liability (thus inducing stratification). We then generated remaining variance in liability from independent normal distributions. We designated individuals with liabilities in the top 1% as cases and sampled 22,177 cases and 27,629 controls in this way, matching our PGC2 data.
To test the robustness of heritability partitioning using REML under case-control ascertainment, we simulated genetic architectures in which SNPs have per-allele average effect size variances proportional to p(1 − p) α (chrNum) chrPow , where p denotes minor allele frequency, chrNum denotes the chromosome (1-22) on which the SNP is located, and α and chrPow are parameters.
We then partitioned heritability across MAF bins (a, c, e) and across chromosomes (b, d, f) using BOLT-REML with 10 principal component covariates and applying the liability threshold transformation of ref. [3]. In panels (a-d), we plot absolute estimates based on the Lee et al. [3] transform; in panels (e-f), we plot relative estimates (minus simulated values). We repeated each simulation 50 times; means and 95% confidence intervals are plotted.
While the absolute h g 2 estimates exhibit downward bias (a-d), as expected due to ascertainment [17, 27], the relative h g 2 contributions of different variance components were correctly estimated, as were their standard errors (e-f).   Figure 19. REML heritability partitioning under case-control ascertainment: Simulations with fewer causal SNPs. The simulation setup and panels are as in Supplementary Fig. 16 (see caption for details) but with 2,000 instead of 20,000 causal SNPs, 5% of liability explained by ancestry, and 10 principal component covariates. We note that the sum of h g,MAF 2 over 6 MAF bins is slightly lower (by a factor of ≈0.95) than h g 2 estimated using all SNPs in one variance component. While this slight difference is not statistically significant (given the standard errors of our REML analyses), it is consistent with the observation of Speed et al. [28] that mismatch between the true MAF-dependence of allele effect sizes (parameterized by α; see Online Methods) and the modeled MAF-dependence (α = −1) leads to slight bias in estimated SNP-heritability. Indeed, the sign and magnitude of the difference are roughly as predicted by the simulations of ref.
[28] (interpolated for our estimated α = −0.28 for schizophrenia). We also note that our MAF-partitioned analyses only approximately allow the variance component model to adjust for MAF-dependent genetic architectures (as all SNPs within a given MAF bin are normalized according to the standard assumption α = −1); however, given the small MAF range within each bin, this procedure effectively estimates h g, MAF 2 [29].  Figure 1 with additional details.) We benchmarked BOLT-REML and GCTA in three heritability estimation scenarios: (a) partitioning across six MAF bins, (b) partitioning across 22 chromosomes, and (c) bivariate analysis. Runs used subsets of the GERA cohort with fixed SNP count M =597,736 and increasing sample size (N ) and analyzed dyslipidemia as the phenotype in the univariate analyses and hypertension as the second phenotype in the bivariate analysis. Reported run times are medians of five identical runs using one core of a 2.27 GHz Intel Xeon L5640 processor. Reported run times for GCTA break down total CPU time into time required for computing the GRM and time required for performing REML analysis. "NA" indicates scenarios in which GCTA required more memory than the 96GB available. Software versions: BOLT-REML, v2.1; GCTA, v1.24. For bivariate analysis (with all SNPs in a single variance component), the default GCTA GRM computation used more memory than GCTA REML analysis at the sample sizes we tested. Given that the memory footprint of GCTA GRM computation can easily be circumvented (e.g., by combining GCTA-computed per-chromosome GRMs or by using PLINK2 [60]), we plotted only the memory use of GCTA REML analysis in Fig. 1b; we report the memory use of GCTA GRM computation in these tables for completeness. We also note that while PLINK2 allows efficient parallelization of the GRM computation, the total CPU time required is similar to that of GCTA GRM computation (data not shown). Schizophrenia cases and controls are from the PGC2 data set [12]; cases and controls for the other 22 diseases are from the GERA data set. For all diseases, the reported h g-cc 2 is the raw REML heritability parameter obtained by running BOLT-REML on the full data set. For schizophrenia, the h g 2 estimate was computed using PCGC regression (assuming a population risk of 1%) to avoid bias induced by over-ascertainment of cases [17, 27] (Supplementary Table 6). For GERA diseases, h g 2 estimates and were computed from h g-cc 2 using the liability transform of ref. [3] (assuming random sample ascertainment). We restricted further analyses to the top nine GERA diseases, eliminating 11 diseases with very low h g-cc 2 ≤0.025 as well as two heterogeneous traits (all types of cancer and all psychiatric diseases).

Supplementary
Supplementary Table 4. Effect of sample size on REML heritability parameter estimates (h g-cc 2 ) in ascertained and non-ascertained traits. We report BOLT-REML h g-cc 2 estimates for schizophrenia and the three GERA traits with highest h g-cc 2 on subsamples of the the PGC2 and GERA data sets. The first row contains estimates using all samples (corresponding to h g-cc 2 values reported in Supplementary Table 3). Subsequent rows show the results of randomly splitting the data into F = 2, . . . , 10 folds and computing h g-cc 2 for each fold. For each value of F , we repeated this experiment F times (using different random splits); values reported are means over F 2 estimates (s.e.m. in parentheses). The analytic REML s.e. estimates for the full-sample (1x) estimates are 0.008, 0.009, 0.009, 0.009 for the four traits; we note that there are no significant differences between the full-sample vs. subsampled h g-cc 2 for the GERA traits. The liability-scale estimates of h g 2 for 1x and 10x subsamples of the schizophrenia data are 0.232 and 0.261, respectively (assuming 1% population risk). We survey estimates of schizophrenia h g 2 from recent studies. All estimates reported above assume population risk of 1% and do not apply an adjustment for incomplete tagging. We observe that our estimate, h g 2 =0.274 (s.e.=0.007) in the top line, is significantly higher than the estimate h g 2 =0.23 (s.e.=0.01) of Lee et al. [5] using the full PGC1 data set. Our analyses suggest that two factors together fully explain this discrepancy: (1) downward bias of large-sample REML h g 2 estimates due to case over-ascertainment [17,27] (Supplementary Table 4), which we accounted for in this study; and (2) LD bias of h g 2 estimates due to LD-dependent genetic architecture [28]. While previous work has shown that REML analysis of LD-pruned data (as in this study) generally produces larger h g 2 estimates than direct REML analysis of all imputed markers [30] (as in ref. [5,11]), the question of which procedure is unbiased (or less biased) remains open. We are aware of preliminary studies suggesting REML h g 2 estimates based on unpruned marker sets incur downward bias due to enrichment of heritability in lower-LD SNPs (B.K.B., H.K.F., V. Anttila; personal communication), but further investigation is needed. We also note that our subsampling experiments (Supplementary Table 4) and simulations (Supplementary Table 6) indicate that at the sample size analyzed by Golan et al.

Supplementary
[17], we would expect REML estimates of schizophrenia h g 2 to be downward biased by <0.01, so our interpretation of the large but non-significant difference between REML and PCGC estimates in ref.
[17] is that it is mostly due to noise. In table (a), we report data from BOLT-REML runs on simulated case-control data sets with case over-ascertainment. We simulated genotypes by generating individuals as mosaics of up to 100 random "ancestors" from the PGC2 controls, resampling ancestors every 500 SNPs. We simulated case-control phenotypes using a liability threshold model in which we first generated continuous phenotypes with h g 2 =0.166-0.31 explained by 20,000 markers and then defined cases as individuals with phenotypes exceeding thresholds corresponding to the desired prevalences K=0.001-0.02. In each simulation, we ascertained 22,177 cases and 27,629 controls, matching the PGC2 schizophrenia data we analyzed. We ran BOLT-REML on the full simulated data (N =50K), obtaining h g-cc 2 , which we transformed to h g 2 (ref. [3]) using the simulated population prevalence; we also averaged 100 h g-cc 2 estimates on 10x subsamples of the data (Supplementary Table 4) and transformed this estimate to h g 2 , and we compared the estimate based on N =5K subsamples to the N =50K estimate. Reported values are the mean and s.d. over 10 independent simulations. In table (b), we report results on the real PGC2 schizophrenia data assuming population risk of 1% and 0.4%.

Supplementary
From (a), we observe that the N =50K h g 2 estimates are substantially (11-14%) downward biased, as expected [17,27], whereas subsampling to N =5K yields h g 2 estimates with at most faint downward bias. We also see that the downward bias h g-cc 2 (50K) / h g-cc 2 (5K) worsens with decreasing disease prevalence (i.e., increasing severity of ascertainment), suggesting that prevalence can be inferred from the extent of the bias. However, the observed h g-cc 2 (50K) / h g-cc 2 (5K) of 0.889 in (b) is within two standard deviations of the mean h g-cc 2 (50K) / h g-cc 2 (5K) in all our simulations (K=0.001-0.02), indicating that a much larger sample size would be required to make precise inference. We filtered out regions with length <80% of the nominal region length (e.g., due to region truncation at centromeres and telomeres); these truncated regions accounted for 3-6% of the ≈2,600 1Mb regions and ≈5,300 0.5Mb regions. The region boundaries and numbers of regions in the PGC2 and GERA analyses were slightly different because PGC2 data was in build 37 coordinates while GERA data was in build 36 coordinates. For each region, we computed GC content based on the human reference, genic content as the complement of the "intergenic" annotation of ref. [6], replication timing as the region-wide average of data from ref.

Supplementary
[33], recombination rate as the ratio of genetic map distance (according to HapMap phase II [34]) to physical distance spanning each region, and background selection as the region-wide average of data from ref.
[35]. We filtered out regions with length <80% of the nominal region length (e.g., due to region truncation at centromeres and telomeres); these truncated regions accounted for 3-6% of the ≈2,600 1Mb regions (a,b) and ≈5,300 0.5Mb regions (c,d). (The region boundaries and numbers of regions in the PGC2 and GERA analyses were slightly different because PGC2 data was in build 37 coordinates while GERA data was in build 36 coordinates.) We then computed correlation coefficients between these annotations and per-region SNP-heritability estimates after regressing out SNP count per region (a,c) or both SNP count and GC content (b,d). We note that the statistical significance threshold after correcting for multiple hypotheses tested is p<0.001. We report correlation coefficients between pairs of annotations across regions after filtering and regressing out SNP count as described in Supplementary Table 10.

Supplementary
Supplementary We compare partitions of schizophrenia h g 2 across GC quintiles and MAF bins computed using the following techniques: • PCGC regression. This method is not biased by case-control ascertainment, so we report both absolute and relative (% of h g 2 ) estimates. We ran PCGC regression with 10 PCs per variance component (as recommended by ref. [17]) as well as sex and 29 study indicators.
• BOLT-REML. We only report relative (% of h g 2 ) estimates, which remain valid despite the downward bias of absolute REML h g 2 estimates under case-control ascertainment ( Supplementary Figures 16-19). REML rows correspond to: -BOLT-REML using the set of covariates (sex, 29 study indicators, and 10 genome-wide PCs) and QC procedures used in our main analyses (Online Methods).
-BOLT-REML using 10 additional PC covariates per variance component.
-BOLT-REML using our usual set of covariates but run on a subset of 321,567 markers remaining after extremely stringent QC (missing rate <0.01).
All analyses produced consistent results.
Supplementary Table 15. Detailed genetic correlations and total correlations of GERA disease liabilities.
(a) BMI not included as a covariate (This table provides numeric data plotted in Figure 5 with additional details.) (a) Correlations from bivariate analyses using only age, sex, 10 principal components, and Affymetrix kit type as covariates. (b) Correlations from bivariate analyses including BMI as an additional covariate. Genetic correlations are above the diagonals; total liability correlations are below the diagonals. We compare partitions of schizophrenia h g 2 across GC quintiles and MAF bins computed by running BOLT-REML on the full PGC2 data set we analyzed (N =50K) as well as one random N =25K subsample and one random N =10K subsample. Standard errors decrease with increasing sample size, as expected. We observed that all methods computed accurate estimates with negligible bias. We expected from theory (Supplementary Note) that BOLT-REML estimates using S Monte Carlo trials would have variance 1+1/S times larger than exact REML estimates (e.g., computed using GCTA), and our simulations confirmed that indeed, to within the statistical resolution available from 50,000 replicates, the RMSE of BOLT-REML was greater than that of GCTA by a factor of 1 + 1/S. (We also observed that in this particular simulation setup, standard errors reported by REML-computed using analytic approximations based on Fisher information-were slightly conservative compared to observed RMSE.)

Supplementary
Supplementary where G chr1 and G chr2 denote genetic relationship matrices computed using SNPs on chromosomes 1 and 2, respectively, and I N denotes the N × N identity matrix (corresponding to environmental effects). We then estimated these variance parameters by running bivariate BOLT-REML (with 15 Monte Carlo trials) using one variance component per chromosome.  We simulated 100 pairs of correlated quantitative traits exactly as in Supplementary Table 18 and then converted each trait to a case-control trait according to a liability threshold model. Specifically, for each pair, we used a case cutoff of z > 0 for trait 1 (corresponding to 50% population prevalence) and a case cutoff of z > 1 for trait 2 (corresponding to 16% population prevalence). We then estimated the variance parameters of the underlying liabilities (i.e., the original quantitative traits) by running bivariate BOLT-REML (with 15 Monte Carlo trials) using one variance component per chromosome and transforming parameter estimates to the liability scale. Explicitly, we used the formula of ref. [3] to transform h g,chr 2 ; we did not transform r g (ref. [7]); and we transformed total phenotypic correlation to total liability correlation as described in Online Methods, after which we computed r e by subtracting the two genetic covariances from the total liability covariance.