Mendelian Randomization using Public Data from Genetic Consortia

2.1 The Mendelian randomization ratio estimator

Suppose that a study or meta-analysis reports the results of regression analyses for each of m genetic variants, G_j, that are associated with their trait, X. These results might be in the form of the estimated regression coefficients, a_{X j}, and their variances, V_{X j}, or other statistics from which these quantities can be derived. It is important that the selection of the genetic variants is not based on the same data that is used to calculate a_{X j} for otherwise the estimated coefficients will be biased away from zero due to the Winner’s curse [10], so a_{X j} and V_{X j} might be taken from a replication study. These estimated regression coefficients will be modelled as,

where αXj is the true regression coefficient and the V_{X j} will be treated as known.

A second study or meta-analysis publishes similar data for the same variants but a different outcome, Y. The variants are unlikely to be top hits for Y so now we will need access to the full set of results in order to look-up the required estimates. The Winner’s curse is no longer a concern because the variants were not chosen for their effect in the second study. Suppose that the regression coefficients and variances from the look-up are b_{Y j} and the V_{Y j}, we can model them as,

where the βYj represent the true coefficients and the V_{Y j} are treated as known.

A Mendelian randomization for a continuous outcome targets the unconfounded regression coefficient, ϕ, of X on Y. Provided that the assumptions for Mendelian randomization hold for every genetic variant, there are m relationships βYj=ϕαXj, each of which creates a ratio, or Wald, estimator ϕˆj=bYj/aXj [11, 12]. These can be averaged with weights inversely proportional to their variances in order to create an overall Mendelian randomization estimate ϕˆ.

When the selected genetic variants are independent, we can estimate the variance of ϕˆ without needing to make assumptions about the unknown pattern of linkage disequilibrium and the coefficients αXj that apply to each variant separately will also be the coefficients in the joint regression of X on the genetic risk score S_Xi. That is,

where the confounder is omitted because it is assumed independent of each G_j and g_i j represents the measured genotype of the i^th subject for variant G_j coded as the number of effect alleles, 0,1 or 2. This genetic risk score S_Xi, assumes a per allele effect of each variant and ignores any interactions. Dominant or recessive genetic effects could be created but the necessary estimates of the coefficients are rarely published. In this model S_Xi represents the ideally weighted combination of the variants for use as a combined instrument in a Mendelian randomization.

We could estimate the variances of the ratio estimates of each variant using a Taylor series [13],

where the covariance term is omitted because X and Y come from different studies and their regression coefficients are independent. To estimate this variance we could just replace αXj and βYj by a_{X j} and b_{Y j}. The weights, w_j, needed for averaging the ϕj would then be the inverse of these variances and the resulting Mendelian randomization estimate of ϕ would have variance,

However, as we will see in the simulations, inverse variance weighting does not work well in this context.

Should we want to test the hypothesis that ϕ=0, we would need the variance when this hypothesis is true, that is when βYj=0. In that case the variance reduces to VYj/αXj2.

2.2 The ICBP estimator

The International Consortium for Blood Pressure Genome-Wide Association Studies [4] considered a subtly different question to Mendelian randomization. Their analysis takes the ratio estimates bYj/aXj and averages them using weights wj=[VYj/aXj2]−1. This produces,

We can think of this estimator either as an approximation to Mendelian randomization based on a simplified variance for ϕˆj that ignores the uncertainty in a_{X j}, or as the correct estimator of the regression coefficient on Y on the genetic risk score,

where SXi∗ differs from S_Xi because the a_{X j} have replaced the αXj. This ICBP estimator correctly addresses the question of what would be the regression coefficient if Y were regressed on the genetic risk score based on the particular study that provided X. Thus instead of estimating ϕ we would in fact be estimating the regression coefficient for the model,

and the actual value of this coefficient is,

where σj2 is the variance of G_j (see supplementary methods). This distinction between regression on S_X and SX∗ is blurred by the fact that, in both cases, the natural estimators for ϕ and ϕ∗ based on the individual variants are bY j/aX j, although the variances of these estimators do differ.

2.3 Bias adjustment

