Empirical comparisons of multiple Mendelian randomization approaches in the presence of assortative mating

Abstract

Background

Mendelian randomization (MR) is widely used to unravel causal relationships in epidemiological studies. Whereas multiple MR methods have been developed to control for bias due to horizontal pleiotropy, their performance in the presence of other sources of bias, like non-random mating, has been mostly evaluated using simulated data. Empirical comparisons of MR estimators in such scenarios have yet to be conducted. Pleiotropy and non-random mating have been shown to account equally for the genetic correlation between height and educational attainment. Previous studies probing the causal nature of this association have produced conflicting results.

Methods

We estimated the causal effect of height on educational attainment in various MR models, including the MR-Egger and the MR-Direction of Causation (MR-DoC) models that correct for, or explicitly model, horizontal pleiotropy.

Results

We reproduced the weak but positive association between height and education in the Netherlands Twin Register sample (P= 3.9 × 10^–6). All MR analyses suggested that height has a robust, albeit small, causal effect on education. We showed via simulations that potential assortment for height and education had no effect on the causal parameter in the MR-DoC model. With the pleiotropic effect freely estimated, MR-DoC yielded a null finding.

Conclusions

Non-random mating may have a bearing on the results of MR studies based on unrelated individuals. Family data enable tests of causal relationships to be conducted more rigorously, and are recommended to triangulate results of MR studies assessing pairs of traits leading to non-random mate selection.

Introduction

Causal inference based on observational data is challenging because of genetic and environmental confounding. Mendelian randomization (MR) has tremendous potential in addressing causality in the presence of confounding, and currently enjoys considerable methodological and substantive interest (for example Hemani et al.¹, Verbanck et al.²). MR employs genetic variants associated with a given exposure as randomization instruments to detect the causal effect of the exposure on disease/outcome in non-experimental settings.^3,4 Specifically, MR exploits the fact that the genetic variants’ meiosis is random with respect to environmental and social confounders,^5,6 and they are not modified by disease outcomes (i.e. are not subject to reverse causality). The MR design yields information about causality, under the assumptions that the (genetic) instrument: (i) is well-associated with the exposure; (ii) is independent of confounders; and (iii) affects the outcome exclusively via the exposure variable (i.e. there is no horizontal pleiotropy). Although assumptions (i) and (ii) are likely to be tenable,^7,8 assumption (iii) is unlikely to hold on account of the ubiquity of pleiotropy.^9–12 As a result, multiple MR methods have been developed to control for bias due to horizontal pleiotropy.^8,13–17

Yet another source of bias that may have a bearing on the results of MR studies is assortative mating¹⁸ (i.e. ‘the tendency of like to marry like’). MR approaches assume that the genetic variants used to instrument the analysis are randomly distributed in the population. This assumption draws heavily on random mating,⁸ the absence of inbreeding and the absence of selection on the genetic instrument. There is extensive evidence for assortative mating, specifically for height and education (for example Courtiol et al.,¹⁹ Hugh-Jones et al.,²⁰ Yengo et al.²¹) as well as cross-trait,²² but also within (and between) psychiatric disorders like substance use.^23,24 Recently, Hartwig and colleagues¹⁸ showed by means of simulations that assortative mating may bias the causal estimate in two-sample methods that are robust to pleiotropy, like MR-Egger or inverse variance weighted regression. Empirical comparisons of MR estimators in such scenarios have yet to be conducted.

Here we assess the performance of several methods robust to horizontal pleiotropy in the presence of assortative mating. Pleiotropy and assortative mating have been shown to account equally for the genetic correlation between height and educational attainment.²² We therefore consider these traits for our illustrative analyses. Empirical evidence has shown that height predicts attained education at different ages throughout the lifetime.^25–28 There are several explanations for the height-educational attainment association.

The first is that stature has a causal effect on education.^22,29–35 The reverse hypothesis that educational attainment causally influences height is arguably too implausible to merit consideration (but see Supplementary material, available as Supplementary data at IJE online, for a sensitivity check). The third ‘underlying factor’ hypothesis is that the association is due to common factors influencing both traits. Twin studies^22,32 and studies in unrelated individuals which used measured genetic variants to disentangle the sources of this correlation³³ showed that the correlation between the two traits is entirely explained by shared genetic factors.^22,32,33 This genetic correlation can arise from vertical pleiotropy (i.e. the genetic effects on height are manifest in educational attainment due to the direct causal effect of height on education). The correlation may also be due to horizontal pleiotropy, i.e. due to overlap in the set of genes affecting brain size, intelligence and height, given the significant bivariate correlations observed between these traits (note that both intelligence and height are genetically correlated with brain volume^36–41; see also Figure S1, available as Supplementary data at IJE online).

