Non-linear Mendelian randomization: detection of biases using negative controls with a focus on BMI, Vitamin D and LDL cholesterol

Mendelian randomisation (MR) is an established technique in epidemiological investigation, using the principle of random allocation of genetic variants at conception to estimate the causal linear effect of an exposure on an outcome. Extensions to this technique include non-linear approaches that allow for differential effects of the exposure on the outcome depending on the level of the exposure. A widely used non-linear method is the residual approach, which estimates the causal effect within different strata of the non-genetically predicted exposure (i.e. the “residual” exposure). These “local” causal estimates are then used to make inferences about non-linear effects. Recent work has identified that this method can lead to estimates that are seriously biased, and a new method—the doubly-ranked method—has been introduced as a possibly more robust approach. In this paper, we perform negative control outcome analyses in the MR context. These are analyses with outcomes onto which the exposure should have no predicted causal effect. Using both methods we find clearly biased estimates in certain situations. We additionally examined a situation for which there are robust randomised controlled trial estimates of effects—that of low-density lipoprotein cholesterol (LDL-C) reduction onto myocardial infarction, where randomised trials have provided strong evidence of the shape of the relationship. The doubly-ranked method did not identify the same shape as the trial data, and for LDL-C and other lipids they generated some highly implausible findings. Therefore, we suggest there should be extensive simulation and empirical methodological examination of performance of both methods for NLMR under different conditions before further use of these methods. In the interim, use of NLMR methods needs justification, and a number of sanity checks (such as analysis of negative and positive control outcomes, sensitivity analyses excluding removal of strata at the extremes of the distribution, examination of biological plausibility and triangulation of results) should be performed. Supplementary Information The online version contains supplementary material available at 10.1007/s10654-024-01113-9.


