CORE GREML for estimating covariance between random effects in linear mixed models for complex trait analyses

As a key variance partitioning tool, linear mixed models (LMMs) using genome-based restricted maximum likelihood (GREML) allow both fixed and random effects. Classic LMMs assume independence between random effects, which can be violated, causing bias. Here we introduce a generalized GREML, named CORE GREML, that explicitly estimates the covariance between random effects. Using extensive simulations, we show that CORE GREML outperforms the conventional GREML, providing variance and covariance estimates free from bias due to correlated random effects. Applying CORE GREML to UK Biobank data, we find, for example, that the transcriptome, imputed using genotype data, explains a significant proportion of phenotypic variance for height (0.15, p-value = 1.5e-283), and that these transcriptomic effects correlate with the genomic effects (genome-transcriptome correlation = 0.35, p-value = 1.2e-14). We conclude that the covariance between random effects is a key parameter for estimation, especially when partitioning phenotypic variance by multi-omics layers.


Supplementary
. Sampling distributions of model parameters for the genome-transcriptome partitioning model. Five-hundred replicates of phenotypic data (n = 10,000) were simulated under each of three parameter settings, where the covariance between the random effects of the genome and the imputed transcriptome was zero (cov=0), positive (cov>0) and negative (cov<0). For each replicate, model parameters were estimated by the traditional method, i.e., GREML, and the proposed method, i.e., CORE GREML. and f influence the variance of SNP-specific effects on phenotypes, i.e., var(β), respectively. Panels differ in the combination of α and γ values, which gives rise to different genetic architectures of the simulated trait. Under each genetic architecture, there are three scenarios for the covariance between the random effects of the genome and those of the transcriptome, i.e., cov=0, cov>0 and cov<0. For each scenario, 500 replicates of phenotypes (each with n = 10,000) were simulated; the wrong estimation model always assumes α = -1 and γ = 0, while the true model assumes the values of α and γ identical to those of the simulation model. Regardless of the estimation model, CORE GREML was applied for parameter estimation. Source data are provided as a Source Data file.
Supplementary Figure 5. Estimation model assuming the true genetic architecture of the simulated trait consistently shows a better fit than the model assuming a wrong genetic architecture. Shown in each row is the distribution of differences in log likelihood between the estimation model (true model) that assumes the correct genetic architecture of the simulated trait and a 'wrong model' that mis-specifies the genetic architecture. Irrespective of the estimation model, CORE GREML was used for parameter estimation. The difference in log likelihood is computed by subtracting the log likelihood of the wrong model from that of the true model, such that values above zero indicate that the true model has a better fit than the wrong model. The genetic architecture of the simulated trait is parameterised by linkage disequilibrium structure, w, and minor allele frequency, f, in forms of For each scenario, 500 replicates of phenotypes (each with n = 10,000) were simulated; the wrong estimation model always assumes α = -1 and γ = 0, while the true model assumes the values of α and γ identical to those of the simulation model. Box-plots elements: center line, median; box limits, upper and lower quartiles; whiskers, 1.5 x interquartile range. Source data are provided as a Source Data file.

Supplementary Figure 6. GREML estimates of model parameters by estimation model and p-values for model comparisons.
Five-hundred replicates of phenotypic data (n = 10,000) were simulated under settings with and without the random effects of the imputed transcriptome, denoted as ! ! " > 0 and ! ! " = 0, respectively. For each replicate, two linear mixed-effects models were fitted, one that did not include a term for random effects of the imputed transcriptome, i.e., y = g + ε, and the other did, i.e., y = g + t + ε. All model parameters were estimated using GREML. Panel a. estimated density of model parameters. histogram of p-values from likelihood ratio tests (df=1) that compared the two models (i.e., y = g + ε & y = g + t + ε.) for replicates simulated under ! ! " > 0. P-values are log transformed with the statistical significance threshold, i.e., -log (0.05), being indicated by the orange vertical dash line. Values above the threshold indicate that y = g + t + ε had a better fit than y = g + ε. Source data are provided as a Source Data file. , and ! *!/&% %&#1*23 " denote phenotypic variances explained by the three functional regions; and ! %&#'()!*%+ ,-. denotes the covariance between random effects of the regulatory regions and of the DHS on phenotypes. Error bars are 95% confidence intervals (based on s.e.m). Model parameters were estimated using the traditional method, i.e., GREML, and the proposed method, i.e., CORE GREML. N = sample size; p = p-values from likelihood ratio tests that compared GREML with CORE GREML to detect ! %&#'()!*%+ ,-. ; and p* = p-values from a sensitivity analysis where a rank-based inverse normal transformation was applied to phenotypic data to check the robustness of signals against the violation of the normality assumption held by both GREML and CORE GREML. Highlighted in orange are significant ! %&#'()!*%+ ,-. after a Bonferroni adjustment for multiple comparisons. Residual variance estimates are omitted for simplicity. Fluid intelligence and years of education are excluded because either the random effects of the regulatory regions or those of the DHS were not significant. Source data are provided as a Source Data file.