Alternatively, cross-trait assortative mating,^18,22 dynastic effects and population structure⁴² may also induce genetic correlation between these traits. We consider these competing explanations in several MR models including MR-Egger¹⁴ and MR-direction of causation (MR-DoC) models^16,17 that correct for, or explicitly model, horizontal pleiotropy. We also use simulations to assess the MR-DoC model’s robustness in the presence of assortative mating and finally, we employ the latent causal variable (LCV) model⁴³ to account for potential bias arising from participant overlap⁴⁴ in the two-sample MR approaches.

Methods

 

Ethical approval

The NTR study was approved by the Central Ethics Committee on Research Involving Human Subjects of the VU University Medical Centre, Amsterdam, an Institutional Review Board certified by the US Office of Human Research Protections (IRB number IRB-2991 under Federal-wide Assurance-3703; IRB/institute codes, NTR 03–180). All subjects provided written informed consent.

Sample

We evaluated the association of height and educational attainment (EA) in a Dutch sample from the Netherlands Twin Register (NTR).^45,46 For the phenotypic association analysis, we selected 5899 Dutch adults (of whom 3184 were twins) with information on height and EA, relevant covariates (i.e. year of birth and sex) and observed genotypic data. EA was measured using a seven-category variable (1 = primary school only; 2 = lower educational schooling; 3 = lower vocational schooling; 4 = intermediate vocational schooling; 5 = intermediate/higher secondary schooling; 6 = higher vocational schooling; and 7 = university). The availability of repeated measures allowed us to run a series of consistency checks (see Supplementary material, available as Supplementary data at IJE online). Given the choice of the scheme used to code EA (1 through 7), and the fact that the underlying assumption of equivalent distances between adjacent EA categories is unlikely to hold in the case of ordinal data,⁴⁷ we converted the seven-category variable to sex- and birth-cohort specific ridit scores^47,48 as described in van Dongen et al.⁴⁹ The resulting scores reflect males’ and females’ hierarchy of attained education in accordance with the societal values and historical conditions that define a given birth cohort. In short, we calculated the sex- and cohort-specific ridit scores for each EA category j by summing the proportion of cases below that category and one half the proportion of cases in the category^47,48  Supplementary Table S1, available as Supplementary data at IJE online, shows descriptive statistics for the sample used in the analyses. Genome-wide single nucleotide polymorphisms (SNP) data were collected in the NTR sample on various genotyping platforms.⁴⁵ Missing data were cross-platform imputed with the Genome of the Netherlands as the reference set.^50–52 MaCH-Admix⁵³ was used for phasing and imputation. We refer to the Supplementary material, available as Supplementary data at IJE online, for details on quality control checks performed before and after imputation.

Statistical analyses

Phenotypic association

Education was regressed on height and on the following covariates: age, sex and birth cohort (see Supplementary material, available as Supplementary data at IJE online for details). To correct for the effect of family clustering, we employed a generalized estimation equation linear model⁵⁴ with an exchangeable model for the background familial covariance matrix and robust standard errors.^55–57

Causality testing: one- and two-sample MR

We tested the causal effect of height on education in one- and two-sample MR models. The one-sample MR was based on the two-stage least squares (TSLS) procedure (but standard MR can also be performed using the Wald estimator) and on the MR-DoC model.¹⁷ As an instrumental variable we used a polygenic score, i.e. a weighed linear combination of SNPs. The weights used in calculating the polygenic score were obtained from the largest to date GWAS of height⁵⁸ (excluding the NTR participants, see Supplementary material, available as Supplementary data at IJE online). Causality was tested in a subset of the sample comprising 3184 twins born between 1926 and 1989. The covariates used in the causality analysis were the dummy coded birth cohorts, age, sex, and 10 principal components (to correct for potential Dutch population structure). To increase the ethnic homogeneity of the sample, we removed ethnic outliers (see⁵⁹ and Supplementary material, available as Supplementary data at IJE online).