Introduction
Mendelian randomization (MR) is now a widely used epidemiological technique contributing to improved causal inference in settings in which unmeasured confounding may bias estimates [1][2][3].MR findings have been consistent with randomised trial data in many settings [4], and have assisted causal inference as a component of triangulation of evidence [5,6] Recently, there have been developments in MR in which attempts are made to extend causal inference from a whole population (e.g. the overall effect of increasing Vitamin D on all-cause mortality), to strata of that population (e.g.only those with low levels of Vitamin D) [7][8][9][10][11].This has the potential to answer important questions around the linearity of exposure-outcome effects.The development of the "residual" method led to MR studies that purported to show striking non-linear effects [12][13][14][15][16][17][18], in particular with respect to Vitamin D and BMI as exposures.However, recent analyses have suggested that this method may be susceptible to serious bias, including that likely related to amplification of selection bias within strata [3,[19][20][21].In this manuscript, we explore a number of issues that suggest caution in using non-linear approaches at present.

Development of the residual method and concerns about inconsistent causal estimates
In the residual method the measured levels of the exposure (e.g.Vitamin D) are regressed onto the genetic instrument (e.g. a polygenic risk score for Vitamin D).The residuals of that regression represent the IV-free exposure (sometimes referred to as the "non-genetic" portion of Vitamin D).These residuals are then split into strata (often 10, or 100 strata) ordered by the mean IV-free exposure value.Therefore, the sample in stratum 1 have the lowest level of "IV-free" Vitamin D, and the sample in stratum 10 (or 100) have the highest value.MR is then run within each stratum, and estimates are combined by parametric statistical methods (such as fractional polynomials) to estimate a non-linear curve.The requirement for stratification based on the residuals rather than the exposure itself is driven by the desire to avoid collider bias [22] as this would be induced by stratification on the exposure itself.
There has been a recent focus on the performance of this method.In particular, one study identified a non-linear predicted causal effect of Vitamin D on a number of outcomes, suggesting that increasing Vitamin D in those who were below the median level might produce considerable benefit (as opposed to the null overall causal effect) [18].In some analyses, despite a precisely estimated null or harmful overall effect, there was a protective effect within each stratum, with some effects substantial.Following the journal and the authors being alerted to this [19] the authors referred to it "as a logical impossibility if all estimates have a causal interpretation" [23] after the journal retracted the paper [24].Several other analyses applying identical methods to the same exposure, with overlapping data and outcomes, have, unsurprisingly produced similar and likely equally spurious findings [15,16].
Following the documentation that the findings were erroneous the authors of the initial non-linear Vitamin D analysis [23] re-evaluated the residual method used in the study.In a recent paper [25] they showed via simulation studies that the residual method was biased under certain situations, in part relating to the "constant genetic effect" assumption-that the genotype-exposure effect is linear and thus constant across strata (i.e.For Vitamin D, the constant genetic effect assumption was shown not to be satisfied, with much larger genotype-exposure effect estimates in top strata than in the bottom strata [25].They also noted that this heterogeneity in IV-exposure association was not identified in the residual method, and suggest that this is a flaw in the residual approach.It is not clear why this would have led to the J-shaped curve observed in the residual method, however.To address this challenge a different way of stratifying the exposure was developed, called the "doubly-ranked" method [11]. In this method, the strata are developed in multiple steps (Fig. 1).In the first step, the population is ranked by the level of the genetic instrument into pre-strata.Subsequently, within each pre-stratum, the participant with the lowest level of the exposure is taken and placed into the lowest final stratum, (e.g.strata 1).The participant with the next lowest level of the exposure in pre-strata 1 is placed into final stratum two.The same process occurs for all pre-strata.The first stratum therefore contains the first participant of every pre-strata, each of which has the lowest level of exposure in the pre-stratum.The number of pre-strata (K) to be generated depends on the number of strata desired (J) and the total sample size (N) using the formula N = J x K.That is, for a sample size of 1000, and to generate a final desired number of strata of 10, 100 pre-strata are generated.
Fig. 1 A worked example of strata generation in the doubly-ranked approach.In the first step, all participants are sorted by their level of the instrumental variable (IV) into pre-strata.In the second step, participants are sorted by their level of exposure within each pre-stratum.In the third step, the first participant from each pre-stratum is placed in the first final stratum.The second participant in each pre-stratum is placed in the second final stratum, and so on.In the orange, a single participant is followed.This participant has the lowest level of IV so is placed in pre-stratum 1.They have the second lowest level of exposure within this pre-stratum so are placed in the second final stratum A more detailed exposition is given elsewhere [11].The doubly-ranked method is suggested (using simulation studies) by the authors to be insensitive to changes in the exposure to genotype association over strata, i.e. does not rely on the constant genetic effect assumption to generate unbiased estimates.The retracted non-linear Vitamin D paper [18]which claimed substantial non-linearity in the causal effect of Vitamin D-has been replaced by a revised paper applying the doubly-ranked approach which suggests no meaningful non-linearity [23,24].
Despite the development of the doubly-ranked method solving some of the issues evident in both empirical and simulation studies of the residual method, there remain a number of uncertainties about the situations in which each method provides unbiased causal estimates.In this paper we evaluate the residual and doubly-ranked methods using "negative control outcome" (NCO) analyses [5,6,26] in UK Biobank, a large and widely researched cohort study.NCO analyses are analyses performed on control outcomes that should lead to a null result, and where identification of a non-null estimate suggests bias.We additionally examined a situation for which there are robust randomised controlled trial estimates of effects-that of the effect of low density lipoprotein cholesterol (LDL-C) reduction on myocardial infarction, where randomised trials have provided strong evidence of the shape of the relationship [27,28].We compared these estimates from trials with those generated by both non-linear MR methods.We present our findings in the context of a recent paper from the principal developer of the doubly-ranked method, and draw attention to highly implausible findings in this report [29].

Data source
This study was performed in UK Biobank [30].UK Biobank recruited around 500,000 participants between 2006-2010 from 22 sites across the UK.Participants were invited via post and had a range of interviews on recruitment to record previous life events, demographics, and medical history.Additionally, participants had blood samples taken for biochemical testing and had a range of physical and anthropometric measures performed (e.g.body mass index, blood pressure).Subsequently, participants had record linkage to electronic health record data for secondary care and national death data (for > 99% of participants).
We analysed data on the subset of participants who were minimally related and of White British ancestry (n = 385,290), defined by relationship analyses [31].Specifically, we used the MRC Integrative Epidemiology Unit genetically quality controlled data, described in detail elsewhere, excluding highly related participants [31].

Exposures and data source
This study focussed on Vitamin D (measured as nmol/L) BMI (measured as kg/m 2 ), and LDL-C and triglycerides (TG) (both measured as mmol/L) as exposures, as all of these were measured in participants on recruitment to UK Biobank, and have been widely used in non-linear MR analyses [17,23,29,32].Details on measurement are available via the UK Biobank website.Our study sample for each analysis included only participants who had these values recorded successfully on recruitment to UK Biobank and had genetic data available.Sample sizes differed slightly due to missing values for each exposure (e.g.due to assay failure, inability to give blood, etc.).For triglycerides, we also extracted data from a subsequent visit, around 2 years later, for around 15,000 participants.
We then identified genetic variants to use in MR for all exposures.To identify variants associated with BMI, we used a large recent meta-analysis of BMI which did not include any participants of UK Biobank [33].We extracted all variants which were associated with BMI (p < 5 × 10 −8 ) and then clumped them by linkage disequilibrium to include only independent variants (r 2 < 0.001, kb = 10,000).LD clumping was performed using the TwoSampleMR package in R and was based on LD from the European ancestry participants from the 1000 Genome's project [34].In total we included 68 variants (Supplementary Table S1).For Vitamin D, to match as close as possible to a previously described non-linear MR analysis, we used the same set of variants that they used [23].These have previously identified to be associated with Vitamin D levels in meta-analysis from four well characterised genomic regions known to be involved in Vitamin D metabolism (GC, DHCR7, CYP2R1, CYP24A1) [35].In that analysis, 21 variants were used but 3 were unavailable in our data so in total 18 variants were used to generate our polygenic risk score (Supplementary Table S2).The three missing variants were all rare (minor allele frequency < 0.005) in 1000 Genomes European subset) and were therefore unlikely to substantially affect our instrument.For LDL-C and TG we used summary statistics from the most recent analysis of the Global Lipids Genetics Consortium (GLGC) excluding UK Biobank participants.Again, we took independent variants (p < 5 × 10 -8 , r 2 < 0.001, kb = 10,000) from across the genome.In total we included 296 variants for LDL-C and 366 variants for TG (Supplementary Table S3).
To generate an individual PRS for each participant we used PLINK v2.0.4 [36], with the summary effect generated by summed weighted allele scores, with each allele weighted by its effect on the exposure.

Genotyping and imputation
This analysis was performed on quality controlled genetic data held within the MRC Integrative Epidemiology Unit.Full details of quality control are available elsewhere [31] but a summary is provided here.The full data release for this study contains the cohort of successfully genotyped samples (n = 488,377).49,979 individuals were genotyped using the UK BiLEVE array and 438,398 using the UK Biobank axiom array.Pre-imputation QC, phasing and imputation are described elsewhere [30].In brief, prior to phasing, multiallelic SNPs or those with MAF ≤ 1% were removed.Phasing of genotype data was performed using a modified version of the SHAPEIT2 algorithm [37].Genotype imputation to a reference set combining the UK10K haplotype and HRC reference panels was performed using IMPUTE2 algorithms [38].The analyses presented here were restricted to autosomal variants using a graded filtering with varying imputation quality for different allele frequency ranges.Therefore, rarer genetic variants are required to have a higher imputation INFO score (INFO > 0.3 for MAF > 3%; INFO > 0.6 for MAF 1-3%; INFO > 0.8 for MAF 0.5-1%; INFO > 0.9 for MAF 0.1-0.5%)with MAF and INFO scores having been recalculated on an in-house derived 'European' subset.Data quality control Individuals with sex-mismatch (derived by comparing genetic sex and reported sex) or individuals with sex chromosome aneuploidy were excluded from the analysis (n = 814).

Degree of relatedness
Estimated kinship coefficients using the KING toolset identified 107,162 pairs of related individuals.An inhouse algorithm was then applied to this list and preferentially removed the individuals related to the greatest number of other individuals until no related pairs remain.These individuals were excluded (n = 79,448).Additionally, 2 individuals were removed due to them being closely relating to a very large number (> 200) of individuals.

Outcomes and covariates
Our primary outcomes for this study were age and sex for BMI and Vitamin D exposures, and myocardial infarction when LDL-C was the exposure.We chose age and sex as key negative control outcomes as these cannot be caused by the exposures in question and any non-null estimates in any strata suggest the analyses are biased [39].However, there are well established selection effects relating to BMI, age, and sex into studies like UK Biobank (see e.g.[40]).It is plausible that these selection effects are amplified by nonlinear methods and could lead to bias.This still, of course, means the apparent findings are unreliable.
Myocardial infarction was extracted from the UK Biobank algorithmic definition (Data Showcase field 42000).MI outcomes were defined as prevalent (occurring before study enrolment), incident (occurring after study enrolment), or both (including both incident and prevalent cases).
We extracted covariates including UK Biobank recruitment centre, smoking status (coded as current or never/ ex) and the first 5 genetic principal components from UK Biobank directly.Smoking status was included for some sensitivity analyses given the known bidirectional association between smoking status and BMI [41].

Negative control outcome analysis -BMI and Vitamin D
To perform our negative control outcome analysis we first performed conventional MR then used the SUMnlmr package in R to perform both the residual and doubly-ranked method on our individual level data [42].Conventional MR was performed using a two stage residual inclusion approach in the OneSampleMR R package [43], using both the whole cohort and then in strata of age and sex.We first tested the effect of Vitamin D and BMI (individually) on age and sex, firstly without covariates (unadjusted), and then age (when sex was an outcome), sex (when age was an outcome), and also including UK Biobank recruitment centre, and the first genetic five principal components.Non-linear MR was performed using the standard settings in the SUMnlmr package, with a Gaussian distribution for linear outcomes, and a binomial distribution for binary outcomes.We chose to use ten strata for most analyses; but performed sensitivity analyses using different numbers of strata.
For BMI, we performed an analysis stratified by smoking status as this had been performed as a key analysis in an as yet uncorrected high-profile paper using non-linear MR [17] As BMI has a bidirectional association with smoking, this could induce collider bias, but we repeated these analyses to determine if estimates were more biased when stratifying by smoking status.

LDL cholesterol on myocardial infarction
In a further analysis, we focussed on the role of LDL-C as an exposure and myocardial infarction as an outcome.We provide the background to this positive control analysis in Box 1.
We performed two analyses.First, for completeness we replicated the negative control analysis discussed above (with age and sex as outcomes) using LDL-C as an exposure.Subsequently, we performed NLMR analyses on the effect of LDL-C on a) incident, b) prevalent, and c) all myocardial infarction to identify if effects matched the pattern shown by randomised trial data.
As non-linear MR results could be altered by statins, which have a large effect on LDL-C levels, we performed sensitivity analyses in a cohort including only statin users, only non-statin users, and those under the age of 50 (where statin use was rare, ~ 5%) and therefore not likely to bias estimates).Finally, we also fitted a model adjusting for statins, age, sex, and the first five genetic principal components.We recognise that the approach of analysing only those who use/do not use statins or adjusting for statins in our models generates collider bias [22] (as statin use is almost entirely downstream of LDL-C levels) but we included this for completeness.