Supplementary Figure 8. Variance component estimates from genome-transcriptome partitioning analyses. ! #
" and ! ! " denote the phenotypic variances explained by the genome and by the imputed transcriptome, respectively; and ! # ! denotes the covariance between genetic effects and effects of imputed gene expressions on phenotypes. Model parameters were estimated using the traditional method, i.e., GREML, and the proposed method, i.e., CORE GREML. Error bars are 95% confidence intervals (based on s.e.m). N = sample size; p = p-values from likelihood ratio tests that compared GREML with CORE GREML to detect ! # ! ; and p* = p-values from a sensitivity analysis where a rank-based inverse normal transformation was applied to phenotypic data to check the robustness of signals against the violation of the normality assumption held by both GREML and CORE GREML. Highlighted in orange are significant ! # ! after a Bonferroni adjustment for multiple comparisons. Residual variance estimates are omitted for simplicity. Fluid intelligence is excluded because the random effects of the imputed transcriptome on this trait was not significant after Bonferroni correction. Source data are provided as a Source Data file. Histograms of phenotypic data for ten selected traits from the UK Biobank. Traits from top left to bottom right are body mass index, standing height, sitting height, heel bone mineral density, weight, fluid intelligence, years of education, diastolic blood pressure, waist circumference and hip circumference. Data were prepared in three sequential steps: 1) adjustment for age, sex, birth year, social economic status, population structure, assessment centre, and genotype batch; 2) standardization; and 3) removal of data points outside +/-3 standard deviations from the mean. The density function of Normal (0, 1) is superimposed as the reference to highlight deviations from normality. A & T are kernel matrices constructed using genotypes of 1,131,002 SNPs and using imputed expression levels of 227,664 genes collapsed across 43 tissues, respectively. I is an identity matrix. Ai is the genomic relationship matrix constructed using SNPs from functional region i of the genome.   Supplementary Table 4. Likelihood of estimation models when GREML and CORE GREML were applied for the genome-transcriptome partitioning of phenotypic variance. *Fitted GCTA model and LDAK model differ in the assumption about the variance of SNP-specific effects on phenotypes, var(%). In general, it is assumed that for a given SNP i, the variance of its effects on phenotypes is proportional to its linkage disequilibrium score, w, and minor allele frequency, f, i.e.,var(% ! ) ∝ ( ! " [* ! (1 − * ! )] #$% , where parameter γ modulates the effect of w and α modulates the effect of f on var(% ! ). In the GCTA model, α and γ were assumed to be -1 and 0, respectively; in the LDAK model, γ was assumed to 1 and α was set to the recommended default, -0.25 1 . **Δ log likelihood was derived by subtracting the log likelihood of the LDAK model from that of the GCTA model, when CORE GREML was used for parameter estimation for both models.***p-values are from likelihood ratio tests that compared GREML with CORE GREML to detect the covariance between random effects of the genome and those of the transcriptome on phenotypes. P-values are not adjusted for multiple comparisons. Notes. 1. y = phenotypes; t = random effects of the imputed transcriptome; g0 = random effects of 1,074,458 SNPs for genome-transcriptome partitioning analyses presented throughout the main text; g1 = random effects of 1,316,391 SNPs used for the transcriptome imputation; ε = residuals. Log-likelihoodx = log-likelihood of model x, where x = a, b, c, or d; Δ log-likelihoodx-y = log-likelihood of model x -log-likelihood of model y, where x = a or b and y = c or d. 2; and p-values are from likelihood ratio tests that compared models a and b. P-values are not adjusted for multiple comparisons. It is important to distinguish the SNP sets involved for different variance components of the models. The kernel matrix for t is based on gene expressions imputed using 1,316,391 SNPs, which is the same set of SNPs used to construct the kernel matrix for g1. This set of SNPs only has ~ 600K in common with the 1,074,458 SNPs used to construct the kernel matrix for g0. In our primary analysis (model a VS. model b), we show that model a has a better fit than model b for eight traits, even though t and g1 are based on the same set of SNPs. Hence, t is orthogonal to g1 for most traits. In our secondary analysis (model c VS. model d), we show that model d has a better fit than model c for seven traits-once again-despite the fact that t and g1 are based on the same set of SNPs. This confirms that t and g1 are distinct for most traits.

Supplementary Note 1
In this note, we report how well GREML can capture the effects of the imputed transcriptome, noting that GREML does not explicitly model the effects of the imputed transcriptome. We simulated 500 replicates of phenotypic data using the genome-transcriptome model (Supplementary Table 1) under settings with and without random effects of the imputed transcriptome (! $ & = 0 vs. ! $ & = 0.4). For both settings, the covariance between random effects of the genome and those of the imputed transcriptome was set to zero. For each replicate, we fitted two models, a 'G model' that breaks phenotypic effects into random effects of the genome and residuals, i.e., y = g + ε, and a 'G-T model' that decomposes phenotypic effects into random effects of genome and of the imputed transcriptome and residuals, i.e., y = g + t + ε; both models were estimated using GREML. We declared the presence of transcriptomic effects when the G-T model had a better fit than the G model via a likelihood ratio test with one degree of freedom. We found that the type I error rate was controlled (0.042) under the zero transcriptomic effects setting. Under both settings, G-T model yielded unbiased estimates of model parameters (Supplementary Figure 6a). On the other hand, the G model produced unbiased estimates only under the null setting, as expected. Interestingly, in the presence of transcriptomic effects, the estimated genetic variance by the G model was equivalent to the sum of estimated variances due to the genome and imputed transcriptome by the G-T model (Supplementary Figure 6a), such that the phenotypic variance explained by the two models were similar. This indicates that SNP-based heritability estimates by GREML can be inflated by phenotypic variance due to the transcriptome. Importantly, although both models captured a similar amount of phenotypic variance, the fit of the G-T model was far better than that of the G model (Supplementary Figure 6b), indicating the partition of phenotypic variance represented by the G-T model is closer to the truth than that represented by the G model.