CORE GREML: Estimating covariance between random effects in linear mixed models for genomic analyses of complex traits

Linear mixed models (LMMs) using genome-based restricted maximum likelihood (GREML) are a key variance partitioning tool, where effects of multiple sources, such as different functional genomic regions, on phenotypes are treated as random. Classic LMMs assume independence between random effects, which can cause biased estimation of variance components. Here, we relax this independence assumption by introducing a generalised GREML, called CORE GREML, that can explicitly estimate the covariance between random effects. Using extensive simulations, we show that CORE GREML outperforms the conventional GREML, providing unbiased estimates of variance and covariance components. Using data from the UK biobank, we demonstrate that CORE GREML is useful for genomic partitioning analyses and for genome-transcriptome partitioning of phenotypic variance. For example, we found that the transcriptome, imputed using genotype data, explained a significant proportion of phenotypic variance for height (0.15, se = 5.4e-3, p-value = 1.5e-283), and that these transcriptomic effects on phenotypes correlated with effects of the genome (r = 0.35, se = 4.6e-2, p-value = 1.2e-14). We conclude that the covariance between random effects is a key parameter that needs to be estimated, especially when partitioning phenotypic variance by omic layer.


Introduction 20
Genome-wide association studies (GWAS) have been incredibly successful in 21 identifying genetic variants associated with complex traits. However, the proportion 22 of phenotypic variance explained by genome-wide significant single nucleotide 23 polymorphisms (SNPs) is far lower than the narrow-sense heritability estimate 1 . This 24 is largely because GWASs typically examine SNP-trait associations one at a time, 25 generating a large number of tests across the genome for which the stringent 26 Bonferroni correction is applied. To overcome this problem, a whole-genome 27 approach that jointly considers all available SNPs has been introduced 2 , allowing 28 estimation of the proportion of phenotypic variance explained by genome-wide SNPs, 29 i.e. SNP-based heritability. Central to this approach is the use of linear mixed 30 models 3 (LMMs)-extensions of random-effects models or variance-component 31 models 4 -that treat SNP effects as random. 32 33 Using genome-based restricted maximum likelihood (GREML) for parameter 34 estimation, LMMs are a key tool not only for SNP-based heritability of complex traits 35 but also for variance partitioning in general. For example, when heritability is 36 partitioned by functional annotation of SNPs using a LMM with multiple random 37 effects, GREML estimation has provided important insights into the latent genetic 38 architecture of complex traits 5 . As multi-omics data become increasingly available 6 , 39 variance partitioning using LMMs will become indispensable to uncover the relative 40 contributions of multiple 'omes' to phenotypic variation. Alongside GREML, linkage 41 disequilibrium score regression (LDSC) provides an alternative way for SNP-based 42 heritability estimation and variance partitioning, using only GWAS summary statistics 43 without the need to access individual genotypes 7,8 . LDSC or stratified LDSC also 44 treats SNP effects as random 8 ; and for the same set of individual-level genotype 45 data, this approach generates similar estimates as GREML [9][10][11] . 46 47 Following classic LMMs 3,4 , GREML assumes independence between random effects 48 when estimating variance components. However, it is questionable if this assumption 49 is always valid especially for genomic analyses of complex traits. For example, gene 50 regulatory networks shared between functional categories may generate non-51 negligible correlations between effects of these categories on phenotypes 12 . In the 52 context of phenotypic variance partitioning by multi-omics layers, effects of genetic 53 variants and their expression levels on phenotypes are likely correlated [13][14][15] , as 54 exemplified by overlaps between GWAS loci and expression quantitative trait loci 55 (eQTL; e.g. 16,17 ). Given these justifiable covariance terms in genomic analyses of 56 complex traits, the naïve assumption of independence between random effects held 57 by GREML can lead to biased partition of phenotypic variance and false inferences 58 on the underlying architecture of complex traits. 59 60 Here, we introduce an alternative GREML, referred to as CORE GREML (CORE for 61 COvariance between Random Effects), that fits the Cholesky decomposition of 62 variance, the kernel matrix for genetic variance estimation was constructed using 105 genotypes of 1,133,273 genome-wide SNPs, and the kernel matrix for the estimation 106 of phenotypic variance explained by the transcriptome was based on imputed 107 expression levels of 227,664 genes from 43 tissues 18 (see Supplementary Table 2 Incorrect assumptions in the estimation model about the genetic architecture of the 142 trait can also bias variance-components estimation in the context of GREML (e.g. 19 ). 143 Therefore, we tested the extent to which CORE GREML estimation is sensitive to a 144 wrong assumption of the genetic architecture in the estimation model. This was 145 achieved by simulating phenotypes under different genetic architectures (see 146 Methods) and comparing CORE GREML estimation from fitting an estimation model 147 that has the correct assumption about the genetic architecture, referred to as the 148 'true model', with that from fitting a 'wrong model' that has an incorrect assumption 149 about the genetic architecture. 150 151 We found that misspecification of genetic architecture in the estimation model in 152 general biased CORE GREML estimation of variance components but not the 153 covariance term for genome-transcriptome partitioning of phenotypic variance (see 154 Supplementary Figure 4). Nonetheless, under any given genetic architecture, 155 misspecification can be feasibly diagnosed by comparing the likelihood of estimation 156 models that assume a wide range of possible genetic architectures. As shown in 157 Supplementary Figure 5, differences in the likelihood of estimation models are highly 158 indicative of deviations from the true underlying genetic architecture. In fact, a grid 159 search approach has been practiced choosing the model closest to the true 160 underlying genetic architecture in the GREML context 19 .

162
In light of the above results, to reduce the chance of misspecification of genetic 163 architecture for real data analyses, we fitted two estimation models, the Genome-164 wide Complex Trait Analysis (GCTA) model 2 Table 4; although the GCTA model seems more conservative than the LDAK model 175 for detecting covariance terms). This is consistent with the previous observation that 176 heritability estimates based on high-quality common SNPs are robust to variation in 177 assumed genetic architecture (more specifically, parameter α values of the LDAK 178 model; see Figure 6 & Supplementary Figure 4 in Speed et al. 19  we fitted two models using GREML, a 'G model' that breaks phenotypic effects into 194 effects of the genome and residuals, i.e., y=g+ε, and a 'G-T model' that decomposes 195 phenotypic effects into effects of the genome and the imputed transcriptome and 196 residuals, i.e., y=g+t+ε. We declared the presence of phenotypic effects of the 197 imputed transcriptome when the G-T model had a better fit than the G model using 198 the likelihood ratio test with one degree of freedom. We found significant phenotypic 199 effects of the imputed transcriptome for all traits except fluid intelligence (Figure 1).

200
Interestingly, although the G-T model had a much better fit than the G model for the 201 nine traits, the two models explained a similar amount of phenotypic variance ( Figure  202 1), which was verified by additional simulations (see Supplementary Note). This 203 suggests that the partition of phenotypic variance represented by the G-T model is 204 closer to the true underlying model than that represented by the G model.

206
To validate the transcriptomic effects on phenotypes revealed by the G-T model, we 207 performed a 5-fold cross-validation, in which the phenotypic prediction accuracy of 208 the G-T model was compared against that of the G model. For each trait, we 209 randomly split the sample into a training set (~80%) and a validation set (~20%) and 210 iterated this process five times in a manner such that validation sets did not overlap 211 across iterations. To derive the prediction accuracy for the two models, we computed 212 the Pearson's correlation coefficient between the observed and predicted 213 phenotypes of each trait in each iteration and averaged correlation estimates across 214 five iterations. Figure 2 a. shows that the gain in the phenotypic prediction accuracy 215 by the G-T model relative to that by the G model grew as the estimated 216 transcriptomic contribution to phenotypic variance increased (p = 1.86e-06), 217 suggesting that the transcriptomic effects on phenotypes of the selected traits are 218 genuine. Taken together, our results thus far indicate that the variance components 219 of interest differ from zero for all traits with fluid intelligence being the exception for 220 the genome-transcriptome partitioning of phenotypic variance. Importantly, these 221 results established the basis for our subsequent covariance estimation.

223
By comparing model fit by GREML and CORE GREML, we detected significant 224 covariance between phenotypic effects of the regulatory regions and DHS for height 225 and sitting height ( Figure 3). Of note, the genomic partitioning model included three 226 covariance terms, but two of these terms were not significant for height and sitting 227 height, based on the Wald test with one degree of freedom. We therefore reduced 228 the model for these two traits by dropping non-significant covariance terms, noting 229 that the fit of the reduced model did not differ from the full model (p = 0.85 & 0.37 for 230 height and sitting height, respectively). We subsequently used estimates from the 231 reduced model for these two traits. In genome-transcriptome analyses, we found 232 significant covariance between phenotypic effects of the genome and those of the 233 imputed transcriptome for height, sitting height, heel bone mineral density, and 234 diastolic blood pressure ( Figure 4). We standardized all estimated covariance terms 235 using respective variance estimates to derive correlation estimates (Figures 3 & 4 far 236 right), noting that all significant estimates were positive and small to moderate in size, 237 ranging from 0.14 (standing height in Figure 3) to 0.58 (heel bone mineral density in 238 Figure 4). In a subsequent sensitivity analysis, all significant covariance terms 239 emerging from the genomic partitioning and genome-transcriptome analyses 240 remained after applying a rank-based inverse normal transformation to phenotypic 241 observations (see p-values in Supplementary Figures 6 & 7). Thus, the estimated 242 covariance terms were robust against the violation of the normality assumption held 243 by GREML and CORE GREML.

245
To show that covariance terms estimated by CORE GREML are genuine biological 246 parameters, we validated the covariance between phenotypic effects of the genome 247 and the imputed transcriptome using the same 5-fold cross-validation procedure as 248 before (i.e., for the validation of transcriptomic effects on phenotypes). In this 249 instance, the phenotypic prediction accuracy of CORE GREML was compared 250 against that of GREML. We chose genome-transcriptome analyses for validation 251 since significant covariance emerged from four traits in contrast to two traits in 252 genomic partitioning analyses. Figure 2 b. shows that the gain in the phenotypic 253 prediction accuracy of CORE GREML relative to that of GREML grew as the 254 magnitude of covariance estimates increased (p = 1.61e-05).

256
To show the impact of neglecting significant covariance terms, we compared Based on these observations, the larger variance estimates by GREML for 272 correlated random effects relative to those by CORE GREML are most likely due to 273 bias from neglecting the correlations between these random effects.

275
We also considered the impact of neglecting covariance between random effects on been standardized prior to analyses). The functional significance of SNPs from a 287 given genomic region can be tested by assessing if the proportion of genetic 288 variance attributable to the region is substantially higher than the the proportion of 289 SNPs from the region 5 . Notably, GREML estimates of these functions were larger 290 than CORE GREML estimates whenever there was a significant covariance term 291 ( Supplementary Figures 8-10); but the two methods agreed with each other 292 otherwise. For example, the relative phenotypic contributions by the genome and the 293 imputed transcriptome inferred from GREML were larger than those inferred from 294 CORE GREML for heel bone mineral density, diastolic blood pressure, sitting height 295 and standing height ( Figure 8). Similarly, the functional significance of SNPs from 296 regulatory regions inferred by GREML was larger than that by CORE GREML for 297 sitting height and standing height ( Figure 6b). Taken together, our results indicate 298 that GREML can lead to incorrect inferences on the underlying architecture of 299 complex traits unless correlations between random effects are properly modelled. 300 301 It is important to note that kernel matrices used for variance-components estimation 302 in a LMM can be similar, for instance, due to linkage disequilibrium in a genomic 303 partitioning analysis. In fact, the correlations between off-diagonal entries of kernel 304 matrices used in our analyses are moderate to high (0.35-0.98; see Supplementary 305 Table 6). This similarity might give rise to the covariance between random effects. 306 However, this possibility is unlikely for at least two reasons. First, if covariance is 307 driven by the similarity between kernel matrices, then we would expect that in the 308 null setting of our simulations, type I error rate is inflated, given the high similarities 309 between kernel matrices in our analyses. Contrary to this, we found that type I error 310 rate is controlled. Second, the kernel matrix constructed using genotypes from DHS 311 is more similar to the kernel matrix constructed using genotypes from other regions 312 than to the one for regulatory regions; but significant covariance was only detected 313 between effects of regulatory regions and those of DHS on standing height and 314 sitting height. Therefore, covariance between random effects is unlikely driven by the 315 similarity between kernel matrices for variance-components estimation. 316 317 318