Triglycerides and cancer mortality
Finally, given a recent paper identifying biologically highly implausible results, we aimed to replicate an analysis of the effect of triglycerides on cancer mortality reported in the above paper [29].We used the same definitions of cancer mortality and performed non-linear MR adjusting for age, sex, age * sex, age*age, age*age*sex, and the first ten genetic principal components.This choice of covariates was to match the prior analysis.We also -using the repeat sampling data available for around 15,000 participants -assessed the stability of strata for triglycerides, as these are known to be highly fluctuant.To do this we simply performed the doubly-ranked method at each time point (initial visit, repeat visit), and showed whether participants remained in the same strata using an alluvial plot.
In line with the authors of the doubly-ranked method recommendations [11] we performed multiple replicates (20) of each doubly-ranked approach and combined estimates using Rubin's rules.

Negative control outcomes as an assessment of the potential risk of bias
For our negative control outcome analysis, we generated instrumental variables (IVs) for two exposures (Vitamin D and BMI) using a summed linear score across single nucleotide polymorphisms (SNPs) and performed conventional and non-linear MR on participant age in the self-reported white British population in UK Biobank.
In total, we included 351,005 participants in our analysis using Vitamin D as an exposure and 383,793 participants in our analysis on BMI.Demographics of the population are in Supplementary Table S4.
We also calculated the effect of BMI and Vitamin D on age in males and females separately, and the effect of BMI and Vitamin D on sex in deciles of age.These are reported in Supplementary Table S5.Estimates across subgroups were similar-accepting the play of chance-across these analyses.
Box 1 Purpose of positive control analysis and background to the relationship between LDL-C and MI The purpose of positive control analyses is to perform an analysis whereby the true result is already known, and a comparison can be made to the results of a presented analysis and this truth.In this paper, we chose the effect of LDL-C on myocardial infarction.We chose this example because this is a well-studied exposure outcome relationship, where conventional MR estimates match in direction to randomised trial estimates, although MR estimates are around 40% larger, presumably due to the lifelong effect of lipid reduction conferred by genotype compared to that of drugs [1,44,45].In addition there are a large number of randomised trials with differing levels of baseline LDL-C [27].These individual trial estimates have been meta-regressed to identify if the effectiveness of LDL-C reduction depends on the baseline LDL-C [27,28], so we can calculate the shape of the LDL-C and MI relationship.In 2010, a meta-analysis of 26 trials (largely statin vs control) identified no clear evidence of the effect of baseline LDL-C on the effectiveness of more intensive (vs less intensive) LDL-C therapy (on all major vascular events), with a summary OR of 0.78 (95% CI 0.76-0.80)per 1 mmol/L reduction in LDL-C [28].In a more recent 2018 meta-analysis of 34 trials which reported specific LDL-C effects on myocardial infarction [27], the overall RR for more intensive vs less intensive therapy was similar (0.76, 95% CI 0.72-0.80),but in meta-regression they identified a small decrease in rate ratio (RR) of 0.90 (95% CI 0.84-0.97)with each 1 mmol/L increase in LDL-C (i.e. the effectiveness of more intense therapy was greater in those with higher LDL-C levels).In those with an LDL-C of more than 4 mmol/L, the RR for more vs less intense therapy was 0.64 (95% 0.53-0.78),while in those trials with a baseline LDL-C of less than 2.5 mmol/L, the RR was 0.85 (95% CI 0.76-0.92).The summed data therefore suggests a mostly linear effect with a moderate increase in effectiveness of more intense therapy in those with higher LDL-C It is possible for estimates from meta-regressions to be confounded, if for example those trials that included participants with higher baseline LDL-C, had also, on average, some other features that led to increased effectiveness of high intensity LDL-C lowering therapy (e.g.male sex, if male sex was associated with increased effectiveness of more intense therapy).However, the estimates from this meta-regression are highly biologically plausible, and consistent with our known understanding of cholesterol biology