A problem arises with both the Mendelian randomization ratio estimator and the ICBP estimator because the sampling distribution of bY j/aX j has a noticeable skewness especially when the study measuring X is small or the true effect size αXj is small. Employing the delta method based on a Taylor series [13],

So we can obtain a less biased estimate of ϕ by replacing bY j/aX j with [bY j/aX j]/[1+VX j/aX j2].

2.4 Improved estimation of αXj

Much of the instability in MR estimates is due to the difficulty of estimating αXj accurately, because when the estimate, aXj, approaches zero the ratio, bY j/aX j, will become large and unstable. Were ϕ known, under the assumptions for Mendelian randomization, we would have two related relationships for each variant,

So a better estimate of αXj would be the weighted average,

Of course, ϕ is not known but we could create a two-stage procedure by estimating its value, for example by using the ICBP estimator, and then using that estimate in place of ϕ to improve the estimates of the αXj. Because Mendelian randomization is so sensitive to the accuracy of the estimates of αXj, basing those estimates on both X and Y may produce less bias and greater stability, provided the assumptions of Mendelian randomization hold for each variant.

2.5 Simulation

To investigate the properties of the different estimators, a simulation study was performed that was based on the model shown diagrammatically in Figure 1. The simulations were conducted twice for each scenario as if there were two independent but identically designed studies with the G_j and X taken from one study and G_j and Y were taken from the other. In each case we considered three sample sizes, 1,000, 5,000 and 20,000.

Figure 1:

The model for the simulations based on genetic variants G₁ to G_m combined into a risk score, S_X that drives and exposure X and possibly an outcome Y. The association between X and Y is confounded by an unobserved U.

In constructing the genetic risk score we used either 5, 10 or 50 independent variants. The minor allele frequencies of the variants were randomly selected to lie between 0.1 and 0.9, and the coefficients were adjusted so that each variant explained the same percentage of the variance in X. In the case of a score based on 5 genes, each gene explained 1 % of the variance and in the cases of 10 and 50 genes, each gene explained 0.5 % of the variance. So the scores based on 5 and 10 variants both explained 5 % of the variance in X and the 50 genes explained 25 %. The unconfounded effect of X on Y, ϕ, was varied between 0, 0.3, 0.6 or 0.9. The variances of X and Y were fixed at one so that for different ϕ,X explained 0 %, 9 %, 36 % or 81 % of the variance in Y. The confounder explained a third of the non-genetic variance of X and a third of the variance of Y that was not due to X. All scenarios were repeated 10,000 times.

3.1 Simulation

Table 1 summarizes the performance of the ICBP estimator [4] when the Mendelian randomization model holds and all genes act on Y through X. When we want to estimate the coefficient of the genetic risk score, SXi∗=∑j=1maXjgij, then the expected value of that coefficient varies depending on the results of the first study. Alternatively, we might want to estimate the regression coefficient for the regression of the genetic risk score, SXi=∑j=1mαXjgij, on Y in which case the true answer is always the value of ϕ used in the simulation. Regression on S_Xi is the appropriate analysis for a Mendelian randomization.

Table 1 shows that the estimation of ϕ∗ is very good, as one might expect since this is the situation that the ICBP estimator is designed to tackle. The coverage stays at its nominal level, the bias is small and RMSE decreases with sample size and with the percent of the variance explained by the genes. The results for 5 and 10 genes are very similar as both explain the same percentage of the variance of X, while the RMSEs for 50 genes are much smaller. The RMSEs for a sample size of 20,000 are about half those for a sample size of 5,000 consistent with negligible bias and an increase of four times in the sample size.

The ICBP results for estimating the MR coefficient, ϕ, used in the simulation, are less impressive. The coverage is correct for ϕ=0 but falls as ϕ increases. The bias is negative and increases in magnitude with ϕ and is especially noticeable with small samples. It decreases with sample size but increases with the number of genes. A bias (x1,000) of −50 for ϕ=0.3 or −100 for ϕ=0.6 or −150 for ϕ=0.9 would represent a 17 % error in estimation. In summary, the ICBP analysis performs well when estimating the regression coefficient for SXi∗ but is only approximate when used in a Mendelian randomization.