320
When applying the classic LMMs for standard heritability estimation, where 321 phenotypic variance is only partitioned into genetic and residual variances, the model 322 assumption of negligible covariance between random effects (i.e., genetic and 323 residual effects) may be met in many cases. However, when phenotypic variance is 324 further partitioned, for example, by functional genomic region or omic layer, using a 325 model with multiple random effects, the covariance terms between these random 326 effects can be substantial, as what we demonstrated using the genomic partitioning 327 analyses and the genome-transcriptome partitioning analyses for complex traits. 328 Unless these non-negligible covariance terms are properly accounted for, variance-329 components estimation would be biased, resulting in misleading inferences on the 330 latent omic architecture of complex traits, as shown by simulation results. Therefore, 331 we recommend that covariance terms between random effects need to be carefully 332 checked and properly modelled for genomic analyses of complex traits. 333 334 CORE GREML can serve as a useful tool for detecting and estimating covariance 335 terms between random effects, as demonstrated using analyses of both simulated 336 and real data. Prior to the proposal of CORE GREML, there have been several 337 attempts to relax the assumption of independence between random effects 20,21 , but 338 they are specific to experimental studies and are not readily applicable to genome-339 wide analyses for human complex traits. To our knowledge, CORE GREML is the 340 first of its kind in variance partitioning analyses, which correctly models the 341 covariance between random effects. 342 343 We demonstrated the use of CORE GREML in genomic partitioning analyses and 344 genome-transcriptome partitioning of phenotypic variance, and found that significant 345 covariance terms mostly emerged from the latter (for 4 traits out of 10 Finally, we showed, using simulations, that misspecification of LD and MAF 400 dependent genetic architecture can cause substantial bias in variance-components 401 estimation by CORE GREML, although covariance estimation seems robust. 402 However, the likelihood of the true estimation model in general is much greater than 403 a wrong model, suggesting that likelihood-based comparisons of models that 404 assume different genetic architectures is a useful way to reduce the chance of mis-405 specifying genetic architecture. For demonstration, we fitted the GCTA model and 406 the LDAK model with the recommended default setting 19 and found the GCTA model 407 in general had a better fit than the LDAK model for our traits. In the absence of the 408 knowledge of the true genetic architecture, however, it is recommended to more 409 systematically vary parameter settings of the LDAK model (as in 19 ) and choose the 410 best fitting model via likelihood comparison before applying CORE GREML.

412
In this study, we introduce a generalised GREML, referred to as CORE GREML, that 413 relaxes the assumption of independence between phenotypic effects of random 414 terms held by classic mixed-effects models for variance component analyses. Using 415 both simulations and real data, we showed that in the presence of non-negligible 416 covariance terms, CORE GREML improved genomic partitioning and multi-omics 417 partitioning analyses by the conventional GREML. We conclude that the covariance 418 between random effects for analysis of complex traits is a key parameter for 419 estimation, and hence, recommend that covariance terms should be carefully 420 checked and properly modelled. 421

Methods 422
Generalizing variance-covariance matrix of phenotypic observations

424
A LMM can be written as where y is a vector of trait phenotypes, b is a vector of fixed effects, g is a vector of 429 additive genetic effects and e is a vector of residual effects. X and Z are incidence where , is the covariance between g and e, i.e. , = ( , ) • √ 2 • 2 .

448
When considering multiple random effects in the LMM (e.g. genomic partitioning 449 approach), the model can be written as where gi is the random genetic effects of the ith pre-defined functional category, e.g. 454 regulatory regions.

456
Such mixed-effects models with multiple random terms typically assume no 457 correlation between gi and gj. However, this assumption can be violated if effects of 458 two categories on phenotypes are associated, for example, through the same gene 459 pathway. We relax this assumption and write the variance-covariance matrix of all 460 observations, var(y), as 461 462 where Ai is the genetic relationship kernel matrix constructed using SNPs from the 465 ith functional category, √ is the Cholesky decomposition of Ai and is the 466 genetic covariance between gi and gj. It is noted that the correlation term between 467 genetic effects (gi) and residuals (e) is not included and hence assumed to be zero 468 in Eq. (4), which is usually valid, although it is possible to parameterise this term.

470
The log likelihood of the proposed model, which can be generally applied to Eq. (1) 471 and (4), is 472 473 where ln is the natural log and | | the determinant of the associated matrices. The 476 projection matrix is defined as = − − − ( ′ − ) − ′ − . By maximising the 477 log likelihood, the direct average information algorithm 29,30 can be used to obtain 478 CORE GREML estimates of parameters including the covariance terms between 479 random effects.

481
This CORE GREML approach can be easily extended to phenotypic variance 482 partitioning using multi-omics data, e.g. genome-transcriptome analyses (see 483 Genome-Transcriptome Partitioning Model section below). where the phenotypic variance is 2 = 2 + 2 in the absence of cor(g, e). When 492 there is non- negligible cor(g, e), the phenotypic variance should be written as σ y 2 = 493 σ g 2 + σ e 2 + 2 • cov( , ). 494 495 For Eq. (4), a general expression of heritability for the ith genetic component is 496 497 where σ y 2 = ∑ σ g i Using the Delta method (see 31  where var(σ g i 2 ), var(σ y 2 ) and cov(σ g i 2 , σ y 2 ) can be obtained from the average 507 information matrix 32 of CORE GREML.

509
Correlation between two random effects 510 511 The correlation between two random (genetic) effects ( ) can be defined as the 512 genetic covariance scaled by the square root of the product of the genetic variances 513 of the two random effects, that is 514 515 Using the Delta method (see 31 ), the sampling variance of genetic correlation can be 518 obtained as 519 520 var (r g i g j ) = r g i g j 2 ( var(σ g i 2 ) 4•σ g i 4 + var(σ g j 2 ) 4•σ g j 4 + var(σ g i g j ) σ g i g j 2 + 2•cov(σ g i 2 ,σ g j 2 ) 4•σ g i 2 •σ g j 2 − 2•cov(σ g i 2 ,σ g i g j ) 2•σ g i 2 •σ g i g j − 521 2•cov(σ g j 2 ,σ g i g j ) 2•σ g j 2 •σ g i g j ) 522 where the variance and covariance terms used are from the information matrix of 523 CORE GREML.

525
Real data analysis

527
Genotype data 528 The UK biobank (project approval number 14575) contains health-related data from 529 ~ 500,000 participants aged between 40 and 69, who were recruited throughout the 530 UK between 2006 and 2010 33 . Prior to data analysis, we applied stringent quality 531 control to exclude unreliable genotypic data. We filtered SNPs with an INFO score < 532 0.6, a MAF < 0.01, a Hardy-Weinberg equilibrium p-value <1e-4, or a call rate < 0.95. 533 We then selected HapMap3 SNPs, which are known to yield reliable and robust 534 estimates of SNP-based heritability 25,34,35 , for downstream analyses. We filtered 535 individuals who had a genotype-missing rate > 0.05, were non-white British ancestry, 536 or had the first or second ancestry principal components outside six standard 537 deviations of the population mean. We also applied relatedness-cut-off quality 538 control to exclude one of any pair of individuals with a genomic relationship > 0.025. 539 From the remaining individuals, we selected those who were included in both the first 540 and second release of UK biobank genotype data. We calculated the discordance 541 rate of imputed genotypes between the two versions and excluded individuals with a 542 discordance rate > 0.05. Eventually, genotypes of 1,131,002 SNPs from 91,472 543 individuals remained for data analysis. 544

545
To preclude negligible heritability as a possibility for negative findings (i.e., no 546 covariance between random effects), we deliberately chose ten UK biobank traits 547 available to us with the largest heritability estimates by an independent open source 548 (https://nealelab.github.io/UKBB_ldsc/), which included standing height, sitting height, 549 body mass index (BMI), heel bone mineral density, fluid intelligence, weight, waist 550 circumference, hip circumference, diastolic blood pressure and years of education 36 .

551
Heritability estimates for all selected traits were at least 20 times greater than their 552 standard errors to ensure that they were significantly different from zero. We further 553 verified SNP-based heritability of these traits using GREML and estimates are shown 554 in Supplementary Figure 1. 555 556 Prior to model fitting, phenotypic data were prepared in three sequential steps: 1) 557 adjustment for age, sex, birth year, social economic status (by Townsend 558 Deprivation Index), population structure (by the first ten principal components of the 559 estimated genomic relationship matrix), assessment centre, and genotype batch 560 using linear regression; 2) standardization; and 3) removal of data points outside +/-561 3 standard deviations from the mean. The distributions of phenotypes of the ten traits 562 are shown in Supplementary Figure 11. We noted mild to strong deviations from 563 normality for traits such as BMI and years of education. This motivated a subsequent 564 sensitivity analysis to test the robustness of our findings against the violation of the 565 normality assumption held by GREML and CORE GREML. Specifically, we applied a 566 rank-based inverse normal transformation to phenotypes of all traits and repeated 567 our analyses on the transformed phenotypes.