Non-linear MR estimates differ substantially across strata
We then used both the residual and doubly-ranked method to generate stratum-specific estimates of the effect of each exposure on age and sex in univariable analyses.We chose to use ten strata, although results were similar using a differing number of strata (data not shown).
In contrast to the null estimate across the whole cohort, estimates across each stratum were markedly different using either method (Fig. 2).Focussing on the residual method, for the predicted effect of Vitamin D on age (Fig. 2A), there were positive effect estimates in lower strata which decreased moving from lower to upper strata for the effect of Vitamin D on age, leading to null effects.
However, the uppermost strata had a positive effect estimate.For BMI on age (Fig. 2C), we saw a similar trend, with positive effect estimates in lower strata and negative effect estimates in higher strata.
For the effect of BMI on sex (Fig. 2D), we saw positive estimates in lower strata and negative effects in upper strata, in line with a previous similar analysis we have performed, although that analyses used different parameters including a different instrumental variable [21].In short the estimates would suggest that BMI increased the odds for being male in those with the lowest non-genetic BMI and increased the odds of being female in upper strata.For the effect of Vitamin D on sex (Fig. 2B) we observed null effects in lower strata but modest positive effects in Fig. 2 The predicted causal effect of Vitamin D on age (A), sex (B) and BMI on age (C) and sex (D).Estimates were generated in each stratum using the residual method (blue) and the doubly-ranked method (red) and were unadjusted for covariates.The black estimate represents the conventional MR estimate generated using the whole cohort upper strata suggesting that Vitamin D increases the odds of being male in those with the highest non-genetic BMI.
For the doubly-ranked method, we identified largely null estimates for the effect of Vitamin D on age (Fig. 2A).For the effect of Vitamin D on sex (Fig. 2B), effect estimates were similar between the residual and doubly-ranked method, with increasingly positive effect estimates when moving from lower to upper strata.For the effect of BMI on age (Fig. 2C), we saw similar estimates to the residual method, with positive estimates in lower strata and negative estimates in upper strata.For the effect of BMI on sex (Fig. 2D), we saw marked non-null estimates with a similar trend but, in general, estimates closer to the null for the doubly-ranked relative to the residual method.We calculated Cochran's Q to formally assess the heterogeneity of strata specific estimates.These results are reported in Supplementary Figure S6.There was evidence-extremely strong in many cases-of strata specific estimate heterogeneity across all negative control outcome analyses except the analysis of Vitamin D on age and sex using the doubly-ranked method.
Estimates were similar in analyses adjusted for covariates: age (for sex as an outcome), sex (for age as an outcome) and the first 5 genetic principal components (Supplementary Figure S1).Log-transforming Vitamin D (as recommended by the creators of the residual method) brought estimates from the residual model closer on average to the doubly-ranked model (Supplementary Figure S2), although these were still non-null across many strata.
For our BMI analysis, we then stratified by smoking status (as was performed in a previous non-linear MR analysis in BMI, producing the headline findings from the paper [17]) and re-ran analyses (Supplementary Figure S3).Smoking had been shown to have a bidirectional relationship with BMI [41] before the non-linear BMJ MR paper was published [17], and therefore it should have been clear that this could have produced collider bias [22].In these analyses, stratification on smoking made considerable difference to some strata-specific estimates, although the general shape of the association remained similar.
To summarise, using both the residual and doubly-ranked method we identified non-null, stratum-specific associations between two exposures (Vitamin D and BMI) and two negative control outcomes (sex and age) across strata of the exposure in which the expected result is null.