The two-sample MR was based on MR-Egger,¹⁴ the inverse variance weighted approach, mode-⁶⁰ and median-based MR methods.^13,15 We assessed the suitability of summary data for conducting MR-Egger using the I²_GX statistics.⁶¹ We also employed heterogeneity tests and Rucker’s framework⁶² to further assess pleiotropy and for model selection. To account for potential sample overlap that may bias the MR results, we used the latent causal variable (LCV) model⁴³ (see Supplementary material, available as Supplementary data at IJE online). The two-sample MR used summary statistics for the SNPs available for both traits, downloaded from Wood et al.⁵⁸ and from Okbay et al.,⁴¹ respectively. We refer to [https://gitlab.com/camish/empiricalcomparisons] for the R-script used in the analyses and to the Supplementary material, available as Supplementary data at IJE online, for more details on the methods employed.

Simulation: checking the robustness of MR-DoC in the presence of assortative mating

To evaluate the effect of assortative mating in the MR-DoC model, we created genotypes based on 150 diallelic loci (all minor allele frequencies = 0.5) in 50 000 parent pairs. Of the 150 loci, 50 affected height, 50 affected education and 50 had pleiotropic effects (i.e. affecting directly both height and EA). All 150 loci were used to compute polygenic scores, which were subsequently sorted ascendingly and used to divide the samples of mothers and fathers into five quintiles. Mate selection was contingent upon genetic similarity: males and females with the same quantile score were allowed to form parent pairs. To control the degree of assortative mating, we used a mover-stayer transition probability matrix to allow the males and females to switch to any other quintiles; the probability of switching to any other quantile equalled 0.05. Based on the resulting matings, we created two offspring per parent pair. The offspring were subsequently randomly assigned to pairs and the assortative mating procedure was repeated. After 10 generations, the familial (co)variances reached the expected equilibrium values⁶³ such that parent-parent genetic correlations equalled about 0.4 and the sibling correlations equaled about 0.7. Half of the resulting sibling pairs were used to create MZ twin pairs, by duplicating the first genotype. Phenotypic data were next created based on the genotypes and environmental variation. Polygenic scores computed using the 20 (horizontally) pleiotropic loci were used as the genetic instrument in the MR-DoC model.

Results

The observational regression analysis in the NTR sample showed that height is associated with educational attainment (Wald χ²(1) = 21.3; P = 3.9 × 10^–6; R² = ∼0.3%). TSLS indicated that height has a causal effect on educational attainment (b_causal = 0.004; standard error (SE) = 0.001; Wald χ²(1) =6.7; P  = 0.0096). Table 1 includes the results obtained using two-sample MR methods. The I²_GX statistic equalled 0.84—a value close to the recommended value of 1—providing no evidence of regression dilution, and supporting the use of MR-Egger. Consistent with the TSLS and with the other pleiotropy-robust methods results (see Table 1), MR-Egger implied a causal effect of height on EA. The MR-Egger intercept provided no evidence of detectable directional pleiotropy (resulting in the mean value of the pleiotropy distribution deviating from zero). In line with the MR-Egger intercept result, the funnel plot suggested the ‘no directional pleiotropy’ assumption is likely to be satisfied (Figure 1).

The Cochran’s Q test of heterogeneity⁶⁴ detected heterogeneity (P <10^–4). Here we assume this is due to balanced pleiotropy (i.e. the pleiotropy distribution has a zero mean); yet we note that heterogeneity may arise due to factors other than pleiotropy.⁶² The Rucker’s model selection framework (see Supplementary material, available as Supplementary data at IJE online, for details) showed there is no detectable directional pleiotropy, and indicated that an inverse variance weighted (IVW) random effects model adequately accommodated the heterogeneity present in the data (Q_diff(1) = 1.4; P = 0.2; where diff stands for the difference between the Cochran’s Q statistic—that tests for heterogeneity about the IVW model, and Rucker’s Q’ statistic—that tests for heterogeneity about the MR-Egger model). Consistent with these results, the MR-Egger intercept indicates that directional pleiotropy does not bias the estimate, and the Steiger test⁶⁵ (which removes SNPs with potential reverse-causal effects) supports the conclusion that causation flows from height to EA. It can be seen from Figure 2 that all approaches based on summary statistics yielded consistent results.

The LCV model (see Supplementary material, available as Supplementary data at IJE online, for details), which is robust to sample overlap, found a genetic correlation of ρ_g = 0.15, and a genetic causality parameter GCP = 0.27 (SE = 0.08; z-score = 2.93; one-tailed P = 0.002), but did not reach the GCP value of at least 0.6 to claim robust positive finding.

MR-DoC: simulation and empirical results

Simulation results showed that non-random mate selection had no detectable effect on the causal effect estimate in the MR-DoC model when the pleiotropic effect is freely estimated (i.e. parameter b₂ in Figure S2, available as Supplementary data at IJE online). The estimate of the causal effect (g₁ in Figure S2, available as Supplementary data at IJE online) was close to its true value of zero (the statistical test of the hypothesis g₁ = 0 was not rejected given alpha = 0.05). Interestingly, in the MR-DoC model with a non-pleiotropic polygenic score computed based on the 20 loci unique to height, non-random mating resulted in a pleiotropic effect (the statistical test of parameter b₂ = 0 was rejected given alpha = 0.05). This suggests that assortative mating is accommodated by the b₂ parameter.

The MR-DoC model¹⁷ explicitly estimates the amount of horizontal pleiotropy, but has to assume, for reasons of identification, that the unique environmental component influencing the exposure does not directly influence the outcome (other than through its effect on the exposure, i.e. re = 0, see Figure S2, available as Supplementary data at IJE online). With the parameter re constrained to 0, we tested the causal effect of height on EA while explicitly modelling horizontal pleiotropy (estimating the direct effect of the instrumental variable on the outcome). MR-DoC yielded a null finding: the causal effect was not different from zero at alpha = 0.05 (g₁ = 0.0019; SE = 0.002; standardized g₁ = 0.017; SE = 0.025; likelihood ratio test χ²(1) = 0.46; P = 0.49). When we assumed there is no horizontal pleiotropy (the parameter b₂ in Figure S1, available as Supplementary data at IJE online, was fixed to zero), we detected a causal effect of height on EA (g₁ = 0.004; SE = 0.001; standardized g₁ = 0.039; SE = 0.015; χ²(1) = 6.9, P = 0.0085). The test of the pleiotropic effect failed to reject the null hypothesis of ‘no pleiotropy’ (b₂ = 0.01; SE = 0.01; χ²(1) = 1.08; P = 0.29). Constraining b₂ to zero and testing whether re = 0 suggested that the assumption is tenable (χ²(1) = 1.08; P = 0.29). However, the null findings concerning absence of non-shared environmental confounding and absence of directional pleiotropy may also reflect a lack of power (see Supplementary material, available as Supplementary data at IJE online). Importantly, constraining the parameter b₂ to zero is not recommended in the light of ubiquity of pleiotropy across the human genome, and given that this parameter also accommodates potential non-random mating effects.

Discussion

We assessed the performance of several MR approaches in empirical data, in the presence of non-random mating. Recently, Hartwig and colleagues¹⁸ showed by means of simulations that non-random mating may bias the two-sample MR results based on summary statistics. Building on this work, the current study presented an empirical comparison of various one- and two-sample MR models when addressing the causal nature of the relationship between height and education (two traits extensively documented to lead to non-random mate selection). Previous studies addressing this relationship produced conflicting results. Tyrell et al.²⁹ suggested that the observational association between height and educational attainment (indexed by a continuous variable representing time spent in full-time education) is causal in nature. Using MR-Egger regression and a TSLS model in the UK Biobank sample, they showed that taller individuals complete a larger number of years of formal education. Conversely, based on the IVW MR,¹⁵ Hagenaars et al.⁶⁶ found no support for a causal effect flowing from height to education (indexed by a four-category variable representing degree level in the UK biobank). However, the authors noted that phenotype definition (i.e. using as the exposure a four-category variable instead of a continuous one as Tyrell et al. did), in conjunction with the use of weak genetic instruments may explain their null findings.⁶⁶

Assuming all SNPs meet the three standard instrumental variable assumptions (see Introduction), the current one-sample TSLS results implied that height has a small positive causal effect on EA. This conclusion was reinforced by several alternative MR approaches robust to certain forms of pleiotropy, like IVW, median- and mode-based approaches and MR-Egger. Whereas the mode- and median-based approaches can each handle certain proportions of pleiotropic SNPs, MR-Egger explicitly tests for directional pleiotropy. The MR-Egger test showed there is no detectable directional pleiotropy. It is, however, important to note that the power of the MR-Egger test to detect pleiotropy may be low (see Table 1 in Bowden et al.¹⁴). We used heterogeneity tests as additional means to assess pleiotropy. Although pleiotropy might be present, we found that this is likely to be balanced, and may be adequately modelled by a random effects IVW model.

Sample overlap may affect the conclusions of summary statistics-based MR approaches, with large overlap biasing the causal effect estimate towards the observed association (because of weak instruments). The LCV model is robust to sample overlap and to other factors that may bias two-sample MR methods, like heritable confounding, differences in polygenicity or unequal sample sizes (see O’Connor and Price⁴³ for details). This model provided weak evidence for the hypothesis that height affects education causally (genetic causality parameter GCP = 0.27; P = 0.002). Height and college completion are two of the traits studied by O'Connor and Price in their empirical applications in the UK Biobank.⁴³ They found a genetic correlation between these traits of ρ_g = 0.17 and a GCP of 0.33. The positive GCP value lends some support to the hypothesis of height affecting education, but these values are rather too small to claim a robust positive finding. O’Connor and Price recommended absolute GCPs values larger than 0.6 as unlikely to be false-positives, and concluded that this low GCP value likely reflects shared developmental pathways rather than a direct causal effect. Further consideration of the GCP statistic supports this recommendation: a heritable latent confounder affecting height with a path coefficient of 0.55 and education with a path coefficient of 0.3 could yield a genetic correlation of ρ_g = 0.17 and a partial genetic causality parameter of GCP = 0.33. If the former coefficient were higher (say, at least 0.8), then one could reasonably argue that, although height is not literally a cause of education, it is nevertheless a very good proxy for some other variable that is a true cause.

The MR-DoC model is also useful when the confounding is predominantly genetic. The MR-DoC model has not detected a causal effect of height on EA (χ²(1) = 0.46; P = 0.49), unless we assumed directional pleiotropy was absent. The test of the pleiotropic effect could not reject the null hypothesis that the effect equals zero—a result consistent with the MR-Egger test result. The test of the parameter b₂ suggests that the effect resulting from non-random mating is probably too weak to be picked up by the current sample. Notwithstanding this result, we argue against fixing this parameter to zero given the simulation results reported in Minică et al.¹⁶ (which showed that the estimate of the causal effect g₁ is unbiased provided the pleiotropic effect b₂ is freely estimated), and given the current simulation results (which showed assortative mating is accommodated by the b₂ parameter). Estimating the b₂ parameter might reduce power (if, indeed, pleiotropy is absent), but we rather err on the conservative side, given the ubiquity of pleiotropy across the human genome and given that non-random mating may presumably bias the causal effect if we do not model the directional pleiotropy path.

The usefulness of family data in addressing the non-random mating assumption was also demonstrated by Hartwig and colleagues.¹⁸ They showed, using data on trios, that height had no causal effect on education when correcting for possible effects of non-random mating (by using TSLS and adding the parental polygenic scores). However, as the authors acknowledge, it is worth noting that the small sample employed in this analysis may have resulted in imprecise estimates. The structural equation model developed by Warrington et al.,⁶⁷ to estimate jointly maternal and fetal effects on birthweight, may also be extended to provide another means to address the potential bias induced by non-random mating.

Family data can additionally address bias in MR studies arising from dynastic/genetic nurture effects that may also confound the association between height and education.^18,42,68 For instance, when available, trios data enable causality tests robust to dynastic effects;^18,68 while providing a straightforward means to address dynastic effects and non-random mating; this approach requires trios data and needs to assume absence of horizontal pleiotropy. The within-family-based approach described by Brumpton et al.⁴² may also address dynastic effects. These within-family estimates of SNP outcome and SNP exposure effects obtained in separate sibling samples can be used, in a second step, with pleiotropy-robust summary data estimators; yet, as the authors note, these family-based estimates are likely to be underpowered to allow for both family structure and pleiotropy. It is important to note that the MR-DoC model in some scenarios has to make the strong assumption that there is no appreciable non-shared environmental confounding¹⁶ However, the model does account for environmental confounding shared by twins. In so doing, MR-DoC also controls for dynastic effects because twins share the environments created by parents. As twins share a womb and a household/neighbourhood, often until age 18, this is non-trivial. Dynastic effects are not controlled for in MR on unrelated individuals.^69,70

To conclude, one should remain cautious when addressing the nature of the relationship between height and education as there are other confounders (beside horizontal pleiotropy) that may yield the observed association. This empirical example showed that non-random mating and, possibly, dynastic effects may generate data that are remarkably consistent with a causal model for the relationship between height and education. These factors—and others like fine-scale population structure^42,71—may affect this relationship, and pleiotropy-robust methods are not necessarily robust to these forms of confounding. Structural equation models⁶⁷ and MR models that exploit family data^16,18,42 should enable tests of causal relationships to be conducted more rigorously, and are recommended to triangulate⁷² results of MR studies assessing pairs of complex traits leading to non-random mate selection.

Funding

This work was supported by the National Institute on Drug Abuse (grant number DA-018673).

Acknowledgements

We are grateful to our reviewers for their excellent comments, which substantially improved the paper. We are grateful to the Netherlands Twin Register participants whose data we analysed in this study. We thank the participants in the GIANT and the SSAG consortia for making the GWAS summary statistics for height and educational attainment publicly available.

Conflict of interest

None declared.

References