569
Functional annotation of the genome 570 The genome was annotated using three pre-defined functional categories 571 ( genotype data were quality-controlled (see above for details). 577

578
Using PrediXcan 18 , we imputed expression levels of 2,028 to 9,630 genes for 43 579 non-sex-specific tissues (Supplementary  588 The phenotypic variance of each selected trait was partitioned using two separate 589 random-effects models (see below for model description). The 'genome-590 transcriptome' model partitions phenotypic variance into variation from the genome, 591 the transcriptome and unknown sources (i.e., residual variance), whereas the 592 'genomic partitioning' model assumes phenotypic variation comes from the genome 593 and residuals and further partitions genetic variance by functional category of the 594 genome. Three functional categories were under consideration, namely, regulatory 595 regions (encompassing coding regions, untranslated regions and promotors), DNase 596 I hypersensitivity sites (DHS) and all other regions. 597 598 Each partitioning model was fitted using the conventional method, i.e., GREML, and 599 the proposed alternative, i.e., CORE GREML. Essentially, GREML sets all 600 covariance terms between random effects to zero, whereas CORE GREML treats 601 these terms as free parameters for estimation. To detect significant covariance terms, 602 we performed likelihood ratio tests to determine whether the model fit by CORE 603 GREML was better than that by GREML. Where y is a n x 1 vector of phenotype data and their deviations from the grand 608 mean, μ, are decomposed into three vectors consisting of phenotypic effects of the 609 genome, g ~ N(0, nxn σ g 2 ), those of the transcriptome, t ~ N (0, nxn σ t 2 ), and 610 residuals, ε ~ N (0, nxn σ ε 2 ). The terms, σ g 2 , σ t 2 and σ ε 2 denote phenotypic variances 611 attributable to the genome, the transcriptome and residuals, respectively. nxn and 612 nxn are relationship kernel matrices and nxn is an identity matrix. nxn is derived by  To validate CORE GREML, we simulated 500 replicates of phenotypic data using the 651 two variance partitioning models shown above under each of three parameter 652 settings: zero (i.e., null setting), positive and negative covariance between 653 phenotypic effects of random terms in the variance partitioning model. Simulations 654 were based on quality-controlled genotype data and imputed transcriptome data 655 from a random sample of 10,000 UK biobank individuals. The genotype data 656 contained a total of 1,131,002 SNPs (see genotype data above), and the imputed 657 transcriptome contained imputed expressions of 227,664 genes collapsed cross 43 658 non-sex specific tissues (see gene expression imputation above).

660
For genome-transcriptome analysis, phenotypes were simulated using eq. 6 661 according to the following variance-covariance structure of random effects: 662 663 where the value of σ gt varied across parameter settings (Supplementary Table 1).

666
For genomic partitioning analysis, phenotypes were simulated using eq. 7 according 667 to the following variance-covariance structure of random effects: 668 669 . Imputed transcriptome contributes to phenotypic variance. Shown are estimated variance components as a proportion of total phenotypic variance from a linear mixed model that includes a random term for the phenotypic effects of the imputed transcriptome and another model that does not, denoted as y = g+t+ε and y = g+ε, respectively. P-values are from likelihood ratio tests (df = 1) that compare the two models for detecting phenotypic effects of the imputed transcriptome. g = phenotypic effects of the genome; t = phenotypic effects of the imputed transcriptome; ε = residuals; σ t 2 = phenotypic variance explained by the imputed transcriptome; σ g 2 = phenotypic variance explained by the genome; σ y 2 = total phenotypic variance. The imputed transcriptome consists of expression levels of 227,664 genes from 43 non-sex specific tissues. Figure 2. Imputed transcriptome and covariance between phenotypic effects of the genome and of the imputed transcriptome improve phenotypic prediction accuracy. Panel a. the more the imputed transcriptome contributes to phenotypic variance, the greater the gain in phenotypic prediction accuracy. σ t 2 denotes the phenotypic variance explained by the imputed transcriptome. Prediction accuracy was computed using the Pearson's correlation coefficient between the observed and the predicted for models with and without a random term for phenotypic effects of the imputed transcriptome, denoted as y = g+t+ε and y = g+ε, respectively. g = phenotypic effects of the genome; t = phenotypic effects of the imputed transcriptome; ε = residuals. Prediction accuracy improvement (i.e., y-axis) was derived by subtracting the prediction accuracy of the model y = g+ε from that of the model y = g+t+ε. The least squares line with 95% confidence band is based on a linear model that regressed prediction accuracy improvement on phenotypic variance explained by the imputed transcriptome. The pvalue is for the t-test statistic (df=8) under the null hypothesis that the slope of the regression line is zero. Panel b. the large the estimated covariance between phenotypic effects of the genome and those of the imputed transcriptome (i.e., σ g t ), the greater the gain in phenotypic prediction accuracy. Prediction accuracy was computed using the Pearson's correlation coefficient between the observed and the predicted for two models of the form y = g+t+ε, but one assuming σ g t = 0, i.e., GREML, and the other setting σ g t as a free parameter for estimation, i.e., CORE GREML. Prediction accuracy improvement was derived by subtracting the prediction accuracy of GREML from that of CORE GREML. The least squares line with 95% confidence band is based on a linear model that regressed prediction accuracy improvement on σ g t estimates. The p-value is for the t-test statistic (df=8) under the null hypothesis that the slope of the regression line is zero. Figure 3. Estimated covariance (left) and correlations (right) between genetic effects attributable to the regulatory regions and those attributable to the DHS on phenotypes. Error bars are 95% confidence intervals. p1 = p-values from likelihood ratio tests that compare GREML with CORE GREML to detect σ regulatory DHS ; p2 = p-values based on the Wald test statistic under the null hypothesis that r regulatory DHS = 0. Significant σ regulatory DHS and r regulatory DHS are highlighted in orange. Fluid intelligence and years of education are excluded because either the phenotypic effects of the regulatory regions or those of the DHS were not significant. Figure 4. Estimated covariances (left) and correlations (right) between phenotypic effects of the genome and those of the imputed transcriptome. Error bars are 95% confidence intervals. p1 = p-values from likelihood ratio tests that compare GREML with CORE GREML to detect σ g t ; p2 = p-values based on the Wald test statistic under the null hypothesis that r g t = 0. Significant σ gt and r g t are highlighted in orange. Fluid intelligence is excluded because the phenotypic effects of the imputed transcriptome on this trait was not significant after Bonferroni correction.