LDL cholesterol and myocardial infarction
To examine both methods in a scenario where we anticipate the shape of the causal relationship we examined the association between LDL-C and a key outcome: myocardial infarction.RCT data from > 30 trials have demonstrated strong and broadly linear effects of LDL-C reduction on both outcomes with meta-analyses identifying slightly larger estimates of the effect of more intensive therapy in those with higher levels of LDL-C at baseline, while a recent MR analysis identified the opposite effect (increased effectiveness of LDL-C reduction in those with lower LDL-C [29], see Box 1 for further background).
First, as with analyses we presented above, we performed negative control outcome analyses of the effect of increasing LDL-C on age and sex.In conventional MR, estimates were essentially null for age (beta on age per 1 mmol/L increase in LDL -0.10; 95% CI -0.23, 0.03, p = 0.15), but suggested increased LDL-C was 'causal' for being more female (OR 0.95; 95% CI 0.92-0.98,p = 0.004).Conventional MR estimates of the effect of LDL-C on age were similar in men and women (Supplementary Figure S6), although there was some evidence of heterogeneity of MR estimates of LDL-C on sex across age strata, p value for heterogeneity 0.03).
Compared to conventional MR estimates, divergence from the null was much greater in non-linear MR (Fig. 3), with estimates on age as high as 5.74 (95% CI 1.79-9.63)years per 1 mmol/L LDL-C increase in one strata from the doubly-ranked method, a ~ 50 fold increase in bias compared to the conventional MR estimate.In contrast to our NCO analyses above, the bias was more extreme for the doublyranked method than the residual method.
We then went on to perform MR of the effect of increased LDL-C on myocardial infarction.For our primary analysis we included both incident and prevalent cases of MI, the effect on incident and prevalent MI separately are shown in Supplementary Figure S4.In conventional MR we saw the expected effect of LDL-C (OR for MI 1.74 per mmol/L increase in LDL again; 95% CI 1.61-1.87).However, in NLMR we saw unexpected effects, particularly with the doubly-ranked method.For the effect on MI (Fig. 4), we saw large differences in effects across strata, with the strongest effect in those with the lowest LDL-C, and the weakest effects in those in strata 5, 6, and 8.
For the residual method, the effect of LDL-C on MI was broadly stable and positive (although the trend still suggested reduction of the size of the effect across strata, p = 0.02).When running analyses adjusting for age, sex, and the first 5 principal components, effects were similar but less extreme for MI ( Supplementary Figure S5), although there was still a clear negative trend (p = 0.001 for using the doubly-ranked method), with effect estimates much larger in those with lower LDL-C than those with higher LDL-C.
We ran sensitivity analyses for LDL-C on MI that included a) adjusting for statin use as a covariate b) in statin users and non-statin users, and c) in under 50 s, where statin use was rare -5.1% ( Supplementary Figure S6).We recognise these estimates adjusting for statin use are highly likely to be biased due to collider bias but include these for interest.As expected, stratifying or adjusting for statin usage altered estimates dramatically.When adding statin as a covariate, estimates of LDL-C favoured protection in upper and lower strata, but were null in the middle strata (u-shaped).When analyses were performed using those only under 50, estimates were similar to our primary analyses but showed less precision.Similarly, analyses restricted to nonstatin users looked similar to our primary analyses.Analyses Fig. 3 The predicted causal effect of increased LDL-C on A age and B sex. Estimates were generated in each stratum using the residual method (blue) and the doubly-ranked method (red) and were unadjusted for covariates Fig. 4 The estimated causal effect of increased LDL-C on MI.Estimates were generated in each stratum using the residual method (blue) and the doubly-ranked method (red) and were unadjusted for covariates in statin users had reversed estimates in most strata, with increased LDL-C associated with reduced risk of MI.

