Abstract
Background Mendelian randomization (MR) uses genetic variants as instrumental variables to investigate the causal effect of a risk factor on an outcome. A collider is a variable influenced by two or more other variables. Naive calculation of MR estimates in strata of the population defined by a variable affected by the risk factor can result in collider bias.
Methods We propose an approach that allows MR estimation in strata of the population while avoiding collider bias. This approach constructs a new variable, the residual collider, as the residual from regression of the collider on the genetic instrument, and then calculates causal estimates in strata defined by quantiles of the residual collider. Estimates stratified on the residual collider will typically have an equivalent interpretation to estimates stratified on the collider, but they are not subject to collider bias. We apply the approach in several simulation scenarios considering different characteristics of the collider variable and strengths of the instrument. We then apply the proposed approach to investigate the causal effect of smoking on bladder cancer in strata of the population defined by bodyweight.
Results The new approach generated unbiased estimates in all the simulation settings. In the applied example, we observed a trend in the stratum-specific MR estimates at different bodyweight levels that suggested stronger effects of smoking on bladder cancer among individuals with lower bodyweight.
Conclusions The proposed approach can be used to perform MR studying heterogeneity among subgroups of the population while avoiding collider bias.
Introduction
Mendelian randomization (MR) is the use of genetic variants as instrumental variables to assess the causal relationship between a risk factor and an outcome [1,2]. A valid instrumental variable (IV), or genetic instrument, must meet the following assumptions [3]: IV1, the instrument is associated with the risk factor; IV2, the instrument cannot affect the outcome directly, only potentially indirectly via the risk factor; and IV3, the instrument is not associated with any measured or unmeasured confounders (Figure 1A). If these assumptions are satisfied, an association of the instrument with the outcome is indicative of a causal effect of the risk factor on the outcome [1,4]. However, if either the IV2 or IV3 assumption are not satisfied, then the instrument could be associated with the outcome in the absence of a causal effect of the risk factor. However, only the IV1 assumption can be verified based on measured data [5].
Directed Acyclic Graphs (DAGs) illustrating relationships between the variables.
A) Mendelian Randomization causal diagram with the instrumental variable assumptions. The dashed lines between G and Y and between G and U, represent violations of the IV2 and IV3 assumptions respectively.
B) DAG considering a collider variable C, being a common child of genetic instrument G and confounders U. When conditioning on C (indicated by the square box on C), G and U become correlated (dashed line between G and U) and a violation of the IV3 assumption occurs.
C) DAG considering a collider variable C, being a common child of risk factor X and outcome Y.
D) DAG illustrating the variables and parameters used for the simulation study. Dash line from Y to C correspond to simulation scenarios B1 to B3
Collider bias can occur when conditioning on a collider, defined as a variable that is a common effect of two or more variables [6,7]. The presence of collider can be recognized in a causal diagram when there are two arrows pointing at the same variable; the node at which the arrowheads “collide” together is a collider. For example, in the standard MR diagram, the risk factor is a collider as it is affected by both the instrument and the confounders. Moreover, any variable that is a causal descendent of collider is also affected by the same variables and so is itself a collider; hence in MR any variable influenced by the risk factor is a collider (Figure 1B). Even if the variables influencing a collider are independent, they will typically become dependent when conditioning on the collider. Hence conditioning on a variable affected by the risk factor will typically generate a conditional association between the instrument and the confounders, violating the IV3 assumption, and biasing Mendelian randomization estimates of the risk factor on the outcome.
Selection bias is a form of collider bias that occurs when selection of individuals into a dataset is dependent on a collider. For example, when disease progression is considered as an outcome, only patients who have already developed the disease would be recruited into the study [6]. If risk of developing the disease is influenced by the risk factor, then it is a collider when considering disease progression as the outcome, and selection of the study sample would result in collider bias. Several papers related to selection bias in the context of IV analysis and MR have been already published [8–10].
Collider bias could also occur when stratifying the population based on a collider. As an example, we consider investigating the causal effect of the risk factor on the outcome for individuals with specific levels of a stratifying variable. Stratification is important for identifying whether there are subgroups of the population for which causal effects of the risk factor are different, and so the outcome would be affected more strongly by an intervention on the risk factor. However, if the stratifying variable is a collider, an association between the instrument and the outcome in strata of the population could arise due to collider bias, invalidating the results. In particular, collider bias could affect some estimates more than others, leading to heterogeneity in the stratum-specific causal estimates even if the true causal effect is the same across strata.
The aim of this paper is to present an MR approach that obtains estimates in strata of the population that do not suffer from collider bias. The structure of this paper is as follows: first, we demonstrate the bias that arises from conditioning on a collider; second, we propose an approach to calculate MR estimates in strata of the population and evaluate heterogeneity between estimates in the different strata; third, we illustrate this new technique in simulation studies and an applied example using the UK Biobank resource; and finally, we discuss the interpretation of estimates and limitations of the approach.
Methods
Illustration of collider bias
The simplest MR method to estimate the causal effect of a risk factor X on outcome Y with a genetic instrument G is the ratio method [2]. With a single instrument, a continuous risk factor and outcome, and under assumptions of linearity and no effect modification, the ratio estimate is defined as: , where
is the coefficient from regressing Y on G, and
is the coefficient from regressing X on G [11].
Collider bias will occur when adjusting for a collider variable C in the regression models for the ratio estimate, since an association between the instrument and the outcome will occur through conditioning on the collider. To demonstrate the impact and magnitude of collider bias, we performed a simulation study in which we compared estimates when no adjustment on C is made versus when the outcome regression is adjusted for C. It is also possible to adjust the risk factor regression for C; however, while this will distort estimates, this adjustment alone will not bias causal estimates when the true causal effect is null.
Stratification in Mendelian randomization
To further illustrate the impact of collider bias, we performed a simulation study in which we calculated causal estimates using the ratio method within strata of the population defined using a variable that is influenced by the risk factor, and hence is a collider. We compared two approaches: first, we stratified directly on the collider C, and second, we stratified on a new variable C0, referred to as the “residual collider”. The residual collider was generated as the residual from regression of the collider on the genetic instruments:
The residual collider C0 is not associated with the instrument, and hence it is not itself a collider. It is influenced by the component of the risk factor that is not a function of G, but not by the component that is a function of G. However, provided that the genetic instrument does not explain much of the variance in the risk factor (as is typical in a MR application), it is likely not to explain much of the variance in the collider, and so the residual collider will be highly correlated with the collider. Hence, while stratifying on the residual collider is important to avoid bias, the strata defined by stratifying on the collider or residual collider are likely to be similar and so any difference in the interpretation of stratum-specific estimates is minimal.
Here we considered estimates in four strata of the population defined by quartiles of the distribution of the collider or residual collider; however, in practice any number of strata could be considered. We estimated genetic associations with the outcome in each stratum separately. We estimated genetic associations with the risk factor in the full dataset, although if it is believed that these associations vary between strata, it would be possible to estimate these within each stratum as well. The stratum-specific estimate is calculated as the ratio of the stratum-specific genetic association with the outcome divided by the genetic association with the risk factor. We also investigated heterogeneity between the stratum-specific estimates using Cochran’s Q statistic [12], and (in the applied example) we examined the presence of a trend in the estimates by meta-regression of the stratum-specific estimates on the median value of the collider in each stratum [13].
Simulation set-up
To investigate the impact of collider bias in realistic scenarios, we generated simulated data using the following data-generating model:
We simulated the instrument G, the confounder U, and the error terms for X, Y and C, εX, εY and εC, as independent normally distributed variables. The risk factor X is defined as a linear combination of the instrument, the confounder, and the error term. The outcome Y and the collider C are both linear combinations of the risk factor, confounder, and their error terms. In each simulated dataset, we also generated the residual collider C0 as the residual from regression of C on G as previously described.
The causal estimate of interest is β1, while α2 and β2 represent the effects of U on X and Y respectively; α1 is the effect of G on X; and μ1 and μ2 are the effects of X and U on C, respectively.
We considered three scenarios based on the parameter β1: Scenario A1, where there is a null causal effect of X on Y (β1 = 0); Scenario A2, where the effect is constant and positive (β1 = 0.5); and Scenario A3, where the effect depends on C (β1 = 0.5 + 0.2C). In Scenario A1, we considered estimates from the ratio method with and without adjustment for the collider. In Scenarios A2 and A3, we consider stratum-specific estimates from stratification on the collider C or the residual collider C0.
We varied the other parameters to consider the impact of different settings on collider bias: i) α1 = (0.05, 0.1, and 0.3), in order to study the impact of the strength of the instrument on estimates; ii) positive confounding (α2 = 0.8, β2 = 0.8) negative (α2 = −0.8, β2 = −0.8) and mixed (α2 = 0.8, β2 = −0.8), to study how the direction of confounding affects the estimates and, iii) μ1 and μ2 = (−1, −0.5, 0, 0.5, 1) to study how the strength of the collider effects influence bias.
We also considered scenarios where the collider is a common effect of X and Y (Figure 1C). In these scenarios, the collider is generated as C = μ0 + μ1X + μ2U + μ3Y + εC, where μ2 = 0.3 and μ3 = (−1, −0.5, 0, 0.5, 1). In Scenario B1, the causal effect of X on Y is null (β1 = 0), in Scenario B2, the causal effect is constant and positive (β1 = 0.5), and in Scenario B3, the causal effect depends on U (β1 = 0.5 + 0.2U), as it is not possible for the causal effect to depend on C when C is a function of Y. Finally, we investigated additional scenarios with a binary outcome Y. We generate Y from a Binomial distribution where the probability is obtained from a logit transformation as: logit(P(Y = 1)) = β0 + β1X + β2U, where β0 = 0.5. In Scenario C1, the causal effect of X on Y is null (β1 = 0), in Scenario C2, the causal effect is constant and positive (β1 = 0.5) and in Scenario C3, the causal effect depends on C (β1 = 0.5 + 0.2C). In the binary outcome scenarios, genetic associations with the outcome were estimated by logistic regression. For these additional scenarios, we only consider α1 = 0.1 and the positive confounding values; otherwise, we consider all parameters as in scenarios A1 to A3.
We considered a sample size of n = 10,000 and m = 500 replications for each set of parameter values. A directed acyclic graph illustrating the simulation parameters is shown in Figure 1D.
Applied example: effect of tobacco smoking on bladder cancer risk across bodyweight strata
We applied the proposed MR stratification approach to investigate the causal effect of tobacco smoking on bladder cancer across strata of the population defined by bodyweight. Tobacco smoking is one of the strongest risk factors for cancer, and it has already been reported to be causally associated with bladder cancer risk in a previous Mendelian randomization study [14]. With our current example, the objective was to investigate whether the effect of smoking on the risk of developing bladder cancer is homogeneous across the bodyweight distribution of the population, while avoiding potential collider bias by applying our new stratification approach.
We performed analyses in the UK Biobank study, a population-based cohort of more than 500,000 United Kingdom residents recruited between 2006 and 2010 [15]. For our analysis, we restricted to unrelated European ancestry participants, resulting in a final sample size of 367,643 individuals following sample selection and quality control procedures as described previously [14]. The risk factor is a binary variable representing the smoking behaviour, defined as being a current smoker versus a former or never smoker; the stratifying variable is bodyweight, measured in kg; and the binary outcome is bladder cancer status, defined based on the data from national registries (International Classification of Diseases 9th edition codes: 188, 189.1, 189.2, V10.51, V10.53; or International Classification of Diseases 10th edition codes: C67, C65, C66, Z85.51, Z85.54, Z85.53), and self-reported information from an interview with a nurse practitioner. The instrument for smoking was a weighted genetic risk score comprising 378 conditionally independent SNPs obtained from a genome-wide association study (GWAS) assessing associations with smoking initiation (i.e., probability of ever smoked regularly), and weighted by the associations with smoking initiation [16]. Genetic associations with the risk factor and outcome were obtained by logistic regression in UK Biobank with adjustment for age, sex, and 10 genomic principal components. While age, sex, and principal components cannot logically be colliders as they are not affected by the risk factor or outcome, bodyweight is likely to be a collider, as it is influenced by smoking status [17].
Results
Illustration of collider bias
Results from Scenario A1 (β1 = 0, null causal effect) are presented in Table 1 for α1 = 0.1 (corresponding to R2=0.006 for the mean proportion of variance in the risk factor explained by the instrument and a mean F statistic of 60.8) and Supplementary Tables 1 and 2 for α1 = 0.3 (corresponding to R2=0.051, mean F statistic of 548.6) and α1 = 0.05 (corresponding to R2=0.001, mean F statistic of 15.3). In each case, we report the median estimate of β1 across simulations, and the empirical type I error rate, representing the proportion of simulated datasets where the 95% confidence interval for the ratio estimate excludes zero. With no adjustment for the collider, median estimates were close to zero and empirical type I error rate was close to the expected value of 5%. When adjusting for the collider in the regression of Y on G, estimates were biased, and type I error rates were substantially above 5%. The only exception was for μ1 = 0; in this case, the variable C is not a function of the risk factor, and so does not act as a collider. Bias and type I error rates generally increased for more extreme values of μ1 and μ2 (both positive and negative values). The direction of bias depended on μ1 and μ2 and the direction of confounding.
Median of β1 estimates and empirical Type I error rates for Scenario A1 (null causal effect, β1 =0) with positive, negative, and mixed confounding, and α1 =0.1
Stratification in Mendelian randomization
Results from Scenario A2 (β1 = 0.5, constant positive effect) are presented in Table 2 for α1 = 0.1 with positive confounding. Supplementary Table 3 shows results for α1 = 0.1 with negative and mixed confounding, and Supplementary Tables 4 and 5 for α1 = 0.3 and α1 = 0.05. We report the median estimate of β1 in four strata of the sample defined by quartiles of the collider C or residual collider C0, and the proportion of simulated datasets for which the heterogeneity test statistic is rejected. When stratifying on the collider, median estimates were somewhat variable between the strata, although the proportion of datasets in which the heterogeneity test rejects the null hypothesis of homogeneity was not much above 5% in any scenario, reaching a maximum of 11% when α1 = 0.3. However, if we considered stronger instruments or larger sample sizes, we would see this proportion considerably exceed 5% (see Supplementary Table 6 where we first set α1 = 0.5 and n=10,000, and then set α1 = 0.1 and n = 50,000, and the type I error rate reached 16% in each case). Median estimates differed substantially from the true value of 0.5 across strata, especially when the collider was strongly affected by the risk factor. In contrast, when stratifying on the residual collider, median estimates of β1 were close to 0.5 throughout, and there was no suggestion in any case that the heterogeneity test rejected the null above the expected 5% rate.
Median of causal estimates in different quartiles, and proportion of datasets in which the homogeneity test was rejected for Scenario A2 (fixed causal effect of β1 =0.5) with positive confounding and α1 =0.1
Results from Scenario A3 (variable effect) are presented in Table 3 for α1 = 0.1 with positive confounding. Supplementary Table 7 shows results for α1 = 0.1 with negative and mixed confounding, and Supplementary Tables 8 and 9 for α1 = 0.3 and α1 = 0.05. Estimates differed somewhat when stratifying on the collider versus the residual collider, although in both cases median estimates increased across the four strata. The proportion of datasets in which the heterogeneity test was rejected, which in this case represents the empirical power to detect heterogeneity in the stratum-specific estimates, was consistently higher when stratifying on the residual collider, indicating that true differences in the stratum-specific estimates were better detected when stratifying on the residual collider.
Median of causal estimates in different quartiles, and proportion of datasets in which the homogeneity test was rejected for Scenario A3 (varying causal effect) with positive confounding and α1 =0.1
Additional scenarios
In Scenarios B1 (β1 = 0), B2 (β1 = 0.5) and B3 (β1 = 0.5 + 0.2U), where the collider was a function of both the risk factor and outcome, similar results were observed, with collider bias evident when conditioning on the collider (Supplementary Table 10) and when stratifying on the collider (Supplementary Table 11). Collider bias in Scenarios B1 and B2 was greater compared with Scenarios A1 and A2 where the collider was a function of the risk factor only. Similarly, bias was not observed when stratifying on the residual collider (Supplementary Table 11). For Scenario B3, the power of the homogeneity test was lower in comparison to Scenario A3 (Supplementary Table 11), as the dependence of effect heterogeneity on the collider was weaker; however, heterogeneity was detected more often when stratifying on the residual collider than on the collider.
For Scenarios C1 (β1 = 0), C2 (β1 = 0.5) and C3 (β1 = 0.5 + 0.2C), where the outcome was binary, again similar results were observed, with collider bias evident when conditioning on the collider in Scenario C1 (Supplementary Table 12) and when stratifying on the collider in Scenarios C2 and C3 (Supplementary Table 13). Bias was smaller than in cases with a continuous outcome, although direct comparison is somewhat unfair as estimates with a binary outcome were obtained from logistic regression and so represent log odds ratios. Estimates when stratifying on the residual collider were slightly attenuated from 0.5 due to the non-collapsibility of the odds ratio [18,19]. Despite this, in Scenario C2 we observed similar estimates across the different strata of C0 for each set of parameter values. Similarly, in Scenario C3 we observed that median stratum-specific estimates increased across the four strata when stratifying on either the collider or residual collider. Power to detect heterogeneity was lower compared with Scenario A3 as the stratum-specific estimates are less precise, although again power was consistently higher when stratifying on the residual collider.
Applied example: effect of tobacco smoking on bladder cancer risk across bodyweight strata
Estimates for the causal effect of smoking on bladder cancer in strata of bodyweight and residual bodyweight are shown in Table 4. Estimates represent the odds ratio for bladder cancer per one unit increase in the log odds of being a current smoker. Estimates were positive in all strata, although larger in strata 1 and 2 for both bodyweight and residual bodyweight, and 95% confidence intervals excluded the null in these strata only. Although the homogeneity test was not rejected for either collider variable (p-value = 0.151 and p-value = 0.084 for bodyweight and residual bodyweight, respectively), there was evidence of trend in the stratum-specific estimates for residual bodyweight from meta-regression on the mean value of bodyweight in each stratum (p-value = 0.019). These results suggest that the effect of smoking on bladder cancer is stronger for subgroups of the population with lower bodyweight.
Applied example using UK Biobank to investigate the effect of smoking status on bladder cancer risk in different bodyweight strata.
Discussion
In this paper, we have demonstrated that conditioning or stratifying on a variable that is a collider can have a serious impact on MR estimates. We have introduced a simple approach that constructs a new variable, the residual collider, which is typically highly correlated with the collider, but is independent of the instrument. Estimates obtained from stratification on the residual collider did not suffer from bias in a range of simulation studies. Stratification on the residual collider allows investigators to explore causal estimation in relevant subgroups of the population. We applied our new approach to demonstrate that MR estimates for the effect of smoking on bladder cancer differ within strata of bodyweight, suggesting that the effect of smoking is stronger for subgroups of the population with lower bodyweight.
The approach of stratifying on the residual collider follows the same logic as a previously proposed method for non-linear MR, in which causal estimates are obtained in strata of the population defined by the “residual risk factor” or “IV-free exposure” [20,21]. This variable is defined similarly to the residual collider, except the collider variable is the risk factor itself. This method has been used previously to estimate the causal effect of blood pressure on coronary heart disease risk within strata of blood pressure, resulting in a curve that represents the shape of the causal relationship between the risk factor and the outcome [22]. This paper extends on that method, showing that the same idea can be used to provide causal estimates stratified on a separate variable even if that variable is a collider.
There are some limitations to this approach. First, while the independence of the residual collider from the instrument is theoretically justified, we demonstrated the validity of our approach through simulation studies. Although we considered a range of different scenarios and parameter values, it is not possible to consider every possible data-generating mechanism by which that a collider could arise. Second, in practice, the relationships between variables are unknown, and so it may be unclear whether a proposed stratifying variable is a collider. However, even if the variable is not a collider, it is unlikely stratification on the residual variable will lead to invalid estimates, suggesting that this approach would be valid for stratifying on variables that are not colliders. This was demonstrated in the simulation study when the effect of the risk factor on the “collider” was zero (μ1 = 0), and so the stratifying variable was not a collider. One exception is if the stratifying variable is on the causal pathway from the risk factor to the outcome. Stratification on such a variable (a “mediator”) will lead to biased estimates even in the proposed approach. Finally, the degree of collider bias depended on the strength of the effects of the risk factor and confounder on the collider, and the direction of confounding. It is possible that collider bias may not be substantial in practice, as observed in the applied example, where estimates were broadly similar when stratifying on bodyweight or residual bodyweight. However, the power to detect heterogeneity in stratum-specific estimates in the simulation study was greater when stratifying on the residual collider, especially when the proportion of variance of the risk factor explained by the instrument was higher. This was also observed in the applied example, where a lower p-value was observed in both the heterogeneity test and the trend test when stratifying on residual bodyweight.
The finding that the effect of smoking on bladder cancer is greater in lower bodyweight subgroups is plausible, because for any given level of cigarette consumption smaller individuals will tend to be exposed to greater concentrations of carcinogens [23]. An alternative explanation is that the genetic variants could associate more strongly with smoking intensity in individuals of lower bodyweight. However, we would be cautious not to interpret estimates in the higher bodyweight quartiles as implying an absence of a causal effect in heavier individuals; it is possible that the null estimates reflect limited power. A limitation of the applied example is overlap between the discovery dataset for the genetic variants, and the dataset used in the MR analysis, which can lead to winner’s curse, and the one-sample setting, which can lead to weak instrument bias.
In conclusion, we recommend that researchers performing MR to investigate causal effects in strata of a population defined by a collider stratify on residual values of the collider rather than stratifying on the collider directly.
Funding
SB is supported by a Sir Henry Dale Fellowship jointly funded by the Wellcome Trust and the Royal Society (204623/Z/16/Z). This research was funded by United Kingdom Research and Innovation Medical Research Council (MC_UU_00002/7) and supported by the National Institute for Health Research Cambridge Biomedical Research Centre (BRC-1215-20014). The views expressed are those of the authors and not necessarily those of the National Institute for Health Research or the Department of Health and Social Care.
The work was partially supported by Fondo de Investigaciones Sanitarias (FIS), Instituto de Salud Carlos III, Spain (#PI18/01347), Ministerio de Ciencia e Innovación, Spain (#PID2019-104681RB-I00). CC received a fellowship from CIBERONC, Insitituto de Salud Carlos III, Spain.
Competing interests
DG is employed part time by Novo Nordisk, outside of the submitted work