Table 2 compares the results of several different estimators of the MR coefficient, ϕ, for the middle sample size of 5,000. First we consider the use of the Taylor series variance estimate in an inverse-variance weighted analysis as described in section 2.1. When compared with Table 1, performance is actually worse than that for ICBP estimator, despite the improved variance estimation for the individual ratios. The problem lies in the correlation between the ratio estimates and the variance estimates. When the estimated ratio for a particular variant is randomly high, its estimated variance will also be larger. This correlation causes randomly large ratios to be down-weighted in an inverse-variance weighted analysis and so creates a bias towards zero. The effect of this correlation is removed if we use weights that do not depend on the estimated variances as in the second block of Table 2 that contains the results for a simple average with equal weights.

The third block of Table 2 shows further improvement by using the Taylor series adjustment to the estimate of the ratio as described in section 2.3 and the final block shows the result of using estimates of αXj that use information on both the gene-X and gene-Y relationships as described in section 2.4. In both cases the performance is improved.

3.2 Birth weight and glucose levels in adulthood

Horikoshi et al. published the results of a meta-analysis of genome-wide studies of birth weight [14]. They identified seven loci that replicated with a p-value below 5×10−8.Table 3 shows the replication results for the seven SNPs. The Meta-Analyses of Glucose and Insulin-related traits Consortium (MAGIC) have made their genome-wide results freely available on the internet [15] and in Table 3 we show the association between the lead SNPs for birth weight and fasting glucose level; similar data can be found in the supplementary table 5 of Horikoshi et al. [14]. The results were used to define a risk score for birth weight, which was used as the instrument in a Mendelian randomization by looking at its association with fasting glucose levels in adults.

Many epidemiological studies have found that low birth weight babies are at increased risk of diabetes [16] and if this relationship is causal we would expect that genes that are negatively related to birth weight would show a positive association with glucose levels and vice versa. The results in Table 3 do show such an inverse relationship.

Using the ICBP estimator ϕˆ=−0.155(se=0.032) with a pvalue=1.0×10−6. This is the correct p-value as it uses the standard error calculated under the null. When we are interested in producing confidence intervals for a MR estimate it is important to incorporate the extra uncertainty due to the estimation of the weights in the risk score. We have already noted that the inverse-variance weighted average performs badly because it is biased towards zero, here it gives, ϕˆ=−0.133(se=0.034). The simple average of the individual ratios is less prone to bias and gives ϕˆ=−0.186(se=0.044). Adjusting for the bias in the individual ratios bring the solution down slightly, ϕˆ=−0.176(se=0.044) and much the same effect is seen if a two-stage procedure is used to improve the individual estimates of αXj giving ϕˆ=−0.155(se=0.034). It is evident that the simple ICBP analysis performed well and although it does underestimate the standard error, that under estimation is slight.

The key question that this analysis does not address is whether the assumptions required by a Mendelian randomization hold for these variants. The evidence that they are truly associated with birth weight is strong. The likely confounders between birth weight and glucose level in adulthood relate to lifestyle and these are unlikely to be associated with these genes, so the only likely confounder is ethnicity. The meta-analysis of birth weight was conducted across populations of European origin and each meta-analysis adjusted for internal population stratification using genomic control, so confounding by ethnicity is unlikely to be a major problem.

The chief concern with the validity of this Mendelian randomization is pleiotropy. Biological knowledge about these genes is limited although, for instance, HMGA2 has previously been associated with height while ADRB1 has been associated with blood pressure and heart failure. Findings such as these suggest that the genes act through different pathways and so some of these genes might have secondary effects with a long-term influence on glucose levels. Genes that exhibit such pleiotropy would give different ratio estimates from those obtained from valid instruments. The ratios b_Y j/a_X j for the seven variants are, –0.05, –0.41, –0.02, –0.19, –0.29, –0.19, –0.15 each with an ICBP standard error of about 0.08. The difference between the second and third genes, ADCY5 and HMGA2, is 0.39 with a p-value of 4×10−4. Adjusting for the possible 21 pairs of genes, the Bonferroni adjusted p-value is 8×10−3. This is still significant and suggests that this analysis might be influenced by pleiotropy.