Triglycerides and cancer mortality
In the recent paper on non-linear effects of lipids [29], one analysis focussed on the effect of triglycerides on cancer mortality.The overall effect was null in both univariable and multivariable MR, but the authors report extremely implausible results: a strong positive effect in strata 1 (of ten): an OR of 2.57 per mmol/L increase of TG (95% CI 1.67 to 3.96); but then a strong negative effect in strata 2: OR 0.56 (95% CI 0.39 to 0.83).All other strata specific estimates are close to the null.We report these estimates in Fig. 5.These results are especially implausible given the known variability in measurements of triglycerides [46,47].Therefore we aimed to replicate this as closely as possible, using the same dataset, exposure, outcome, and covariates.In our analysis (Fig. 6), we were unable to replicate this finding, despite the correlation between our strata specific mean triglyceride levels being > 0.99, suggesting we are analysing the same strata.
To further assess whether the prior analysis was unreliable, we used the repeat blood sampling data from UK Biobank, that was taken approximately 2 years after the original visit.Data on triglyceride levels at both time points was available for 13,535 participants.Correlation of TG measures between each time point was moderate (Pearson's R 0.60).As expected, when generating strata of TG using the doubly-ranked method on the original and repeat sample, participants were often not classified in the same strata.In fact, only 37% of those participants originally in strata 1 remain in strata 1, with 22% of them in strata 2, and the rest in higher strata.For those participants originally in strata 2, only 20% remain in strata 2, with 22% now in strata 1, and the rest in higher strata.These results are visualised in the alluvial plot below (Fig. 7).
This variability is not consistent with the reported estimates from Yang et al. [29] which would suggest those in strata 1 are at greatly increased risk of cancer mortality, and those in strata 2 at greatly reduced risk, with no effect elsewhere.If this were the case, more than half of the participants in the lower two strata dramatically change their risk of cancer mortality over 2 years, with some estimates flipping from greatly increased risk to protective.To present analyses which implicitly assumes this is possible is simply not credible.
Fig. 5 Estimates of the effect of triglycerides on cancer mortality from Yang et al. [29].Strata specific estimates generated using the doublyranked method [29]

Discussion
Mendelian randomization has become a commonly used approach in biomedical research [3], with the original intention being to strengthen causal inference regarding the effect of modifiable exposures on health outcomes [1].Quantified estimates of various average benefits that would be seen from treatment are often presented in MR papers [2], with the more recent development of methods that allowed for estimation of the effect of such treatments in groups at differing levels of the exposure [7][8][9][10].The residual method [10,7] is currently the most widely used approach.However recent findings have identified that this method can be seriously biased [19] leading to estimates that are inconsistent with causal interpretation.
In particular, the central estimate for some outcomes was precisely estimated and close to the null, while all strata specific estimates favoured protection.The authors of this paper and the residual method accept this is not possible and referred to it "as a logical impossibility if all estimates have a causal interpretation" [23] after the journal retracted the paper [24] following being alerted to it being obvious seriously problematic [19].However, whether the recently developed doubly-ranked method also leads to biased results remains unknown, leading us to perform our negative control outcome analyses.
In our negative control outcome analyses we identified that both methods showed bias, with effect estimates on age and sex differing across generated strata, when estimates should be null across all strata.It is worth considering some of the potential reasons for the distorted estimates in our negative control outcomes.One explanation is that these relate to a form of selection bias into strata which lead to associations between the genetic variants used to proxy for the exposure in the MR analyses-the so-called "instruments"-and other factors, including age and sex.This essentially reintroduces bias in MR analyses that were advanced with the intention of producing unbiased evidence on effects of exposures [1,2].
For body mass index overall selection bias (on entering studies, not on entering strata) is known to distort estimates, and thus the non-linear analysis of this exposure may be particularly liable to bias.For example, a large (~ 3.3 million person) GWAS of chromosomal sex identified over 150 autosomal loci associated with male sex, including the body mass increasing allele at FTO (OR 1.02, P = 4.4 × 10 −36 ) [40].Our Fig. 6 Estimates of the effect of triglycerides on cancer mortality using both the doubly-ranked (red) and residual method (blue).Estimates were generated adjusting for age, sex, age * sex, age * age, age * age * sex, and the first ten principle components findings may in part reflect the strong associations between age, sex, and body mass index in the selected population of UK Biobank.This would be supported by the weaker association with Vitamin D. However, this explanation is potentially less appealing for the effect of LDL-C-with the most extreme non-null estimates.Although there is likely selection onto LDL-C, there is little evidence that adjusting for participant bias in UK Biobank affects MR estimates of LDL-C [48].
Secondly, these NCO results could represent population stratification differing across strata, which can be uncovered using negative controls [26].Thirdly, these could represent methods related issues that occur on generation of the strata unrelated to selection.As we have demonstrated, stratification by both methods can lead to substantial alteration of estimates of IV-exposure associations.Regardless of the reason for the bias, these findings support more widespread use of negative control outcomes in Mendelian randomization studies, alongside a closer look at model assumptions.
Our analysis on the effect of LDL-C on myocardial infarction raises further concerns.The evidence base supporting the effectiveness of LDL-C reduction for MI is among the strongest in modern medicine, with > 30 randomised trials to date showing benefit [27].As such, meta-regressions of these trials using their baseline LDL-C have confidently suggested increased benefit of LDL-C reduction in those with higher LDL-C.This is also consistent with decades of non-RCT evidence, and our understanding of the pathogenesis of atherosclerosis (see Box 1 for further background).In non-linear MR (especially with the doubly-ranked method), we identified the opposite non-linear effect, suggesting a decreased effectiveness of LDL-C reduction in those with higher LDL-C.Importantly, this bias was not ameliorated by the addition of age, sex, and genetic principal components into the model.Estimates were also similar whether using incident, prevalent, or other MI, so these results are not likely to be generated by survivorship bias.
In fact, these results do not depend on our particular analytical approach and/or outcome definition.A recent paper from the authors of the doubly-ranked method performed a similar analysis on the same dataset [29] (using coronary artery disease as an outcome, not MI), and found similar findings, and with a strong negative trend with increasing LDL-C reduction (Fig. 8 compares both analyses).This consistent disagreement between trial data and MR data suggests a fundamental issue with the method, or with the application of this method to this data.This may reflect selection, but importantly, estimates were broadly similar in our primary analysis (unadjusted for age, sex and genetic principal components) and their analysis (which adjusts for the above), and suggests that if selection bias is occurring, it is not fixed by adjusting for age and sex.
It is therefore plausible these findings relate to other issues such as statin use (which has a large effect on LDL-C).However, adjusting for statins did not recover sensible estimates, nor would this be advised as this would lead to collider bias.Even if sensible estimates were recovered, this would suggest that non-linear MR estimates require understanding and modelling of the covariate structure of the exposure and the outcome, which would bring in all the issues of traditional observational analyses that MR is designed to remove [1,2].This particular example of LDL-C on CAD is particularly challenging for non-linear MR, as the strength of evidence supporting the effectiveness of conventional MR for the effect of LDL-C on CAD is very high, with multiple drugs having strong support from MR analyses [49].
Further implausible analyses have been generated in the same paper with an analysis of the predicted causal effect of triglycerides and cancer mortality.It is not possible to think of a plausible biological model where those in the bottom strata have very robustly estimated increased cancer mortality with higher triglycerides, whilst shifting into the 2nd lowest strata have robustly estimated decreased cancer mortality, with there being no causal effect for all other participants.In addition, given the considerable fluctuations demonstrated in triglyceride levels individuals would shift their classification in one into the other strata depending on exactly when the measurement for the particular individuals were made.Of note, we were unable to replicate the findings in this paper [29] despite being able to generate strata specific mean estimates of TG that exactly matched theirs, and using extremely similar definitions and analytical approaches, and despite being able to replicate similar findings for the effect of LDL-C on MI as they report on CAD.It is likely that researchers using these methods will identify similarly implausible results for an exposure-outcome relationship, and willbecause the data make little sense-selectively report only those exposure outcome relationships that make biological sense, leading to a publication bias of 'believable' non-linear MR papers, with a large number of performed but not reported unbelievable non-linear MR papers.This would lead to a body of published evidence where non-linear MR seems to 'work', but would not reflect the reality of numerous unpublished examples.
Therefore, our findings-and the results of other, independently performed analyses-suggest that we should be cautious when applying either the residual or doubly-ranked MR methods, as both methods seem potentially more susceptible to bias in the generation of strata specific estimates than conventional MR estimates are.Whether this entirely relates to selection bias or to other methodological issues remains unknown.Corrections should be considered for papers which have produced likely misleading estimates of non-linear effects [20,21,24], especially when their misleading findings gained widespread publicity due to the efforts of authors [50].
That these biases relate only to the examples of body mass index, Vitamin D, and lipid traits is unlikely, although it is of course challenging (and beyond the scope of this manuscript) to demonstrate why this method produces incorrect estimates.We urge authors of the relevant methods (and any future methods) to produce both further simulation work-and ideally empirical demonstrations of the reliability of any method before widespread use of it.

Limitations
This paper provides an initial assessment of the doubly doubly-ranked and residual methods for performing non-linear MR.In particular, our analysis focussed on three exposures (BMI, Vitamin D and lipids) a small number of outcomes and included only one dataset (UK Biobank), which we limited to White British participants.It is possible that our findings do not extend to all exposures, or only occur in certain settings.It may be that the methods perform better in other datasets, and particularly ones where selection bias is less prevalent than in UK Biobank.However, nearly all published non-linear analyses use UK Biobank, and it remains the largest easily accessible dataset with linked phenotypic and genetic data to perform these sort of analysis.Additionally, the authors of the doubly-ranked and residual method use worked examples from UK Biobank in the papers presenting the methods, suggesting this is an appropriate dataset in which to perform analyses.It would be useful to further examine these methods on diverse datasets, exposures, and outcomes, with proof of principle examples, to identify whether these issues are present in all datasets.

Conclusion
Using a negative control outcome approach, we identified bias in non-linear mendelian randomization estimates of the effect of BMI, Vitamin D and LDL-C using two commonly used methods.Estimates of the effect of LDL-C on cardiovascular outcomes were inconsistent with RCTs and biological understanding.Until reliable evidence is presented that the methods are generating sensible findings there should be a pause on further publication of non-linear MR findings using the two methods we examine in this paper.

Fig. 7 Fig. 8
Fig. 7 Classification of 13,535 participants who have repeat measurements into doubly-ranked strata based on triglyceride levels.The left hand Y axis represents the original strata, the right hand Y-axis the strata at repeat sampling ◂