# Statistical Methods for cis-Mendelian Randomization

Apostolos Gkatzionis1,2111Corresponding author. Email:
<EMAIL_ADDRESS>Stephen Burgess1,3, Paul J Newcombe1

(1MRC Biostatistics Unit, University of Cambridge, UK.
2MRC Integrative Epidemiology Unit, University of Bristol, UK.
3Department of Public Health and Primary Care, School of Clinical Medicine,
University of Cambridge, UK.)

###### Abstract

Mendelian randomization is the use of genetic variants to assess the existence
of a causal relationship between a risk factor and an outcome of interest. In
this paper we focus on Mendelian randomization analyses with many correlated
variants from a single gene region, and particularly on cis-Mendelian
randomization studies which uses protein expression as a risk factor. Such
studies must rely on a small, curated set of variants from the studied region;
using all variants in the region requires inverting an ill-conditioned genetic
correlation matrix and results in numerically unstable causal effect
estimates. We review methods for variable selection and causal effect
estimation in cis-Mendelian randomization, ranging from stepwise pruning and
conditional analysis to principal components analysis, factor analysis and
Bayesian variable selection. In a simulation study, we show that the various
methods have a comparable performance in analyses with large sample sizes and
strong genetic instruments. However, when weak instrument bias is suspected,
factor analysis and Bayesian variable selection produce more reliable
inference than simple pruning approaches, which are often used in practice. We
conclude by examining two case studies, assessing the effects of LDL-
cholesterol and serum testosterone on coronary heart disease risk using
variants in the HMGCR and SHBG gene regions respectively.

## Introduction

Mendelian randomization (MR) is the use of genetic data in order to assess the
existence of a causal relationship between a modifiable risk factor and an
outcome of interest (Davey Smith and Ebrahim, 2003; Burgess and Thompson,
2015). It is an application of instrumental variables analysis in the field of
genetic epidemiology, where genetic variants are used as instruments. The
approach has received much attention in recent years and has been used to
identify a large number of causal associations in the literature (Boef et al.,
2015). For example, MR studies have strengthened the evidence for a causal
link between lipoprotein(a) and coronary heart disease (CHD) (Kamstrup et al.,
2009), but have weakened the evidence for an effect of C-reactive protein on
CHD risk (CCGC (CRP CHD Genetics Collaboration), 2011).

A Mendelian randomization analysis requires a set of genetic variants that
satisfy the instrumental variable assumptions: genetic variants should be (i)
associated with the risk factor of interest, (ii) independent of confounders
of the risk factor-outcome association, and (iii) they should only influence
the outcome through their effect on the risk factor (the no-pleiotropy
assumption). Under these three assumptions, MR offers a framework for
assessing whether the risk factor is causally related to the outcome; under
additional modelling assumptions, one can also estimate the magnitude of the
risk factor-outcome causal effect.

XGYU✓xx(1)(2)(3) Figure 1: A causal diagram representation of the three
assumptions of Mendelian randomization. Here, $X$ represents the risk factor,
$Y$ the outcome, $G$ the genetic instrument and $U$ denotes confounders of the
$X-Y$ relationship.

Traditionally, Mendelian randomization studies are implemented using
independent genetic variants from across the whole genome. This is especially
the case when studying complex polygenic traits such as body mass index,
cholesterol or blood pressure. Such genome-wide analyses rely on published
results from large-scale GWAS studies to identify large numbers of genetic
regions associated with the risk factor studied. Genetic variants in each
region are pruned for independence and only the variant with the smallest
p-value, or a small number of weakly correlated variants with small p-values
are selected. Variants from different regions are then combined to create a
genome-wide set of instruments for MR, thus increasing the power of the
analysis to detect a causal relationship.

Conversely to the more traditional type of genome-wide MR described above,
recent years have seen an increased number of MR studies based on a single
gene region, and particularly cis-Mendelian randomization studies (Schmidt et
al., 2019; Gill et al., 2020). Cis-MR uses a protein as a risk factor, and
investigates the gene region encoding the studied protein in order to identify
variants that regulate protein expression. Mendelian randomization is then
conducted using these variants as instruments, in order to assess the effects
of protein expression on one or more disease outcomes. Cis-MR studies have
grown in popularity because they can be used for drug target discovery and
validation (Gill et al., 2020). Since proteins are the drug targets of most
medicines, an MR investigation of the effects of protein expression on a
disease can provide insights on whether a drug targeting the protein may
reduce the risk of suffering from the disease.

Cis-MR studies typically rely on a single gene region, and therefore have to
select genetic instruments among a potentially large number of correlated
variants. Through the existing cis-MR literature, a variety of techniques have
been employed to accomplish this task. However, to date no review or
comparison of the different methods has been published, and applied cis-MR
analyses typically use some form of LD-pruning for variable selection. The aim
of this paper is to provide an overview of statistical methods for cis-MR
studies with summary-level data. After reviewing commonly used approaches such
as LD-pruning and conditional analysis, we discuss the use of principal
components analysis (PCA), factor analysis and stochastic-search variable
selection. The use of PCA was proposed by Burgess et al. (2017) for MR with
correlated instruments and can be directly extended to cis-MR. The use of
factor analysis was recently proposed by Patel et al. (2020), building on
similar methods for instrumental variables analysis with individual-level data
(Bai and Ng, 2002). Stochastic search can be implemented using the JAM
algorithm, originally proposed by Newcombe et al. (2016) for fine-mapping
densely genotyped regions and recently adapted for MR by Gkatzionis et al.
(2020). We discuss how to implement each method in the context of cis-MR and
compare their performance in simulation studies.

We also provide two real-data applications, investigating associations of
variants in the SHBG and HMGCR regions with coronary heart disease risk. The
HMGCR gene encodes the protein HMG-coenzyme A reductase, which is known to
play an important role in the cholesterol-biosynthesis pathway. Its inhibition
with statins is known to reduce LDL-cholesterol levels, therefore statin
treatment is often prescribed to individuals in high CHD risk. The SHBG region
encodes the protein sex-hormone binding globulin, and is known to contain
genetic associations with testosterone and other sex hormones; there is some
uncertainty on whether testosterone levels are associated with CHD risk,
although a recent MR study has suggested no causal relationship (Schooling et
al., 2018).

## An Overview of Cis-Mendelian Randomization

Drug development is a challenging and costly process with uncertain results,
and the failure rate for drugs in clinical development can be as high as
$90\%$ (Hay et al., 2014). Since proteins are typically the targets of common
drugs, cis-Mendelian randomization uses expression data for druggable proteins
as risk factors, in order to aid in drug target discovery and validation. The
wealth of available genetic data, combined with the ability to routinely
implement MR analyses using existing software (Walker et al., 2019; Yavorska
and Burgess, 2017) makes cis-MR a useful approach to complement and inform
drug development. This can be evidenced by the increased success rate of drug
development programs with genomic support (Nelson et al., 2015).

Cis-MR studies can be used in order to assess the suitability of proteins as
potential drug targets, predict and inform the results of clinical trials, or
investigate repurposing opportunities and potential side-effects of existing
drugs. For example, a Mendelian randomization analysis of LDL cholesterol and
cardiovascular disease risk based on the ACLY region (Ference et al., 2019)
was recently used alongside a clinical trial (Ray et al., 2019) to assess the
suitability of the ACLY region as a potential drug target for lowering LDL-
cholesterol. Similar analyses have been conducted using the CETP region to
instrument lipid fractions and blood pressure (Sofat et al., 2010), using the
IL6R region to assess the effect of interleukin 6 reduction on coronary heart
disease risk (The Interleukin-6 Receptor Mendelian Randomisation Analysis
(IL6R MR) Consortium, 2012), and using the HMGCR (Swerdlow et al., 2015) and
PCSK9 regions (Schmidt et al., 2017) to explore potential effects of LDL-
cholesterol lowering treatments on body mass index and type 2 diabetes risk.

Ideally, cis-MR should be conducted using protein expression data for the
target protein risk factor (Gill et al., 2020). However, data on protein
expression are not always available, and researchers often use either gene
expression data for the associated gene region or genetic associations with a
downstream biomarker instead. For example, Swerdlow et al. (2015) used LDL-
cholesterol as a downstream risk factor in order to study the effects of HMGCR
inhibition on the risk of type 2 diabetes.

Genetic variants for cis-MR are selected from a region containing the protein-
encoding gene. The gene region is defined to contain the target gene, as well
as variants within a few hundred thousand base pairs on either side of the
gene, in order to include possible cis-eQTL variants. Within the limits of a
gene region, it is often possible to identify hundreds of correlated variants
marginally associated with protein expression. Researchers must then select
which of these variants to incorporate in a subsequent MR analysis; including
all available variants is usually avoided as it can introduce numerical
approximation errors in estimating the MR causal effect (Burgess et al.,
2017). Some authors have suggested prioritizing variants based on functional
annotations; however, Schmidt et al. (2019) investigated this strategy in a
simulation study and found little to no difference in causal inferences
conducted using functional variants compared to non-coding variants.
Therefore, the selection of genetic variants is usually performed based on the
strength of their associations with the risk factor, while accounting for LD
correlation to reduce numerical approximation errors. In the next section, we
review statistical approaches for implementing this variable selection
procedure.

Finally, similar to standard Mendelian randomization analyses, genetic
variants used for cis-MR must satisfy the instrumental variables assumptions.
Violations of the instrumental variable assumptions are often a concern
because they have the potential to devalidate tests of causal association
between the risk factor and the outcome and to bias causal effect estimates.
However, violations of the no-pleiotropy assumption are less likely to occur
when studying protein expression, because proteins are causally upstream of
metabolites and biomarkers, which constitute the most common risk factors used
in traditional MR analyses (Swerdlow et al., 2016; Schmidt et al., 2019).

## Statistical Methods for cis-Mendelian Randomization

### Notation

We now present various approaches for selecting genetic variants and computing
causal effect estimates in cis-MR studies. Let $X$ denote the risk factor of
interest, $Y$ denote the outcome, $G_{1},\dots,G_{P}$ denote a set of $P$
putative genetic instruments from a single gene region, and let $U$ represent
confounders of the risk factor-outcome association. Also, let $\theta$ denote
the causal effect of the risk factor on the outcome. Our task is to obtain a
subset of the genetic variants to use as instruments, an estimate
$\hat{\theta}$ of the causal effect parameter and an associated standard error
$\text{se}(\hat{\theta})$.

We work under the framework of two-sample summary-data MR and assume the
availability of two sets of summary statistics, one for SNP-risk factor and
one for SNP-outcome associations, obtained from different datasets. We use
$\hat{\beta}_{Xj}$, $\hat{\beta}_{Yj}$ to denote the marginal associations
between genetic variant $G_{j}$ and the two traits, and
$\hat{\sigma}_{Xj}^{2}$, $\hat{\sigma}_{Yj}^{2}$ for the corresponding
standard errors. In addition, $p_{Xj}$ denotes the p-value of association
between variant $G_{j}$ and the risk factor. Finally, we assume that we have
access to a reference dataset, such as $1000$ genomes or the UK Biobank, from
which to estimate correlations between genetic variants.

### Single-Variant MR

Some studies in the cis-MR literature have used only the variant with the
smallest univariate p-value in a region as an instrument (e.g. Sofat et al.,
2010; Swerdlow et al., 2015). With a single genetic variant, the MR causal
effect can be estimated using a ratio formula. This avoids the need to account
for genetic correlations and to perform complex numerical tasks such as
inverting the genetic covariance matrix. On the other hand, if the region
contains multiple variants with independent associations with the risk factor,
a single-variant MR analysis will not be able to model all the genetic effects
in the region. This typically translates into a loss of power for the cis-MR
analysis. Therefore, most recent cis-MR papers have utilised variable
selection or dimension reduction approaches in order to include several
variants per region.

### LD-Pruning

Perhaps the most common approach in the literature for selecting genetic
variants for inclusion into a cis-MR study is LD-pruning (see e.g. Dudbridge
and Newcombe, 2015; Schmidt et al., 2019). Given a set of candidate genetic
variants $G_{1},\dots,G_{P}$, the stepwise pruning procedure iterates between:

  * •

Selecting the variant with the smallest p-value and including it in the
analysis.

  * •

Identifying all variants whose correlation with the previously selected
variant exceeds a correlation threshold $\rho$ and removing them from further
consideration.

This is repeated until all the candidate variants have either been selected
for inclusion or discarded. Alternatively, the process may be stopped when
there are no remaining variants with p-values lower than a significance
threshold $\tau$, in order to guard against the inclusion of weak instruments
into the analysis. The threshold $\tau$ is often taken to be the GWAS
significance threshold $\tau=5\times 10^{-8}$, although there is some evidence
that using more relaxed thresholds may be beneficial (Dudbridge, 2013).

Once LD-pruning is implemented and a suitable set of genetic instruments is
selected, the MR causal effect can be estimated using an inverse variance
weighted (IVW) estimate. If the correlation threshold $\rho$ is small enough,
the selected variants can be considered approximately independent and the
standard IVW estimate,

$\displaystyle\hat{\theta}_{ivw}$ $\displaystyle=$
$\displaystyle\frac{\sum_{j=1}^{P}\hat{\beta}_{Xj}\hat{\beta}_{Yj}\hat{\sigma}_{Yj}^{-2}}{\sum_{j=1}^{P}\hat{\beta}_{Xj}^{2}\hat{\sigma}_{Yj}^{-2}}$
$\displaystyle\text{se}(\hat{\theta}_{ivw})^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{\sum_{j=1}^{P}\hat{\beta}_{Xj}^{2}\hat{\sigma}_{Yj}^{-2}}$

can be used. This approach is typically used in genome-wide MR analyses.
However, for single-region MR, the use of a small correlation threshold may
result in the exclusion of causal variants, especially when studying gene
regions with multiple independent associations with the risk factor. On the
other hand, if a larger threshold is used, the standard IVW estimate will
underestimate the causal standard error, since it does not adjust for
correlations. In this case it is better to use a modified version of of the
IVW estimator that takes genetic correlations into account (Burgess et al.,
2016):

$\displaystyle\hat{\theta}_{cor}$ $\displaystyle=$
$\displaystyle\left(\hat{\beta}_{X}^{T}\Omega^{-1}\hat{\beta}_{X}\right)^{-1}\left(\hat{\beta}_{X}^{T}\Omega^{-1}\hat{\beta}_{Y}\right)$
(1) $\displaystyle\text{se}(\hat{\theta}_{cor})^{2}$ $\displaystyle=$
$\displaystyle\left(\hat{\beta}_{X}^{T}\Omega^{-1}\hat{\beta}_{X}\right)^{-1}$
(2)

where $\Omega$ is (an estimate of) the genetic covariance matrix:
$\Omega_{i,j}=\text{Cov}(G_{i},G_{j})$. The estimator (1) can be motivated
from the literature on meta-analysis. It is sometimes referred to as the
generalized least squares method (Schmidt et al., 2019) because it can be
derived using generalized linear regression to model the genetic associations
with the outcome in terms of the genetic associations with the risk factor
(Burgess et al., 2016; Burgess et al., 2017).

The correlated-instruments IVW estimator requires the genetic covariance
matrix $\Omega$ to be computed and inverted, a process which can become
numerically unstable if highly correlated variants are used (Burgess et al.,
2017). This is the main reason why cis-MR analyses cannot use all available
variants within a region. Numerical instabilities are less likely to occur
when the IVW estimator is used in combination with pruning, since the pairwise
correlations of variants selected by LD-pruning are all lower than the
specified correlation threshold $\rho$, although problems can still occur if
the genetic covariance matrix is near-singular.

A potential criticism for the LD-pruning approach is that there is no
consensus in the literature on how to select the correlation threshold.
Implementations of the method using different correlation thresholds can give
rise to different results, especially when large thresholds are used (Burgess
et al., 2017). To address this issue, it has been proposed that stepwise
pruning is implemented with a variety of pruning thresholds, as a form of
sensitivity analysis (Schmidt et al., 2019). In addition, stepwise methods for
variable selection have been criticised for being unstable and getting stuck
in local maxima of the model space (Hocking, 1976; Tibshirani, 1996; Miller,
2002) and more flexible variable selection approaches are often preferred in
some applications of biostatistics (Newcombe et al., 2016; Benner et al.,
2016; Asimit et al., 2019), although it is not clear whether the added
flexibility provides a substantial benefit in Mendelian randomization studies.

### Conditional Analysis

Conditional analysis is a process similar to LD-pruning, except variant
selection is based on genetic associations with the risk factor conditional on
any previously selected variants instead of univariate marginal p-values.
Perhaps the most established algorithm for conditional analysis is the
Conditional and Joint (CoJo) algorithm proposed in Yang et al. (2012). The
algorithm starts by selecting the genetic variant with the smallest univariate
p-value, and then in each iteration it iterates between the following steps:

  * •

Calculate the p-values of all candidate genetic variants conditional on any
variants that have already been selected in previous steps. As with stepwise
pruning, the p-values are only computed for genetic variants whose
correlations with previously selected variants do not exceed a pre-specified
correlation threshold $\rho$.

  * •

Select the variant with the smallest conditional p-value, provided that its
p-value is smaller than a significance threshold $\tau$.

  * •

Regress the risk factor on all previously selected variants and compute the
multivariate p-values for each selected variant.

  * •

If any of the p-values in the multivariate model is larger than the
significance threshold $\tau$, drop the variant with the largest p-value from
the model.

This variable selection algorithm is run until it can no longer add or remove
any variants. Upon finishing, it provides a set of genetic variants which can
subsequently be used to conduct Mendelian randomization. A causal effect
estimate based on the selected variants is then obtained using correlated-
instruments IVW (1)-(2).

Yang et al. (2012) provide details on how to compute the conditional and
joint-model p-values using marginal summary statistics. The fact that the
algorithm only requires summary-level data, along with its implementation as
part of the GCTA software, have contributed to conditional analysis becoming a
popular option for selecting variants from single genetic regions.

In practice, conditional analysis often has a similar performance to LD-
pruning (Dudbridge and Newcombe, 2015). Both approaches rely on a correlation
threshold parameter $\rho$ and can produce inconsistent results for different
values of the correlation parameter, especially for large values (e.g.
$\rho\geq 0.9$), due to numerical instabilities.

### Principal Component Analysis

Principal components analysis (PCA) is widely used in GWAS studies to adjust
for population stratification. In the context of correlated-instruments
summary-data MR, the use of principal components was proposed by Burgess et
al. (2017). The method aims to identify linear combinations of genetic
variants that are orthogonal to each other and explain as much of the variance
in the genetic dataset as possible. These linear combinations are then used as
genetic instruments to assess the causal relationship between the risk factor
and outcome, using an IVW estimate.

In Burgess et al. (2017), principal components analysis is applied to a
transformed version $\Psi$ of the covariance matrix, with elements

$\Psi_{ij}=\hat{\beta}_{Xi}\hat{\beta}_{Xj}\hat{\sigma}_{Yi}^{-1}\hat{\sigma}_{Yj}^{-1}\rho_{ij}=\hat{\beta}_{Xi}\hat{\beta}_{Xj}\hat{\sigma}_{Yi}^{-2}\hat{\sigma}_{Yj}^{-2}\Omega_{ij}$

where $\rho_{ij}=\text{Cor}(G_{i},G_{j})$. The matrix $\Psi$ is preferred over
the untransformed matrix of genetic correlations because its entries depend on
the precisions $\hat{\sigma}_{Yj}^{-1}$ of univariate estimates and thus, if
two genetic variants are almost perfectly correlated, the variant with the
highest precision will be prioritized.

Principal components analysis uses an eigendecomposition for the matrix $\Psi$
and identifies the diagonal matrix of eigenvalues $\Lambda=(\lambda_{ii})$ and
the corresponding matrix of eigenvectors $W$ so that $\Psi=W^{T}\Lambda W$. By
transforming the summary statistics

$\bar{\beta}_{X}=W^{T}\hat{\beta}_{X}\;\;\;,\;\;\;\bar{\beta}_{Y}=W^{T}\hat{\beta}_{Y}$

and the genetic covariance matrix

$\bar{\Omega}=W^{T}\Omega W$

one can rewrite the IVW estimate (1) in terms of these transformed vectors:

$\hat{\theta}_{cor}=\left(\bar{\beta}_{X}^{T}\bar{\Omega}^{-1}\bar{\beta}_{X}\right)^{-1}\left(\bar{\beta}_{X}^{T}\bar{\Omega}^{-1}\bar{\beta}_{Y}\right)$
(3)

Evaluating (3) is subject to the same numerical difficulties as (1). However,
one can construct a numerically stable approximation for $\hat{\theta}_{cor}$
by using only the first $k$ principal components. If $\tilde{W}$ is the sub-
matrix of $W$ containing only the first $k$ columns and
$\tilde{\beta}_{X}=\tilde{W}^{T}\hat{\beta}_{X}$,
$\tilde{\beta}_{Y}=\tilde{W}^{T}\hat{\beta}_{Y}$,
$\tilde{\Omega}=\tilde{W}^{T}\Omega\tilde{W}$, then the causal effect estimate
based on the first $k$ principal components is

$\hat{\theta}_{cor}\approx\left(\tilde{\beta}_{X}^{T}\tilde{\Omega}^{-1}\tilde{\beta}_{X}\right)^{-1}\left(\tilde{\beta}_{X}^{T}\tilde{\Omega}^{-1}\tilde{\beta}_{Y}\right)$

The number $k$ of principal components to use is a tuning parameter for the
method. In practice, it is often specified so that the selected principal
components explain a specific proportion of variation in the genetic dataset.
The criterion used in Burgess et al. (2017) is to select principal components
that explain either $99\%$ or $99.9\%$ of the variation in the data.

### The JAM Algorithm

Joint analysis of marginal summary statistics (JAM) is a Bayesian stochastic-
search variable selection algorithm that was introduced by Newcombe et al.
(2016) for the purpose of statistical fine-mapping from summary GWAS results.
The algorithm aims to identify genetic variants robustly and independently
associated with a trait of interest, among a large set of correlated candidate
variants. JAM only requires the availability of GWAS summary-level data, and
genetic correlations which can be estimated from an external reference
dataset. Using these summary data, JAM identifies sets of genetic variants
that are most plausibly associated with the trait of interest.

The algorithm uses linear regression to model trait $X$ in terms of genetic
variants $G=(G_{1},\dots,G_{P})$:

$X\sim N\left(Gb,\sigma_{X}^{2}\right)$ (4)

where $b$ denotes the vector of multivariate associations between each genetic
variant and the trait (in contrast to $\beta_{X},\beta_{Y}$, which represent
marginal associations). In order to implement this regression using summary
statistics, JAM employs a transformation $z=G^{T}X$, for which (4) implies
that

$z\sim N\left(G^{T}Gb,G^{T}G\sigma_{X}^{2}\right)$ (5)

Assuming Hardy-Weinberg equilibrium and that the trait measurements and
genetic data have been centered, each element $z_{j}$ of the vector $z$ can be
approximated using the marginal summary statistics $\hat{\beta}_{Xj}$ and
effect allele frequencies $f_{j}$,

$z_{j}=2f_{j}(1-f_{j})N_{1}\hat{\beta}_{Xj}$

where $N_{1}$ is the sample size from which the G-X associations are inferred.
In practice, an additional Cholesky decomposition is implemented to avoid
modelling under the correlated error structure of (5); we refer to Newcombe et
al. (2016) for the technical details. In addition, note that $G^{T}G$ is
($N_{1}-1$ times) the genetic covariance matrix, which can be estimated from
reference data. Hence, equation (5) can be used to construct a summary-data
likelihood $\ell(z,G_{ref}|b,\sigma_{X}^{2})$. A similar likelihood is used to
conduct the joint analysis step in Yang et al. (2012). JAM employs Bayesian
analysis instead, and uses conjugate normal-inverse-gamma g-priors
$p(b,\sigma_{X}^{2})$ to facilitate Bayesian inference for the genetic
associations with the trait.

If the trait studied is binary, $z$ can be derived by mapping univariate log-
odds ratios to the univariate effects that would have been estimated if the
binary trait was modelled by linear regression, as has been done in other
linear-based summary data algorithms (Vilhjalmsson et al., 2015).

The likelihood (5) relies on a fixed set of genetic variants. For variable
selection, we assume that the likelihood has been conditioned on a specific
set of variants $\gamma$ and use a Beta-Binomial prior $p(\gamma)$ for the
model coefficient in order to obtain the posterior

$p(\gamma,b,\sigma_{X}^{2}|z,G_{ref})\propto\ell(z,G_{ref}|b,\sigma_{X}^{2},\gamma)p(b,\sigma_{X}^{2}|\gamma)p(\gamma)$
(6)

This posterior is difficult to evaluate analytically, but JAM implements
stochastic-search model selection via a reversible-jump Markov Chain Monte
Carlo (RJMCMC) algorithm to sample from (6) and identify the most suitable
sets of genetic variants. The stochastic search procedure starts from a set
$\gamma_{(0)}$ containing only the variant with the smallest p-value. Then in
each iteration, given a current set of genetic variants $\gamma_{(k)}$, JAM
randomly proposes a new set $\gamma_{(k+1)}$ by either

  * •

Adding a new variant to the current set $\gamma_{(k)}$.

  * •

Removing a variant from $\gamma_{(k)}$.

  * •

Swapping a variant in $\gamma_{(k)}$ for a variant not in $\gamma_{(k)}$.

The algorithm then decides whether to accept the new set by evaluating how
well it fits the trait, according to the posterior (6). The process is
repeated for a large number of iterations, and the various sets of variants
are assigned posterior probabilities according to how often they were visited.
The stochastic search procedure allows JAM to evaluate a wide range of
different models and efficiently explore the parameter space of genetic
configurations.

For Mendelian randomization, the above procedure can be used to identify
variants strongly and independently associated with the risk factor of
interest (Gkatzionis et al., 2020). An IVW causal effect estimate can then be
obtained for each set of variants visited by JAM, and these model-specific
estimates can be combined into an aggregate causal effect estimate by model
averaging.

Since it jointly models all available genetic variants and can account for
genetic correlations, JAM is naturally suited for cis-MR. The algorithm
discards variants which are associated with the risk factor only through their
correlation with other present variants, and its posterior models include
variants independently associated with the risk factor. JAM is not completely
free of numerical issues when implemented on near-collinear variants, but
these issues can be overcome by employing pruning as a first step before
calling the algorithm. An advantage over the LD-pruning approach is that JAM
is more robust to the choice of pruning thresholds, since it implements a
second layer of variable selection after pruning is finished. The pruning step
can be implemented using a high correlation threshold (e.g. $\rho=0.9$) and no
significance threshold. As an alternative, JAM can avoid the need for pruning
by using a ridge term for the genetic covariance matrix to make its inversion
more stable. The choice of a pruning threshold is then replaced by the tuning
of the ridge term parameter. Here we focus on implementations using pruning as
a first step and do not explore the use of a ridge term.

Outside of the context of cis-MR, Gkatzionis et al. (2020) extended JAM to
handle violations of the no-pleiotropy assumption when working with multiple
independent genetic variants, by augmenting its variable selection with a
heterogeneity loss function to penalize and downweight pleiotropic variants.
Furthermore, Jiang et al. (2020) proposed a hierarchical version of the
algorithm that is useful for multivariable MR, as well as transcriptome-wide
association studies.

### Factor-based Methods

Recently, Patel et al. (2020) proposed the use of factor analysis for MR with
correlated weak instruments. The authors model the genetic variants in terms
of a factor model,

$G=\Lambda f+\epsilon_{f}$ (7)

where $f$ is a vector of $k<P$ latent factors, $\Lambda$ is a matrix of factor
loadings and $\epsilon_{f}$ is a vector of idiosyncratic errors. The matrix
$\Lambda$ can be estimated using the first $k$ eigenvectors in the
eigendecomposition of the (rescaled) genetic covariance matrix $G^{T}G$. This
covariance matrix and its eigendecomposition can be approximated if one has
access to a reference dataset. The number $k$ of latent factors can be
selected either empirically, as is common in factor analysis, or in a more
disciplined way, for example using the factor penalization method of Bai and
Ng (2002).

The factors are then used as genetic instruments in a standard MR model:

$\displaystyle X$ $\displaystyle=$ $\displaystyle b_{X}^{T}f+\epsilon_{X}$
$\displaystyle Y$ $\displaystyle=$ $\displaystyle b_{Y}^{T}f+\epsilon_{Y}$
$\displaystyle\beta_{Y}$ $\displaystyle=$ $\displaystyle\theta\beta_{X}$

where $\epsilon_{X},\epsilon_{Y}$ are mean-zero error terms. An estimate of
the causal effect parameter $\theta$ can be obtained by limited-information
maximum likelihood (LIML), minimizing the criterion

$\hat{g}(\theta)=\frac{1}{P}\hat{\Lambda}^{T}\hat{D}_{Y}\hat{\beta}_{Y}-\theta\frac{1}{P}\hat{\Lambda}^{T}\hat{D}_{X}\hat{\beta}_{X}$

where $\hat{\Lambda}$ is the estimated matrix of factor loadings and
$\hat{D}_{X},\hat{D}_{Y}$ are the diagonal matrices whose diagonal terms are
the sample variances of each genetic instrument; these can be approximated
from the reference dataset. Minimization of (a weighted version of)
$\hat{g}(\theta)$ yields the factor LIML (F-LIML) estimate.

Patel et al. (2020) describe conditions that guarantee the consistency and
asymptotic normality of the F-LIML estimator. Unfortunately, these conditions
are often violated under a weak instruments framework and F-LIML may suffer
from weak instrument bias. In this scenario, Patel et al. (2020) propose a
range of hypothesis tests, using the factors as genetic instruments, in order
to assess the existence of a causal effect. These tests are derived from the
instrumental variables literature and adapted to work with summary-level data.
They include the Anderson-Rubin (Anderson and Rubin, 1949), Lagrange
multiplier (Kleibergen, 2002) and conditional likelihood ratio (Moreira, 2003)
tests. The various tests are more robust to weak instrument bias than the
F-LIML estimator, and have higher power than the principal components method
under a weak instruments setting.

## Simulation Study

### Simulation Design

To compare the various cis-MR approaches, we implemented a simulation study.
We generated data under the following simulation model:

$\displaystyle X$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{P}b_{Xj}G_{j}+\alpha_{X}U+\epsilon_{X}$ (8)
$\displaystyle Y$ $\displaystyle=$ $\displaystyle\theta
X+\alpha_{Y}U+\epsilon_{Y}$ (9) $\displaystyle U,\epsilon_{X},\epsilon_{Y}$
$\displaystyle\sim$ $\displaystyle
N(0,\sigma_{0}^{2})\;\;\;\;\text{independently of each other.}$ (10)

Our simulation design was informed by the real-data applications in the next
section. We obtained genetic variants from two gene regions: the SHBG region,
which encodes the protein sex-hormone binding globulin and is known to be
associated with sex hormone traits, and the HMGCR region, which encodes the
protein HMG coenzyme-A reductase, plays an important role in the metabolic
pathway that produces cholesterol and is the drug target of statins. We refer
to the real-data applications for more information about these regions. SHBG
and HMGCR represent two distinctly different genetic correlation structures;
correlation heat-maps for the two regions are provided in Figures 2 and 3
below. By basing simulations on two different regions we guarded against
spurious observations on comparative performance of the methods that may arise
due to subtle properties of the correlation structure in any one region.

We defined the two gene regions as spanning 100Kb pairs on either side of the
corresponding gene and used all variants in these regions for which
individual-level data were available in the UK Biobank and the minor allele
frequencies were higher than $1\%$. This resulted in $P_{1}=717$ genetic
variants for the SHBG region and $P_{2}=590$ variants for the HMGCR region. We
extracted individual-level genetic observations for the selected variants from
the UK Biobank, using a sample of $N_{ref}=367643$ non-related individuals of
European origin. Aiming to replicate a two-sample Mendelian randomization
analysis, we then bootstrapped the UK Biobank data and obtained two genetic
matrices $G_{1},G_{2}$ of sizes $N_{1},N_{2}$ respectively. Using the genetic
matrix $G_{1}$, we simulated risk factor measurements according to (8) and
used them to obtain $G-X$ summary statistics. Using the matrix $G_{2}$ and
(8)-(9), we simulated outcome measurements and used them to generate $G-Y$
summary statistics. Finally, the reference correlation matrix was computed
using the entire UK Biobank dataset.

To generate risk factor measurements, we used information from our real-data
applications and the literature. For the SHBG region, we assumed that six
independent genetic effects were present in the region. This was inspired by
Coviello et al. (2012) who investigated the associations of SHBG variants with
concentration of the SHBG protein, and suggested that as many as nine
independent signals may be present in the region. Three of the nine causal
variants identified by Coviello et al. (2012) were not included in our
dataset, as they were located more than 100Kb pairs away from the SHBG gene,
therefore we used the remaining six variants in our simulations. For these six
variants, the effects of the risk-increasing allele on the risk factor were
generated according to $b_{Xj}\sim|N(0,0.2)|+0.1$. This created a pattern of
univariate SNP-exposure associations similar to that we observed in our real-
data application.

For the HMGCR region, we also assumed the existence of six variants
independently associated with the risk factor. The six causal signals were
placed at genetic variants used by Ference et al. (2019) to construct a
genetic score for LDL-cholesterol based on the HMGCR region. Effects were
generated according to $b_{Xj}\sim|N(0,0.03)|+0.03$ for the risk-increasing
allele; again, this closely resembled the univariate summary statistics
obtained in our real-data application.

We considered three different values for the causal effect parameter:
$\theta=0,0.05,0.1$. We set the $G-X$ sample size to be $N_{1}=10000$, which
is fairly small for modern MR standards but may be reasonable for a cis-MR
study using protein expression data for the risk factor. In supplementary
material, we report simulation results using a larger sample size of
$N_{1}=100000$. The sample size for SNP-outcome associations was set equal to
$N_{2}=180000$, similar to that of the CARDIoGRAMplusC4D dataset
(CARDIoGRAMplusC4D Consortium, 2015) used in our real-data applications.
Finally, we assumed that all the genetic variants are valid instruments and
did not generate pleiotropic effects.

One of our goals in the simulations was to assess the effect of weak
instruments bias on the performance of the various methods. This form of bias
can be of particular concern in cis-MR studies, as focusing on a single region
means that usually there will be far fewer genetic instruments to use. A
common diagnostic for weak instrument bias is to compute the F statistic for
the regression of the risk factor on all the genetic variants; a rule-of-thumb
is that $F\leq 10$ is indicative of bias.

Therefore, we considered two simulation scenarios. In the first scenario, the
proportion of variation in the risk factor explained by the genetic variants
in each region was specified according to what has been empirically observed
for their protein products, and was set equal to $3\%$ for the SHBG region
(Jin et al., 2012) and $2\%$ for the HMGCR region (Krauss et al., 2008). This
is a “strong instruments” setting, in which an “oracle” MR analysis using only
the six causal variants in each region attained an average F statistic of $52$
for the SHBG region and $35$ for the HMGCR region. In this scenario, the
effect of weak instrument bias on the performance of cis-MR methods should be
small.

In our second scenario, we reduced the proportion of variation in the risk
factor that is explained by the genetic variants to $10\%$ of its value in the
previous scenario. For the SHBG region, we assumed that the genetic variants
only explain $v_{G}=0.3\%$ of variation in the risk factor, while for the
HMGCR region, the genetic variants explained $0.2\%$ of variation. This
resulted in a “weak instruments” scenario in which the “oracle” MR analysis
had an average F statistic of $6.0$ for the SHBG region and $4.2$ for the
HMGCR region.

For each scenario, the simulation was replicated $1000$ times per region and
per value of the causal effect.

### Methods

In our simulations we implemented the following cis-MR methods:

  * •

Single-instrument MR using only the minimum p-value variant in the region.

  * •

LD-pruning.

  * •

Principal components analysis.

  * •

The JAM algorithm.

  * •

The F-LIML algorithm.

  * •

The factor-based conditional likelihood ratio (CLR) test.

The performance of the conditional method is often quite similar to that of
stepwise pruning (Dudbridge and Newcombe, 2015), hence the method was not
implemented. A brief comparison between conditional analysis, stepwise pruning
and PCA can be found in Burgess et al. (2017).

Many of the above methods depend on additional tuning parameters. To assess
the sensitivity of each method to the specification of its tuning
parameter(s), we used a range of different values. For stepwise pruning, we
set a correlation threshold of $\rho=0.1,0.3,0.5,0.7,0.9$. For the
significance threshold, we use the standard GWAS threshold ($\tau=5\times
10^{-8}$) in simulations with “strong instruments”; with “weak instruments”
there were simulations in which no variants passed the GWAS significance
threshold, so we used a lower value ($\tau=10^{-3}$) instead. For the PCA
method, we used principal components that explained $k=0.99$ or $k=0.999$ of
the variation in the genetic data. For JAM, we implemented a pre-pruning step
with a correlation threshold of $\rho=0.6,0.8,0.9,0.95$ and no significance
threshold - note that JAM does not need a significance threshold for the pre-
pruning step because its variable selection discards variants weakly
associated with the risk factor anyway. The algorithm was run for 1 million
iterations in each instance. Finally for the F-LIML method, we allowed the
algorithm to determine the number of latent factors to use.

The simulations were implemented in the statistical software R. Stepwise
pruning and correlated-instruments IVW were implemented manually. To implement
the PCA approach, we used code provided in the appendix of Burgess et al.
(2017). The JAM algorithm was implemented using the R package R2BGLiMS. The
F-LIML algorithm and the summary-data CLR test were implemented using code
provided to us by the authors of the relevant publication (Patel et al.,
2020). Stepwise pruning was the fastest algorithm to implement and the JAM
algorithm was the slowest, although the differences were small on an absolute
scale and none of the methods required more than a few seconds to run for each
dataset.

The various MR methods in our simulations were subject to two sources of bias.
The first source is weak instrument bias, which is known to act towards the
direction of the observational association in one-sample MR analyses and
towards the null in two-sample MR (Burgess et al., 2011; Burgess et al.,
2016). As discussed earlier, weak instrument bias should be a concern in our
second simulation scenario, where the F statistic is below $10$. Note however
that even in the first scenario, weak instrument bias may still affect methods
using many nuisance variants to obtain a causal effect estimate, as using such
variants reduces the value of the F statistic.

Second, the performance of various methods in our simulations can be affected
by numerical issues, in particular related to inverting the genetic covariance
matrix $G^{T}G$. Inaccurate estimation of $(G^{T}G)^{-1}$ may lead to spurious
increases or decreases in genetic correlations. This form of bias can be
detected by computing the condition number of the genetic covariance matrix,
but its direction is not straightforward to assess. It is more likely to
affect methods using highly correlated variants, such as stepwise pruning with
high correlation thresholds.

### Results

Simulation results are reported in Tables 1-2. For each method and each value
of the causal effect parameter, we report median causal effect estimates and
estimated standard errors. For simulations with $\theta=0$ we also report the
empirical Type I error rates, defined as the proportion of simulations in
which the method rejected the null causal hypothesis. For simulations with
$\theta\neq 0$ we also report the empirical coverage of confidence intervals
and the empirical power to reject the null causal hypothesis. All methods had
very high power in the “strong instruments” scenario with $\theta=0.1$, and
very low power in the “weak instruments” scenario with $\theta=0.05$,
therefore we do not report power in these two scenarios. Finally, the
conditional likelihood ratio test does not provide a causal effect estimate or
an associated standard error, hence we only report its coverage, power and
Type I error rate.

|  | $\theta=0$ | $\theta=0.05$ | $\theta=0.1$
---|---|---|---|---
Method | $\hat{\theta}$ | $se(\hat{\theta})$ | Type I | $\hat{\theta}$ | $se(\hat{\theta})$ | Cov | Power | $\hat{\theta}$ | $se(\hat{\theta})$ | Cov
SHBG Region
Top SNP | —– | 0.000 | 0.019 | 0.051 | 0.050 | 0.019 | 0.956 | 0.731 | 0.096 | 0.020 | 0.927
Pruning | $\rho=0.1$ | -0.002 | 0.016 | 0.051 | 0.048 | 0.017 | 0.938 | 0.837 | 0.096 | 0.017 | 0.936
$\rho=0.3$ | -0.001 | 0.014 | 0.049 | 0.048 | 0.014 | 0.949 | 0.923 | 0.096 | 0.015 | 0.932
$\rho=0.5$ | -0.001 | 0.013 | 0.051 | 0.047 | 0.014 | 0.943 | 0.941 | 0.095 | 0.014 | 0.927
$\rho=0.7$ | -0.001 | 0.013 | 0.050 | 0.046 | 0.013 | 0.939 | 0.940 | 0.092 | 0.014 | 0.901
$\rho=0.9$ | 0.000 | 0.012 | 0.047 | 0.043 | 0.013 | 0.912 | 0.928 | 0.086 | 0.013 | 0.796
PCA | $k=0.99$ | -0.001 | 0.015 | 0.063 | 0.049 | 0.015 | 0.954 | 0.914 | 0.099 | 0.016 | 0.938
$k=0.999$ | 0.000 | 0.014 | 0.050 | 0.048 | 0.014 | 0.941 | 0.933 | 0.096 | 0.015 | 0.935
JAM | $\rho=0.6$ | -0.001 | 0.014 | 0.048 | 0.048 | 0.014 | 0.954 | 0.935 | 0.095 | 0.015 | 0.938
$\rho=0.8$ | 0.000 | 0.014 | 0.049 | 0.047 | 0.014 | 0.957 | 0.934 | 0.095 | 0.015 | 0.939
$\rho=0.9$ | 0.000 | 0.014 | 0.041 | 0.048 | 0.014 | 0.955 | 0.932 | 0.094 | 0.015 | 0.935
$\rho=0.95$ | 0.000 | 0.014 | 0.045 | 0.047 | 0.014 | 0.956 | 0.933 | 0.094 | 0.015 | 0.933
F-LIML | —– | 0.000 | 0.014 | 0.066 | 0.051 | 0.015 | 0.942 | 0.936 | 0.100 | 0.016 | 0.953
CLR | —– | —– | —– | 0.054 | —– | —– | 0.947 | 0.926 | —– | —– | 0.958
HMGCR Region
Top SNP | —– | 0.001 | 0.018 | 0.048 | 0.049 | 0.019 | 0.953 | 0.776 | 0.099 | 0.019 | 0.924
Pruning | $\rho=0.1$ | 0.001 | 0.018 | 0.048 | 0.049 | 0.019 | 0.953 | 0.776 | 0.099 | 0.019 | 0.924
$\rho=0.3$ | 0.001 | 0.017 | 0.048 | 0.048 | 0.018 | 0.948 | 0.795 | 0.097 | 0.018 | 0.928
$\rho=0.5$ | 0.000 | 0.017 | 0.056 | 0.047 | 0.017 | 0.953 | 0.825 | 0.096 | 0.017 | 0.932
$\rho=0.7$ | 0.000 | 0.016 | 0.052 | 0.046 | 0.016 | 0.938 | 0.820 | 0.093 | 0.017 | 0.917
$\rho=0.9$ | 0.000 | 0.015 | 0.057 | 0.042 | 0.015 | 0.924 | 0.766 | 0.084 | 0.016 | 0.802
PCA | $k=0.99$ | 0.000 | 0.017 | 0.049 | 0.049 | 0.017 | 0.949 | 0.830 | 0.099 | 0.018 | 0.930
$k=0.999$ | 0.000 | 0.016 | 0.052 | 0.048 | 0.017 | 0.953 | 0.826 | 0.098 | 0.017 | 0.934
JAM | $\rho=0.6$ | 0.001 | 0.017 | 0.047 | 0.048 | 0.018 | 0.957 | 0.789 | 0.096 | 0.018 | 0.935
$\rho=0.8$ | 0.000 | 0.017 | 0.039 | 0.047 | 0.018 | 0.961 | 0.777 | 0.096 | 0.018 | 0.942
$\rho=0.9$ | 0.000 | 0.017 | 0.039 | 0.047 | 0.018 | 0.957 | 0.777 | 0.096 | 0.018 | 0.939
$\rho=0.95$ | 0.001 | 0.017 | 0.035 | 0.047 | 0.018 | 0.957 | 0.777 | 0.095 | 0.018 | 0.939
F-LIML | —– | 0.000 | 0.017 | 0.064 | 0.050 | 0.017 | 0.945 | 0.850 | 0.101 | 0.019 | 0.946
CLR | —– | —– | —– | 0.055 | —– | —– | 0.947 | 0.834 | —– | —– | 0.945
Table 1: Simulation results for the “strong instruments” scenario.

Table 1 contains simulation results from the “strong instruments” scenario. In
this scenario, all methods performed quite well in simulations with a null
causal effect. When a positive causal effect was used, PCA, JAM and F-LIML
managed to identify the true value of the causal parameter with decent
accuracy in both regions. LD-pruning did the same in most cases, but the
method’s performance deteriorated for large $\rho$ values, exhibiting bias
towards the null. The bias was more pronounced for $\theta=0.1$ than for
$\theta=0.05$.

The good performance of all methods for $\theta=0$ suggests that any issues in
their performance are due to weak instrument bias, which acts towards the null
and hence would only affect simulations with $\theta\neq 0$. Accordingly, any
biases observed for $\theta\neq 0$ were towards the null.

For LD-pruning, large values of $\rho$ make the inclusion of weak instruments
more likely and the numerical computation and inversion of the correlation
matrix more challenging. For small values of $\rho$ the method is more robust
to the inclusion of weak instruments. On the other hand, when the genetic
region studied contains multiple causal signals, using a small correlation
threshold may discard some of the causal variants from the analysis. In our
simulations, this translated into a fairly small increase in causal standard
errors and a decrease in the method’s power to detect a causal effect. We also
note that in the simulations of Table 1, LD-pruning was implemented after
discarding genetic variants whose p-values did not reach genome-wide
significance. This is a “best-case scenario” for the method in terms of
avoiding weak instrument bias; in practice, pruning is often implemented using
less stringent thresholds (Dudbridge, 2013), and the effects of weak
instrument bias can be more severe.

As expected, an MR analysis using only the variant with the smallest p-value
in the region was unbiased, but had larger standard errors and lower power
than other methods. This was more pronounced for the SHBG region and less so
for the HMGCR region, since genetic correlation were stronger in the HMGCR
region and using only a single variant could partly account for the effects of
other variants through correlation. We also note that in our simulations, all
causal variants had $G-X$ effects in the same direction and of similar
magnitude, and there were no heterogeneous effects towards the outcome. This
is a best-case scenario for single-variant analysis; differences between it
and other methods are likely to be larger in practice.

The performance of the PCA method was similar for $k=99\%$ and $k=99.9\%$ and
was quite accurate. With a causal effect of $\theta=0.1$, the algorithm
exhibited a small reduction in coverage due to weak instrument bias but still
performed better than LD-pruning. In general, the effect of weak instruments
bias on the algorithm is more pronounced for larger values of the tuning
parameter $k$. Here, the standard values of $k=99\%$ and $k=99.9\%$ worked
reasonably well in all simulation scenarios. The power to reject the null
causal hypothesis was greater for $k=99.9\%$, at least for the SHBG region.

The JAM algorithm exhibited similar performance to PCA, with a small
attenuation of causal effect estimates for $\theta=0.1$. JAM requires a
correlation threshold to be specified for the pruning step before running
Bayesian variable selection, but the algorithm’s performance was quite robust
to the value of that tuning parameter, certainly more so than that of LD-
pruning. Its empirical power was slightly higher than PCA for the SHBG region,
but lower for he HMGCR region.

The F-LIML method does not depend on a tuning parameter, as it can
automatically determine the number of latent factors to use. Compared to JAM
and PCA, the algorithm yielded slightly more accurate causal effect estimates
and slightly better calibrated confidence intervals for $\theta=0.1$ but had
slightly inflated Type I error rates for $\theta=0$. The latter issue was
addressed by the conditional likelihood ratio test, at the expense of no
causal effect estimates and a slightly lower power than F-LIML.

|  | $\theta=0$ | $\theta=0.05$ | $\theta=0.1$
---|---|---|---|---
Method | $\hat{\theta}$ | $se(\hat{\theta})$ | Type I | $\hat{\theta}$ | $se(\hat{\theta})$ | Cov | $\hat{\theta}$ | $se(\hat{\theta})$ | Cov | Power
SHBG Region
Top SNP | —– | -0.002 | 0.051 | 0.045 | 0.033 | 0.052 | 0.926 | 0.072 | 0.054 | 0.891 | 0.259
Pruning | $\rho=0.1$ | -0.001 | 0.047 | 0.045 | 0.031 | 0.048 | 0.922 | 0.069 | 0.050 | 0.868 | 0.290
$\rho=0.3$ | -0.002 | 0.040 | 0.045 | 0.032 | 0.041 | 0.921 | 0.066 | 0.042 | 0.843 | 0.359
$\rho=0.5$ | -0.003 | 0.038 | 0.051 | 0.031 | 0.039 | 0.916 | 0.067 | 0.040 | 0.829 | 0.380
$\rho=0.7$ | -0.003 | 0.037 | 0.047 | 0.032 | 0.038 | 0.915 | 0.066 | 0.039 | 0.818 | 0.387
$\rho=0.9$ | -0.003 | 0.035 | 0.049 | 0.031 | 0.036 | 0.907 | 0.062 | 0.037 | 0.773 | 0.359
PCA | $k=0.99$ | 0.000 | 0.040 | 0.050 | 0.035 | 0.041 | 0.924 | 0.068 | 0.042 | 0.876 | 0.392
$k=0.999$ | -0.002 | 0.033 | 0.035 | 0.027 | 0.034 | 0.902 | 0.054 | 0.036 | 0.736 | 0.363
JAM | $\rho=0.6$ | -0.002 | 0.061 | 0.013 | 0.031 | 0.063 | 0.972 | 0.069 | 0.066 | 0.950 | 0.121
$\rho=0.8$ | -0.001 | 0.066 | 0.008 | 0.033 | 0.069 | 0.982 | 0.070 | 0.072 | 0.956 | 0.109
$\rho=0.9$ | -0.002 | 0.070 | 0.008 | 0.034 | 0.071 | 0.983 | 0.071 | 0.074 | 0.959 | 0.101
$\rho=0.95$ | -0.002 | 0.070 | 0.005 | 0.035 | 0.071 | 0.984 | 0.071 | 0.076 | 0.968 | 0.096
F-LIML | —– | -0.003 | 0.037 | 0.183 | 0.051 | 0.039 | 0.782 | 0.101 | 0.041 | 0.802 | 0.667
CLR | —– | —– | —– | 0.041 | —– | —– | 0.939 | —– | —– | 0.948 | 0.404
HMGCR Region
Top SNP | —– | -0.001 | 0.051 | 0.055 | 0.037 | 0.054 | 0.939 | 0.077 | 0.054 | 0.916 | 0.309
Pruning | $\rho=0.1$ | -0.001 | 0.050 | 0.055 | 0.037 | 0.053 | 0.939 | 0.077 | 0.054 | 0.914 | 0.308
$\rho=0.3$ | -0.002 | 0.047 | 0.049 | 0.039 | 0.050 | 0.942 | 0.072 | 0.051 | 0.893 | 0.309
$\rho=0.5$ | -0.001 | 0.045 | 0.054 | 0.038 | 0.047 | 0.941 | 0.072 | 0.048 | 0.884 | 0.326
$\rho=0.7$ | -0.001 | 0.043 | 0.051 | 0.036 | 0.046 | 0.941 | 0.072 | 0.047 | 0.876 | 0.346
$\rho=0.9$ | 0.000 | 0.040 | 0.053 | 0.031 | 0.043 | 0.916 | 0.062 | 0.044 | 0.794 | 0.293
PCA | $k=0.99$ | 0.003 | 0.047 | 0.054 | 0.040 | 0.049 | 0.946 | 0.080 | 0.051 | 0.919 | 0.359
$k=0.999$ | 0.001 | 0.043 | 0.054 | 0.033 | 0.045 | 0.936 | 0.067 | 0.046 | 0.854 | 0.303
JAM | $\rho=0.6$ | -0.002 | 0.058 | 0.022 | 0.037 | 0.064 | 0.975 | 0.073 | 0.064 | 0.958 | 0.166
$\rho=0.8$ | 0.000 | 0.062 | 0.017 | 0.038 | 0.066 | 0.980 | 0.075 | 0.068 | 0.969 | 0.140
$\rho=0.9$ | 0.000 | 0.063 | 0.017 | 0.039 | 0.069 | 0.984 | 0.079 | 0.069 | 0.974 | 0.135
$\rho=0.95$ | 0.000 | 0.064 | 0.011 | 0.040 | 0.070 | 0.985 | 0.079 | 0.071 | 0.975 | 0.131
F-LIML | —– | 0.002 | 0.047 | 0.165 | 0.051 | 0.050 | 0.859 | 0.105 | 0.054 | 0.862 | 0.524
CLR | —– | —– | —– | 0.049 | —– | —– | 0.955 | —– | —– | 0.952 | 0.331
Table 2: Simulation results for the “weak instruments” scenario.

Table 2 reports simulation results from the “weak instruments” scenario. Weak
instruments bias had a much higher impact in these simulations, with several
methods facing attenuation of their causal effect estimates. As expected, any
bias occurred only for $\theta\neq 0$ and was towards the null, while all
methods also had increased standard errors and low power.

Top-SNP analysis, LD-pruning, PCA and JAM all suffered from weak instrument
bias in this scenario. The magnitude of bias was similar for these methods.
For pruning and PCA, the bias caused poor coverage properties for confidence
intervals and Type I error rates above nominal levels. JAM only selected a
small number of genetic variants (it selected an average of $1.6$ variants per
run) and attempted to adjust for the presence of weak instruments by producing
wider confidence intervals; in fact it was rather conservative in our
simulations, with coverage rates above $95\%$, and consequently its power was
quite low. In terms of their dependence on tuning parameters, the methods
behaved similar to the “strong instruments” scenario: LD-pruning was more
susceptible to bias for large correlation thresholds, while having smaller
standard errors. PCA and JAM remained robust to the specification of their
tuning parameters.

The F-LIML method had the opposite performance compared to JAM: the algorithm
provided quite accurate causal effect estimates, but underestimated standard
errors and resulted in confidence intervals with inflated Type I error rates
and below nominal coverage. The inefficiency of the algorithm in weak
instrument scenarios has been noted in Patel et al. (2020). On the other hand,
the conditional likelihood ratio test was the method least affected by weak
instruments bias and offered a reliable way of assessing the existence of a
causal effect. Its power was worse than the (overprecise) F-LIML method but
better than that of the JAM algorithm.

In summary, our simulations suggest that simple methods such as LD-pruning and
single-variant MR analysis can be reliable in simulations where weak
instrument bias is of lesser concern, but should be used with caution in cis-
MR analyses where weak instrument bias is suspected. Even when LD-pruning is
used, MR causal effect estimates should be computed and reported for a range
of different correlation thresholds as a form of sensitivity analysis. Between
JAM, PCA and F-LIML, the differences were small in the “strong instruments”
simulation. With weak instruments, F-LIML provided accurate causal effect
estimates but poor uncertainty quantification, JAM yielded biased causal
effect estimates, reasonable confidence intervals but low power, the CLR test
offered nominal confidence intervals and decent power at the expense of no
point estimates and PCA was the method most affected by weak instruments bias.
A theoretical advantage of JAM compared to the other approaches is that JAM’s
variable selection gives an indication of which genetic variants are more
likely to be causally associated with the risk factor, although it is not
clear whether this additional information would be useful in an MR study where
the objective is to perform inference about the $X-Y$ causal effect. Finally,
it may be possible to improve the performance of some of the methods by
removing very weak variants (e.g. with $p>0.05$) from the analysis prior to
implementing the methods, but we have not considered that here.

### Additional Simulations

In supplementary material, we report results from a range of additional
simulation scenarios. Specifically, we conducted three additional simulations,
on which we comment here briefly.

First, we used a larger sample size of $N_{1}=100000$ from which to obtain
risk factor-outcome summary statistics. This could be reminiscent of a cis-MR
analysis using a downstream biomarker as a proxy for protein expression, and
obtaining $G-X$ summary data from a large-scale GWAS. Our results
(Supplementary Table 1) suggested that the various methods perform similarly
in this case and issues such as weak instrument bias and robustness to the
choice of the pruning threshold are less concerning with large sample sizes.

In the second simulation, we experimented with different numbers of causal
variants. We used the SHBG region, simulated under the “strong instruments”
scenario, and assumed that the region contains either one or three genetic
variants causally associated with the risk factor. The performance of the
various methods (Supplementary Table 2) was similar to that observed in our
baseline simulations with six causal variants, suggesting that the number of
causal signals in the region has little impact on their performance relative
to each other, except for the single-variant approach which did expectedly
better in simulations with only one causal variant.

In our third supplementary simulation, we assessed the impact of only having
access to a small reference dataset on the various methods. In each iteration,
we bootstrapped the UK Biobank data in order to obtain a small reference
sample of size $N_{ref}=1000$ and implemented cis-MR using that reference
dataset (Supplementary Table 3). The JAM algorithm was the most sensitive
method to this change. The algorithm selected more variants on average,
exhibited slight differences in its performance for different values of the
correlation threshold and performed poorly for $\rho\geq 0.9$. In such
datasets, we recommend using a more stringent correlation threshold for the
algorithm. The other methods were less affected in comparison.

We refer to the supplement for more details on these simulations.

## Applications

### Causal Effect of LDL-cholesterol on CHD Risk Using Variants in the HMGCR
Region

#### Introduction

We now compare the various cis-MR methods in two real-data applications. In
the first application, we use Mendelian randomization to investigate the
relationship between low density lipoprotein (LDL) cholesterol concentration
and the risk of coronary heart disease (CHD) using genetic variants from the
HMGCR region. The HMGCR region is located in chromosome $5$ and encodes the
protein HMG coenzyme-A reductase. The protein plays an important role in the
metabolic pathway that produces cholesterol, and its inhibition by statins is
a common treatment to reduce LDL-cholesterol levels.

This analysis is performed mainly for illustrative purposes: the causal
connection between LDL-cholesterol and CHD risk has already been explored in
several papers in the literature. Genome-wide MR studies of the effect of LDL-
cholesterol on CHD risk have been performed by Linsel-Nitschke et al. (2008);
Waterworth et al. (2010); Burgess et al. (2014); Holmes et al. (2014); Allara
et al. (2019) among others. In the field of cis-MR, Swerdlow et al. (2015)
considered associations between HMGCR variants and a range of biomarkers in
order to investigate whether statin treatment increases the risk of type 2
diabetes. Ference et al. (2016) compared the effects on type 2 diabetes for
variants in the HMGCR and PCSK9 regions, while Ference et al. (2019) compared
the LDL-CHD causal effects obtained using variants in the HMGCR and ACLY
regions in order to explore the suitability of ACLY as a drug target. The
HMGCR region has also been used as an applied example in methodological papers
(Schmidt et al., 2019).

#### Datasets and Methods

Since a clear link between HMGCR protein expression and lowering LDL-
cholesterol has been established, we used LDL-cholesterol as a risk factor
instead of protein expression. We estimated genetic associations between
variants in the HMGCR region and LDL-cholesterol using data from the UK
Biobank. Genetic associations were computed based on a sample of
$N_{1}=349795$ unrelated individuals of European ancestry. Associations with
coronary heart disease risk were obtained from the CARDIoGRAMplusC4D
consortium (CARDIoGRAMplusC4D Consortium, 2015), based on a sample of
$N_{2}=184305$ individuals. Finally, genetic correlations were computed using
individual-level data for $367643$ unrelated Europeans from the UK Biobank
(the sample size for genetic associations with LDL-cholesterol was slightly
smaller than the reference sample due to missing values).

We used a wider region than in our simulations and extracted information on
variants within 500Kb pairs from the HMGCR gene. In total, $2915$ variants
were present in both the UK Biobank and the CARDIoGRAMplusC4D dataset. Of
those $2915$ variants, $20$ were discarded because they either had missing
associations with LDL-cholesterol in the UK Biobank dataset or concerned
multiple alleles on the same locus that were not detected in both datasets. We
did not discard variants with low effect allele frequencies. Our analysis was
therefore based on $P=2895$ genetic variants. In Figure 2 we visualize the
genetic correlations in the HMGCR region, and present a Manhattan plot of
genetic association with LDL-cholesterol levels.

Figure 2: Left: genetic correlations in the HMGCR region. Right: Manhattan
plot of p-values for $G-X$ associations.

We then conducted cis-MR using the minimum p-value variant, LD-pruning,
principal components analysis, the JAM algorithm and the factor-based methods.
The various methods were implemented using a range of different parameter
specifications, similar to those used for our simulation study. For stepwise
pruning, we used the values $\rho=0.1,0.3,0.5,0.7,0.9$ for the correlation
threshold and a GWAS significance threshold of $\tau=5\times 10^{-8}$. Of the
$2895$ variants in the region, $1424$ had p-values below the GWAS threshold.
For the PCA method, we used principal components that explained $k=99\%$ or
$k=99.9\%$ of variation in the genetic data. For JAM, we implemented a pre-
pruning step with a correlation threshold of $\rho=0.6,0.8,0.9,0.95$ and no
significance threshold. For F-LIML, we allowed the algorithm to determine the
number of latent factors to use.

#### Results

Method | $\hat{\theta}$ | $se(\hat{\theta})$ |  $95\%$ CI
---|---|---|---
Top-SNP | —– | 0.555 | 0.167 | (0.228 , 0.882)
Pruning | $\rho=0.1$ | 0.389 | 0.136 | (0.122 , 0.656)
$\rho=0.3$ | 0.358 | 0.128 | (0.108 , 0.608)
$\rho=0.5$ | 0.305 | 0.120 | (0.069 , 0.541)
$\rho=0.7$ | 0.275 | 0.107 | (0.066 , 0.485)
$\rho=0.9$ | 0.378 | 0.069 | (0.242 , 0.513)
PCA | $k=0.99$ | -0.030 | 0.026 | (-0.081 , 0.020)
$k=0.999$ | 0.000 | 0.022 | (-0.043 , 0.042)
JAM | $\rho=0.6$ | 0.401 | 0.168 | (0.072 , 0.731)
$\rho=0.8$ | 0.355 | 0.128 | (0.104 , 0.606)
$\rho=0.9$ | 0.358 | 0.126 | (0.111 , 0.606)
$\rho=0.95$ | 0.358 | 0.117 | (0.129 , 0.586)
F-LIML | —– | 0.501 | 0.165 | (0.178 , 0.824)
CLR | —– | —– | —– | (0.169 , 0.839)
Table 3: Cis-MR analysis of the effect of LDL-cholesterol on CHD risk, using
genetic variants in the HMGCR region.

Results of our analysis are presented in Table 3. We report causal effect
estimates, standard errors and $95\%$ confidence intervals obtained by each
method. The reported estimates represent log-odds ratios of increase in CHD
risk per standard deviation increase in LDL-cholesterol levels.

Genetically elevated LDL-cholesterol concentration based on variants in the
HMGCR region is known to be associated with increased risk of coronary heart
disease Cholesterol Treatment Trialists (CTT) Collaboration (2010); Ference et
al. (2019). This was confirmed by the pruning, JAM, F-LIML and CLR methods.
For example, the JAM algorithm using a pre-pruning threshold of $\rho=0.8$
suggested a log-odds ratio of $0.355$ (odds ratio $1.426$, $95\%$ confidence
interval $(1.110,1.833)$). Interestingly, the principal components method
suggested a null effect. The main difference in the implementation of the
method compared to our simulation study was the larger number of genetic
variants used here. Since the presence of many non-causal variants can
exacerbate bias due to weak instruments, we conjecture that this was the issue
behind the method’s poor performance. To confirm our suspicions, we
implemented the PCA method using only the $1424$ GWAS-significant variants in
the dataset. The method produced much more reasonable results, in line with
other methods: the log-odds ratio estimates were $0.423$ for $k=99\%$ and
$0.448$ for $k=99.9\%$ and the null causal hypothesis was rejected on both
occasions.

The remaining methods were more consistent in their results. LD-pruning
exhibited the usual attenuation towards the null for $\rho=0.7$ but yielded a
larger causal effect estimate for $\rho=0.9$, as a result of numerical errors.
Bias due to numerical errors is more likely to be a concern here than in our
simulation study, due to the larger number of variants used. For the JAM
algorithm, there were small differences between implementations for $\rho=0.6$
and larger values of the correlation threshold, but overall the algorithm was
more robust to the choice of that threshold than LD-pruning. The algorithm
suggested the existence of $3-5$ genetic variants independently associated
with LDL-cholesterol in the region. The F-LIML method suggested a slightly
higher causal effect than the other methods, although differences were within
the margin of statistical error. As an additional form of sensitivity
analysis, we implemented both JAM and F-LIML using only GWAS-significant
variants; results were similar to those of Table 3.

### Causal Effect of Testosterone on CHD Risk Using Variants in the SHBG
Region

#### Introduction

In our second application, we apply Mendelian randomization using variants in
the SHBG gene region in order to assess the causal relationship between
testosterone levels and coronary heart disease risk. The SHBG region is
located in chromosome $17$ and encodes sex-hormone binding globulin, a protein
that inhibits the function of sex hormones such as testosterone and estradiol.
The region has been shown to contain strong genetic associations with
testosterone levels, as well as a number of sex hormone traits (Jin et al.,
2012; Coviello et al., 2012; Schooling et al., 2018; Ruth et al., 2020). In
addition, previous research has suggested that the region is likely to contain
several genetic variants independently associated with testosterone; Coviello
et al. (2012) claimed that as many as nine independent signals may be present.
This implies that a naive approach using only the variant with the smallest
p-value would underestimate the genetic contributions to testosterone levels.
The causal relationship between sex hormone traits and CHD risk using variants
in the SHBG region has been studied by Burgess et al. (2017) and Schooling et
al. (2018) with evidence mostly suggesting no causal relationship. Here, we
aimed to replicate the analysis of Burgess et al. (2017) using larger
datasets.

#### Datasets and Methods

We used serum testosterone levels as the exposure for our MR analysis. We
obtained genetic associations with testosterone using summary-level data from
the Neale Lab
website222https://docs.google.com/spreadsheets/d/1kvPoupSzsSFBNSztMzl04xMoSC3Kcx3CrjVf4yBmESU/edit?ts=5b5f17db#gid=227859291.
These summary data were derived from the UK Biobank, using a sample of
$N_{1}=312102$ unrelated individuals of European descent. We defined the
genetic region to include genetic variants within 500Kb pairs on either side
of the SHBG gene.

As in the previous application, we obtained genetic associations with CHD risk
from the CARDIoGRAMplusC4D dataset (CARDIoGRAMplusC4D Consortium, 2015), using
a sample of $N_{2}=184305$ individuals. A total of $3053$ genetic variants
were present in both datasets and were therefore included in our analysis. We
did not conduct separate analyses on males and females, since we did not have
access to sex-specific associations with CHD risk. Finally, we used the UK
Biobank as a reference dataset from which to extract LD correlations between
genetic variants. Figure 3 presents a Manhattan plot of associations with
testosterone, as well as the genetic correlations in the region.

Figure 3: Left: genetic correlations in the SHBG region. Right: Manhattan plot
of p-values for $G-X$ associations.

We then implemented the various cis-MR methods in order to select genetic
instruments and assess whether SHBG variants are causally associated with CHD
risk. We used the same setting as in the HMGCR application for the tuning
parameters of each method. The stepwise pruning method was implemented using
only variants with GWAS-significant associations with serum testosterone
levels. A total of $1156$ variants had p-values lower than the $5\times
10^{-8}$ threshold. The remaining methods were implemented using all $3053$
variants in the region.

#### Results

Method | $\hat{\theta}$ | $se(\hat{\theta})$ |  $95\%$ CI
---|---|---|---
Top-SNP | —– | -0.071 | 0.036 | (-0.141 , -0.001)
Pruning | $\rho=0.1$ | -0.040 | 0.025 | (-0.088 , 0.008)
$\rho=0.3$ | -0.031 | 0.023 | (-0.076 , 0.014)
$\rho=0.5$ | -0.040 | 0.021 | (-0.081 , 0.001)
$\rho=0.7$ | -0.033 | 0.015 | (-0.062 , -0.003)
$\rho=0.9$ | 0.032 | 0.009 | (0.014 , 0.051)
PCA | $k=0.99$ | -0.052 | 0.029 | (-0.109 , 0.004)
$k=0.999$ | -0.018 | 0.024 | (-0.065 , 0.028)
JAM | $\rho=0.6$ | -0.042 | 0.024 | (-0.089 , 0.005)
$\rho=0.8$ | -0.037 | 0.024 | (-0.084 , 0.010)
$\rho=0.9$ | -0.040 | 0.024 | (-0.087 , 0.006)
$\rho=0.95$ | -0.035 | 0.024 | (-0.081 , 0.012)
F-LIML | —– | -0.037 | 0.025 | (-0.087 , 0.013)
CLR | —– | —– | —– | (-0.087 , 0.014)
Table 4: Mendelian randomization analysis of the effect of serum testosterone
levels on CHD risk, using genetic variants in the SHBG region.

Table 4 reports causal effect estimates, standard errors and $95\%$ confidence
intervals for each method. Once again, estimates are reported in the log-odds
ratio scale and represent increases in CHD risk per standard deviation
increase in serum testosterone levels. All methods with the exception of top-
SNP and LD-pruning at $0.7$ or $0.9$ suggested no causal relationship between
testosterone and CHD risk based on the SHBG region, although point estimates
were consistently in the risk-decreasing direction.

In this example, the performance of the various methods was less seriously
affected by weak instruments bias because there was apparently no causal
relationship between serum testosterone levels and CHD risk. This is
consistent with our simulation design, where no bias was observed under the
“null causal effect” scenario.

A notable inconsistency was that an MR analysis using only the genetic variant
with the smallest p-value in the region suggested a risk-decreasing,
statistically significant causal effect. This may be suggestive of
heterogeneity in evidence from different variants within the SHBG region.
Therefore, this example empirically demonstrates the pitfalls of using
simplistic single variant analyses when in fact multiple signals exist within
a region. The top variant in our analysis was rs1799941, which is known to be
associated with testosterone levels (Ruth et al., 2020).

The LD-pruning method was rather inconclusive in this application.
Implementations with a low correlation threshold suggested no causal effect.
However, the method suggested an effect in the risk-decreasing direction for
$\rho=0.7$ and in the risk-increasing direction for $\rho=0.9$. This was
combined with a rather sharp increase in precision around the causal
estimates. As in our previous application, this is indicative of numerical
issues in computing the correlated-instruments IVW estimate and its standard
error.

The principal components approach suggested no causal association for both
values of its tuning parameter. JAM did the same, and was once again
consistent with respect to the value of its tuning parameter. The algorithm
suggested the existence of about $8-9$ independent signals in the region,
confirming claims by Coviello et al. (2012). The F-LIML estimate and the CLR
confidence interval were in line with the results obtained by JAM.

Our results obtained in this application are similar to those reported in
Burgess et al. (2017). Similar to the results reported in Table 4, pruning
estimates in Burgess et al. (2017) suggested no causal effect for small
correlation thresholds bt were unstable for large $\rho$ (in fact, they were
more unstable than in our results, possibly due to the smaller sample sizes
used in that paper). Estimates from the PCA method were in the risk-decreasing
direction but did not achieve statistical significance. Overall, both Burgess
et al. (2017) and our analysis suggested no causal association between
testosterone and CHD risk based on variants in the SHBG region; further
evidence of no causality was provided by Schooling et al. (2018), which also
used data from the UK Biobank.

## Discussion

Cis-Mendelian randomization is emerging as a widely applicable approach to
support drug development and inform clinical trials. In this manuscript we
have outlined and compared a range of methods from the cis-MR literature to
select variants from a genetic region of interest and estimate the MR causal
effect. In particular, we have considered MR using only the minimum p-value
variant, and also LD-pruning, which is one of the most commonly used
approaches in practice, principal components analysis, factor-based methods
and stochastic-search variable selection via the JAM algorithm. We compared
the various methods across a range of simulation scenarios, which were based
on real genetic correlation structures of two different gene regions drawn
from the UK Biobank. We also assessed the performance of the different methods
in two cis-MR case studies, investigating the effect of LDL-cholesterol on CHD
risk using variants in the HMGCR region and the effect of testosterone on CHD
risk using variants in the SHBG region.

We have found that LD-pruning can be reliable when using large sample sizes
and strong genetic instruments, but can be susceptible to weak instrument
bias. Moreover, results obtained using LD-pruning can be sensitive to the
correlation threshold used, and using a high correlation threshold can result
in numerical instabilities. Based on these results, we would therefore
recommend that when the pruning method is used, it is implemented using a
range of correlation thresholds as a form of sensitivity analysis.

Although not completely free from weak instrument bias, methods such as PCA,
F-LIML, the JAM algorithm and the conditional likelihood ratio test were less
affected by such bias than LD-pruning. JAM was also more robust to the choice
of its tuning parameter than pruning, while PCA and F-LIML achieved good
performance in a range of different scenarios with the default values for
their tuning parameters. Therefore, we recommend that some of these methods
are implemented as sensitivity tools in applied cis-MR analyses.

Compared to each other, JAM, PCA, F-LIML and CLR exhibited similar performance
in simulations with strong genetic instruments. With weak instruments, the CLR
test was least affected by weak instruments bias and provided a valid test for
the causal null hypothesis, albeit no point estimates. F-LIML obtained
accurate causal effect estimates but yielded confidence intervals with poor
coverage, while JAM had biased causal effect estimates but better uncertainty
quantification. The principal components IVW approach was affected by weak
instruments bias to a greater extent than the other methods, although our
real-data applications suggest that PCA can still perform well if very weak
variants are removed from the analysis prior to its implementation.

The use of multiple methods as a form of sensitivity analysis can increase the
reliability of results of a cis-MR study. This can be further reinforced by
triangulating evidence from different study designs (Gill et al., 2020). This
includes colocalization, an approach which has been shown to yield benefits
when used in conjunction with MR to prioritize proteins for drug development
(Zheng et al., 2019). Here, we have only focused on Mendelian randomization
and not considered alternative approaches.

Our analysis has a number of limitations. Our simulations were by no means
exhaustive; for example, they were based on only two genetic regions. In
addition, our simulations have not considered other forms of bias that may
arise in cis-MR applications, such as selection bias, population
stratification or pleiotropy. Bias due to pleiotropy is worth discussing in
more detail, as it is one of the main concerns in many MR studies. Pleiotropy
is less of a concern in cis-MR compared to traditional MR studies (Schmidt et
al., 2019), but would arise if different variants from the same region were
involved in different biological pathways relevant to the outcome, resulting
in heterogeneous causal effect estimates. When pleiotropy does occur, standard
algorithms for pleiotropy-robust MR (including the JAM algorithm) can be used
to address it. However, care must be taken when extrapolating the assumptions
made by such algorithms to the framework of many correlated genetic variants.
For example, suppose that a region contains a single pleiotropic variant, but
that variant is strongly correlated with half of the other variants in the
region. This will induce bias in the marginal causal effect estimates of the
correlated variants, and hence pleiotropy-robust methods requiring more than
$50\%$ of the genetic variants in the region to be valid instruments will be
biased.

Cis-Mendelian randomization is an active area of research, both methodological
and applied. As the field develops further, it is likely that new statistical
approaches will be created to aid applied investigators in their work. By
providing an overview and evaluation of existing methods, we hope that our
current paper will contribute to that direction.

## Acknowledgements

We would like to thank Ashish Patel for many helpful discussions and comments.
This work was supported by the UK Medical Research Council (Core Medical
Research Council Biostatistics Unit Funding Code: MC UU 00002/7). Apostolos
Gkatzionis and Paul Newcombe were supported by a Medical Research Council
Methodology Research Panel grant (Grant Number RG88311). Stephen Burgess was
supported by a Sir Henry Dale Fellowship jointly funded by the Wellcome Trust
and the Royal Society (Grant Number 204623/Z/16/Z). This research has been
conducted using the UK Biobank Resource under Application Number 30931.

## References

  * Allara et al. (2019) Allara, E., G. Morani, P. Carter, A. Gkatzionis, V. Zuber, C. N. Foley, J. M. Rees, A. M. Mason, S. Bell, D. Gill, S. Lindstrom, A. S. Butterworth, E. Di Angelantonio, J. Peters, S. Burgess, and INVENT consortium (2019).  Genetic determinants of lipids and cardiovascular disease outcomes: A wide-angled Mendelian randomization investigation.  Circulation. Genomic and precision medicine 12(12), e002711.
  * Anderson and Rubin (1949) Anderson, T. W. and H. Rubin (”1949, 3).  Estimation of the parameters of a single equation in a complete system of stochastic equations.  Annals of Mathematical Statistics 20(1), 46–63.
  * Asimit et al. (2019) Asimit, J. L., D. B. Rainbow, M. D. Fortune, N. F. Grinberg, L. S. Wicker, and C. Wallace (2019).  Stochastic search and joint fine-mapping increases accuracy and identifies previously unreported associations in immune-mediated diseases.  Nature Communications 10.
  * Bai and Ng (2002) Bai, J. and S. Ng (2002).  Determining the number of factors in approximate factor models.  Econometrica 70(1), 191–221.
  * Benner et al. (2016) Benner, C., C. C. A. Spencer, A. S. Havulinna, V. Salomaa, S. Ripatti, and M. Pirinen (2016).  FINEMAP: efficient variable selection using summary data from genome-wide association studies.  Bioinformatics 32(10), 1493–1501.
  * Boef et al. (2015) Boef, A. G. C., O. M. Dekkers, and S. le Cessie (2015).  Mendelian randomization studies: a review of the approaches used and the quality of reporting.  International Journal of Epidemiology 44(2), 496–511.
  * Burgess et al. (2016) Burgess, S., N. M. Davies, and S. G. Thompson (2016).  Bias due to participant overlap in two-sample Mendelian randomization.  Genetic Epidemiology 40(7), 597–608.
  * Burgess et al. (2016) Burgess, S., F. Dudbridge, and S. G. Thompson (2016).  Combining information on multiple instrumental variables in Mendelian randomization: comparison of allele score and summarized data methods.  Statistics in Medicine 35(11), 1880–1906.
  * Burgess et al. (2014) Burgess, S., D. F. Freitag, H. Khan, D. N. Gorman, and S. G. Thompson (2014, 10).  Using multivariable Mendelian randomization to disentangle the causal effects of lipid fractions.  PLOS ONE 9(10), 1–8.
  * Burgess and Thompson (2015) Burgess, S. and S. G. Thompson (2015).  Mendelian randomization: methods for using genetic variants in causal estimation.  Chapman & Hall, Boca Raton, FL.
  * Burgess et al. (2011) Burgess, S., S. G. Thompson, and CRP CHD Genetics Collaboration (2011, 03).  Avoiding bias from weak instruments in Mendelian randomization studies.  International Journal of Epidemiology 40(3), 755–764.
  * Burgess et al. (2017) Burgess, S., V. Zuber, E. Valdes-Marquez, B. B. Sun, and J. C. Hopewell (2017).  Mendelian randomization with fine‐mapped genetic data: Choosing from large numbers of correlated instrumental variables.  Genetic Epidemiology 41(4), 714–725.
  * CARDIoGRAMplusC4D Consortium (2015) CARDIoGRAMplusC4D Consortium (2015).  A comprehensive 1000 Genomes-based genome-wide association meta-analysis of coronary artery disease.  Nature Genetics 47(10), 1121–1130.
  * CCGC (CRP CHD Genetics Collaboration) (2011) CCGC (CRP CHD Genetics Collaboration) (2011).  Association between C reactive protein and coronary heart disease: Mendelian randomisation analysis based on individual participant data.  British Medical Journal 342.
  * Cholesterol Treatment Trialists (CTT) Collaboration (2010) Cholesterol Treatment Trialists (CTT) Collaboration (2010).  Efficacy and safety of more intensive lowering of LDL cholesterol: a meta-analysis of data from 170 000 participants in 26 randomised trials.  The Lancet 376(9753), 1670–1681.
  * Coviello et al. (2012) Coviello, A. D., R. Haring, M. Wellons, D. Vaidya, T. Lehtimäki, S. Keildson, and Other Authors (2012).  A genome-wide association meta-analysis of circulating sex hormone–binding globulin reveals multiple loci implicated in sex steroid hormone regulation.  PLOS Genetics 8(7), 1–12.
  * Davey Smith and Ebrahim (2003) Davey Smith, G. and S. Ebrahim (2003).  ‘Mendelian randomization’: can genetic epidemiology contribute to understanding environmental determinants of disease?  International Journal of Epidemiology 32(1), 1 – 22.
  * Dudbridge (2013) Dudbridge, F. (2013, 03).  Power and predictive accuracy of polygenic risk scores.  PLOS Genetics 9(3), 1–17.
  * Dudbridge and Newcombe (2015) Dudbridge, F. and P. J. Newcombe (2015).  Accuracy of gene scores when pruning markers by linkage disequilibrium.  Human Heredity 80(4), 178 – 186.
  * Ference et al. (2019) Ference, B. A., K. K. Ray, A. L. Catapano, T. B. Ference, S. Burgess, D. R. Neff, C. Oliver-Williams, A. M. Wood, A. S. Butterworth, E. Di Angelantonio, J. Danesh, J. J. P. Kastelein, and S. J. Nicholls (2019).  Mendelian Randomization Study of ACLY and Cardiovascular Disease.  New England Journal of Medicine 380(11), 1033–1042.
  * Ference et al. (2016) Ference, B. A., J. G. Robinson, R. D. Brook, A. L. Catapano, M. J. Chapman, D. R. Neff, S. Voros, R. P. Giugliano, G. Davey Smith, S. Fazio, and M. S. Sabatine (2016).  Variation in PCSK9 and HMGCR and risk of cardiovascular disease and diabetes.  New England Journal of Medicine 375(22), 2144–2153.
  * Gill et al. (2020) Gill, D., M. K. Georgakis, V. M. Walker, A. F. Schmidt, A. Gkatzionis, D. F. Freitag, C. Finan, A. D. Hingorani, J. C. M. Howson, S. Burgess, D. I. Swerdlos, M. Dichgans, B. M. Psaty, G. Davey Smith, M. V. Holmes, R. A. Scott, J. Zheng, and N. M. Davies (2020).  Mendelian randomization for studying the effects of perturbing drug targets: a practical guide.  In preparation.
  * Gkatzionis et al. (2020) Gkatzionis, A., S. Burgess, D. V. Conti, and P. J. Newcombe (2020).  Bayesian variable selection with a pleiotropic loss function in Mendelian randomization.  bioRxiv.
  * Hay et al. (2014) Hay, M., D. W. Thomas, J. L. Craighead, C. Economides, and J. Rosenthal (2014).  Clinical development success rates for investigational drugs.  Nature Biotechnology 32, 40–51.
  * Hocking (1976) Hocking, R. R. (1976).  A biometrics invited paper. the analysis and selection of variables in linear regression.  Biometrics 32(1), 1–49.
  * Holmes et al. (2014) Holmes, M. V., F. W. Asselbergs, T. M. Palmer, F. Drenos, M. B. Lanktree, C. P. Nelson, and Other Authors (2014, 01).  Mendelian randomization of blood lipids for coronary heart disease.  European Heart Journal 36(9), 539–550.
  * Jiang et al. (2020) Jiang, L., S. Xu, N. Mancuso, P. J. Newcombe, and D. V. Conti (2020).  A hierarchical approach using marginal summary statistics for multiple intermediates in a Mendelian randomization or transcriptome analysis.  bioRxiv.
  * Jin et al. (2012) Jin, G., J. Sun, S.-T. Kim, J. Feng, Z. Wang, S. Tao, Z. Chen, L. Purcell, S. Smith, W. B. Isaacs, R. S. Rittmaster, S. L. Zheng, L. D. Condreay, and J. Xu (2012).  Genome-wide association study identifies a new locus JMJD1C at 10q21 that may influence serum androgen levels in men.  Human Molecular Genetics 21(23), 5222–5228.
  * Kamstrup et al. (2009) Kamstrup, P. R., A. Tybjærg-Hansen, R. Steffensen, and B. G. Nordestgaard (2009).  Genetically elevated lipoprotein(a) and increased risk of myocardial infarction.  JAMA 301(22), 2331–2339.
  * Kleibergen (2002) Kleibergen, F. (2002).  Pivotal statistics for testing structural parameters in instrumental variables regression.  Econometrica 70(5), 1781–1803.
  * Krauss et al. (2008) Krauss, R. M., L. M. Mangravite, J. D. Smith, M. W. Medina, D. Wang, X. Guo, M. J. Rieder, J. A. Simon, S. B. Hulley, D. Waters, M. Saad, P. T. Williams, K. D. Taylor, H. Yang, D. A. Nickerson, and J. I. Rotter (2008).  Variation in the 3-hydroxyl-3-methylglutaryl coenzyme A reductase gene is associated with racial differences in low-density lipoprotein cholesterol response to simvastatin treatment.  Circulation 117(12), 1537–1544.
  * Linsel-Nitschke et al. (2008) Linsel-Nitschke, P., A. Götz, J. Erdmann, I. Braenne, P. Braund, C. Hengstenberg, and Other Authors (2008, 08).  Lifelong reduction of LDL-cholesterol related to a common variant in the ldl-receptor gene decreases the risk of coronary artery disease—a Mendelian randomisation study.  PLOS ONE 3(8), 1–9.
  * Miller (2002) Miller, A. J. (2002).  Subset Selection in Regression.  Monographs on Statistics and Applied Probability. London: Chapman and Hall.
  * Moreira (2003) Moreira, M. J. (2003).  A conditional likelihood ratio test for structural models.  Econometrica 71(4), 1027–1048.
  * Nelson et al. (2015) Nelson, M. R., H. Tipney, J. L. Painter, J. Shen, P. Nicoletti, Y. Shen, A. Floratos, P. C. Sham, M. J. Li, J. Wang, L. R. Cardon, J. C. Whittaker, and P. Sanseau (2015).  The support of human genetic evidence for approved drug indications.  Nature Genetics 47(8), 856–860.
  * Newcombe et al. (2016) Newcombe, P. J., D. V. Conti, and S. Richardson (2016).  JAM: A scalable Bayesian framework for joint analysis of marginal SNP effects.  Genetic Epidemiology 40(3), 188 – 201.
  * Patel et al. (2020) Patel, A., S. Burgess, D. Gill, and P. J. Newcombe (2020).  Inference with many correlated weak instruments and summary statistics.  arXiv:2005.01765.
  * Ray et al. (2019) Ray, K. K., H. E. Bays, A. L. Catapano, N. D. Lalwani, L. T. Bloedon, L. R. Sterling, P. L. Robinson, and C. M. Ballantyne (2019).  Safety and efficacy of bempedoic acid to reduce LDL cholesterol.  New England Journal of Medicine 380(11), 1022–1032.
  * Ruth et al. (2020) Ruth, K. S., F. R. Day, D. J. T. Jessica Tyrrel and, A. R. Wood, A. Mahajan, R. N. Beaumont, L. Wittemans, S. Martin, A. S. Busch, A. M. Erzurumluoglu, B. Hollis, T. A. O’Mara, The Endometrial Cancer Association Consortium, M. I. McCarthy, C. Langenberg, D. F. Easton, N. J. Wareham, S. Burgess, A. Murray, K. K. Ong, T. M. Frayling, and J. R. B. Perry (2020).  Using human genetics to understand the disease impacts of testosterone in men and women.  Nature Medicine 26(2), 252–258.
  * Schmidt et al. (2019) Schmidt, A. F., C. Finan, M. Gordillo-Marañón, F. W. Asselbergs, D. F. Freitag, R. S. Patel, B. Tyl, S. Chopade, R. Faraway, M. Zwierzyna, and A. D. Hingorani (2019).  Genetic drug target validation using Mendelian randomization.  bioRxiv.
  * Schmidt et al. (2017) Schmidt, A. F., D. I. Swerdlow, M. V. Holmes, R. S. Patel, Z. Fairhurst-Hunter, D. M. Lyall, and Other Authors (2017).  PCSK9 genetic variants and risk of type 2 diabetes: a Mendelian randomisation study.  The Lancet Diabetes and Endocrinology 5(2), 97–105.
  * Schooling et al. (2018) Schooling, C. M., S. Luo, S. L. A. Yeung, D. J. Thompson, S. Karthikeyan, T. R. Bolton, A. M. Mason, E. Ingelsson, and S. Burgess (2018).  Genetic predictors of testosterone and their associations with cardiovascular disease and risk factors: A Mendelian randomization investigation.  International Journal of Cardiology 267, 171 – 176.
  * Sofat et al. (2010) Sofat, R., A. D. Hingorani, L. Smeeth, S. E. Humphries, P. J. Talmud, J. Cooper, and Other Authors (2010).  Separating the mechanism-based and off-target actions of cholesteryl ester transfer protein inhibitors with cetp gene polymorphisms.  Circulation 121(1), 52–62.
  * Swerdlow et al. (2016) Swerdlow, D. I., K. B. Kuchenbaecker, S. Shah, R. Sofat, M. V. Holmes, J. White, J. S. Mindell, M. Kivimaki, E. J. Brunner, J. C. Whittaker, J. P. Casas, and A. D. Hingorani (2016, 06).  Selecting instruments for Mendelian randomization in the wake of genome-wide association studies.  International Journal of Epidemiology 45(5), 1600–1616.
  * Swerdlow et al. (2015) Swerdlow, D. I., D. Preiss, K. B. Kuchenbaecker, M. V. Holmes, J. E. L. Engmann, and Other Authors (2015).  HMG-coenzyme A reductase inhibition, type 2 diabetes, and bodyweight: evidence from genetic analysis and randomised trials.  The Lancet 385(9965), 351–361.
  * The Interleukin-6 Receptor Mendelian Randomisation Analysis (IL6R MR) Consortium (2012) The Interleukin-6 Receptor Mendelian Randomisation Analysis (IL6R MR) Consortium (2012).  The interleukin-6 receptor as a target for prevention of coronary heart disease: a Mendelian randomisation analysis.  The Lancet 379(9822), 1214–1224.
  * Tibshirani (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society. Series B (Methodological) 58(1), 267–288.
  * Vilhjalmsson et al. (2015) Vilhjalmsson, B. J., J. Yang, H. K. Finucane, A. Gusev, S. Lindström, S. Ripke, G. Genovese, P.-R. Loh, G. Bhatia, R. Do, T. Hayeck, H.-H. Won, Schizophrenia Working Group of the Psychiatric Genomics Consortium, DRIVE study, S. Kathiresan, M. Pato, C. Pato, R. Tamimi, E. Stahl, N. Zaitlen, B. Pasaniuc, G. Belbin, E. E. Kenny, M. H. Schierup, P. De Jager, N. A. Patsopoulos, S. McCarroll, M. Daly, S. Purcell, D. Chasman, B. Neale, M. Goddard, P. M. Visscher, P. Kraft, N. Patterson, and A. L. Price (2015).  Modeling linkage disequilibrium increases accuracy of polygenic risk scores.  The American Journal of Human Genetics 97(4), 576 – 592\.
  * Walker et al. (2019) Walker, V. M., N. M. Davies, G. Hemani, J. Zheng, P. C. Haycock, T. R. Gaunt, G. Davey Smith, and R. M. Martin (2019).  Using the MR-Base platform to investigate risk factors and drug targets for thousands of phenotypes.  Wellcome Open Research.
  * Waterworth et al. (2010) Waterworth, D. M., S. L. Ricketts, K. Song, L. Chen, J. H. Zhao, S. Ripatti, and Other Authors (2010).  Genetic variants influencing circulating lipid levels and risk of coronary artery disease.  Arteriosclerosis, Thrombosis, and Vascular Biology 30(11), 2264–2276.
  * Yang et al. (2012) Yang, J., T. Ferreira, A. P. Morris, S. E. Medland, G. I. of ANthropometric Traits (GIANT) Consortium, D. G. Replication, M. analysis (DIAGRAM) Consortium, P. A. F. Madden, A. C. Heath, N. G. Martin, G. W. Montgomery, M. N. Weedon, R. J. Loos, T. M. Frayling, M. I. McCarthy, J. N. Hirschhorn, M. E. Goddard, and P. M. Visscher (2012).  Conditional and joint multiple-SNP analysis of GWAS summary statistics identifies additional variants influencing complex traits.  Nature Genetics 44(4), 369–375.
  * Yavorska and Burgess (2017) Yavorska, O. O. and S. Burgess (2017).  MendelianRandomization: an R package for performing Mendelian randomization analyses using summarized data.  International Journal of Epidemiology 46(6), 1734–1739.
  * Zheng et al. (2019) Zheng, J., V. Haberland, D. Baird, V. Walker, P. Haycock, A. Gutteridge, T. G. Richardson, J. Staley, B. Elsworth, S. Burgess, B. B. Sun, J. Danesh, H. Runz, J. C. Maranville, H. M. Martin, J. Yarmolinsky, C. Laurin, M. V. Holmes, J. Liu, K. Estrada, L. McCarthy, M. Hurle, D. Waterworth, M. R. Nelson, A. S. Butterworth, G. D. Smith, G. Hemani, R. A. Scott, and T. R. Gaunt (2019).  Phenome-wide Mendelian randomization mapping the influence of the plasma proteome on complex diseases.  bioRxiv.

## Supplementary Material

This file contains Supplementary material for the manuscript “Statistical
Methods for cis-Mendelian randomizaton”.

### Proportion of Genetic Variation Explained

To fix the proportion $v_{G}$ of variation in the risk factor explained by the
genetic variants, we adjusted the residual variance $\sigma_{0}^{2}$ in our
simulations accordingly. We first specified a value for $v_{G}$ and generated
SNP-risk factor associations $\beta_{Xj}$ as described in the simulation
design section of the manuscript. The value of $\sigma_{0}^{2}$ was then set
equal to

$\sigma_{0}^{2}=\frac{1-v_{G}}{2v_{G}}\frac{1}{N_{ref}}\beta_{X}^{T}G_{ref}^{T}G_{ref}\beta_{X}$

where $G_{ref}$ is the reference genetic matrix. This can be justified by
recalling that in our simulation design,

$X=\sum_{j=1}^{P}\beta_{Xj}G_{j}+\alpha_{X}U+\epsilon_{X}$

with $U,\epsilon_{X}\sim N(0,\sigma_{0}^{2})$. Taking variances, we have that
$\text{Var}(X)=\text{Var}(G\beta_{X})+2\sigma_{0}^{2}$ and therefore,

$v_{G}=\frac{\text{Var}(G\beta_{X})}{\text{Var}(X)}=\frac{\text{Var}(G\beta_{X})}{\text{Var}(G\beta_{X})+2\sigma_{0}^{2}}=\frac{\frac{1}{N_{ref}}\beta_{X}^{T}G^{T}G\beta_{X}}{\frac{1}{N_{ref}}\beta_{X}^{T}G^{T}G\beta_{X}+2\sigma_{0}^{2}}$

which yields the previous expression for $\sigma_{0}^{2}$. A similar formula
was derived in the appendix of Yang et al. (2012).

### Additional Simulations - Large G-X Sample Size

To augment the analysis presented in the paper, we compared the various cis-MR
methods in a range of additional simulations. This section contains simulation
results for a simulation using a larger sample size of $N_{1}=100000$ to
compute summary-level risk factor-outcome associations. This scenario may
represent an applied analysis using a downstream biomarker as a proxy for
protein expression, as this would allow researchers to obtain summary-level
data from existing large-scale GWAS studies whose sample size often ranges in
the hundreds of thousands. This was the case for the two real-data
applications presented in the main part of our manuscript.

We conducted simulations both for the SHBG and for the HMGCR region, with six
causal variants per region, as in our main simulations. For brevity, we only
considered the “strong instruments” scenario, where $v_{G}=3\%$ for the SHBG
region and $v_{G}=2\%$ for the HMGCR region. Otherwise, the simulations
reported here were set up in the same way as those reported in the main part
of our paper. Results are reported in Table 5.

|  | $\theta=0$ | $\theta=0.05$ | $\theta=0.1$
---|---|---|---|---
Method | $\hat{\theta}$ | $se(\hat{\theta})$ | Type I | $\hat{\theta}$ | $se(\hat{\theta})$ | Cov | Power | $\hat{\theta}$ | $se(\hat{\theta})$ | Cov
SHBG Region
Top SNP | —– | -0.001 | 0.019 | 0.063 | 0.050 | 0.020 | 0.951 | 0.710 | 0.101 | 0.020 | 0.952
Pruning | $\rho=0.1$ | -0.001 | 0.016 | 0.054 | 0.050 | 0.016 | 0.943 | 0.868 | 0.100 | 0.016 | 0.954
$\rho=0.3$ | -0.001 | 0.014 | 0.039 | 0.050 | 0.014 | 0.946 | 0.940 | 0.099 | 0.015 | 0.958
$\rho=0.5$ | -0.001 | 0.014 | 0.041 | 0.050 | 0.014 | 0.946 | 0.949 | 0.099 | 0.014 | 0.956
$\rho=0.7$ | -0.001 | 0.013 | 0.033 | 0.049 | 0.014 | 0.947 | 0.944 | 0.098 | 0.014 | 0.952
$\rho=0.9$ | -0.001 | 0.013 | 0.036 | 0.048 | 0.013 | 0.949 | 0.941 | 0.095 | 0.014 | 0.941
PCA | $k=0.99$ | 0.000 | 0.015 | 0.038 | 0.051 | 0.015 | 0.950 | 0.909 | 0.100 | 0.016 | 0.957
$k=0.999$ | 0.000 | 0.014 | 0.040 | 0.051 | 0.015 | 0.944 | 0.939 | 0.100 | 0.015 | 0.959
JAM | $\rho=0.6$ | -0.001 | 0.014 | 0.037 | 0.050 | 0.014 | 0.947 | 0.951 | 0.100 | 0.014 | 0.959
$\rho=0.8$ | -0.001 | 0.013 | 0.038 | 0.051 | 0.014 | 0.949 | 0.948 | 0.100 | 0.014 | 0.953
$\rho=0.9$ | 0.000 | 0.013 | 0.040 | 0.051 | 0.014 | 0.947 | 0.952 | 0.100 | 0.014 | 0.960
$\rho=0.95$ | 0.000 | 0.013 | 0.038 | 0.051 | 0.014 | 0.952 | 0.951 | 0.100 | 0.014 | 0.960
F-LIML | —– | -0.001 | 0.014 | 0.046 | 0.051 | 0.015 | 0.952 | 0.919 | 0.100 | 0.015 | 0.961
CLR | —– | —– | —– | 0.044 | —– | —– | 0.949 | 0.919 | —– | —– | 0.958
HMGCR Region
Top SNP | —– | 0.000 | 0.019 | 0.055 | 0.050 | 0.019 | 0.950 | 0.738 | 0.100 | 0.020 | 0.951
Pruning | $\rho=0.1$ | 0.000 | 0.018 | 0.048 | 0.050 | 0.018 | 0.956 | 0.769 | 0.100 | 0.019 | 0.956
$\rho=0.3$ | 0.000 | 0.017 | 0.051 | 0.049 | 0.018 | 0.947 | 0.795 | 0.099 | 0.018 | 0.951
$\rho=0.5$ | 0.000 | 0.017 | 0.048 | 0.050 | 0.017 | 0.950 | 0.818 | 0.099 | 0.018 | 0.941
$\rho=0.7$ | 0.000 | 0.016 | 0.048 | 0.049 | 0.017 | 0.947 | 0.836 | 0.098 | 0.017 | 0.942
$\rho=0.9$ | 0.000 | 0.016 | 0.048 | 0.048 | 0.017 | 0.956 | 0.821 | 0.096 | 0.017 | 0.938
PCA | $k=0.99$ | 0.000 | 0.017 | 0.047 | 0.050 | 0.017 | 0.956 | 0.815 | 0.100 | 0.018 | 0.947
$k=0.999$ | 0.000 | 0.017 | 0.049 | 0.050 | 0.017 | 0.956 | 0.819 | 0.100 | 0.018 | 0.948
JAM | $\rho=0.6$ | 0.000 | 0.017 | 0.049 | 0.050 | 0.017 | 0.952 | 0.816 | 0.099 | 0.018 | 0.937
$\rho=0.8$ | 0.000 | 0.017 | 0.048 | 0.050 | 0.017 | 0.955 | 0.827 | 0.100 | 0.017 | 0.942
$\rho=0.9$ | 0.000 | 0.017 | 0.045 | 0.050 | 0.017 | 0.957 | 0.829 | 0.099 | 0.017 | 0.938
$\rho=0.95$ | 0.000 | 0.017 | 0.049 | 0.050 | 0.017 | 0.951 | 0.825 | 0.099 | 0.018 | 0.942
F-LIML | —– | 0.000 | 0.017 | 0.049 | 0.051 | 0.017 | 0.952 | 0.825 | 0.100 | 0.018 | 0.949
CLR | —– | —– | —– | 0.049 | —– | —– | 0.947 | 0.823 | —– | —– | 0.947
Table 5: Simulation results for the “strong instruments” scenario and a large
$G-X$ sample size ($N_{1}=100000$).

A notable difference in the results compared to our main simulations was that
a larger sample size resulted in stronger instruments. A tenfold increase in
the sample size resulted in a similar, tenfold increase in the values of the F
statistics. For a regression using only the six causal variants, the average F
statistic was $514$ for the SHBG and $340$ for the HMGCR region. As a result,
the various cis-MR methods were less affected by weak instruments bias and
performed quite well. The performance of the stepwise pruning method was was
much more consistent for different values of $\rho$ than in the corresponding
simulation with $N_{1}=10000$. PCA and JAM were both unaffected by weak
instrument bias for $\theta=0.1$. Confidence intervals constructed using these
methods practically attained nominal coverage, and so did the confidence
intervals based on F-LIML and CLR. A small reduction in power for single-SNP
analysis and LD-pruning with low correlation thresholds were the only real
difference between the various methods. This was the case both for the SHBG
and for the HMGCR region.

These results imply that weak instrument bias is less likely to affect
analyses based on large $G-X$ sample sizes, and stepwise pruning does not
underperform other methods in such analyses.

### Additional Simulations - Fewer Causal Variants

In our second set of additional simulations, we modified the number of causal
variants that were present in each region. Focusing on the SHBG region, we
considered two additional scenarios.

First, we assumed the existence of only a single causal variant in the region,
placing the causal signal at the variant that had the smallest univariate
p-value in our real-data application (rs1799941). That variant was assumed to
have an effect of $\beta_{Xj}=0.38$ on the risk factor, the same as that
observed in our real SHBG-testosterone dataset.

In the second scenario, we generated three independent genetic effects on the
risk factor. The effects were placed at genetic variants suggested as
independently causal by Jin et al. (2012), who analysed genetic associations
of variants in the SHBG region with serum testosterone levels. Genetic effects
for the causal variants were drawn randomly according to
$\beta_{Xj}\sim|N(0,0.2)|+0.1$ for the risk-increasing allele, same as what we
did for the six-causal-variants simulation in the main part of our paper.

Otherwise, the simulations were set up as previously described. We used the
“strong instruments” scenario with $v_{G}=3\%$ and a $G-X$ sample size of
$N_{1}=10000$. Simulation results are reported in Table 6.

|  | $\theta=0$ | $\theta=0.05$ | $\theta=0.1$
---|---|---|---|---
Method | $\hat{\theta}$ | $se(\hat{\theta})$ | Type I | $\hat{\theta}$ | $se(\hat{\theta})$ | Cov | Power | $\hat{\theta}$ | $se(\hat{\theta})$ | Cov
One Causal Variant
Top SNP | —– | 0.001 | 0.013 | 0.046 | 0.051 | 0.014 | 0.952 | 0.966 | 0.100 | 0.014 | 0.937
Pruning | $\rho=0.1$ | 0.001 | 0.013 | 0.046 | 0.051 | 0.014 | 0.952 | 0.966 | 0.100 | 0.014 | 0.937
$\rho=0.3$ | 0.001 | 0.013 | 0.047 | 0.050 | 0.014 | 0.953 | 0.961 | 0.099 | 0.014 | 0.932
$\rho=0.5$ | 0.001 | 0.013 | 0.045 | 0.049 | 0.014 | 0.949 | 0.956 | 0.096 | 0.014 | 0.920
$\rho=0.7$ | 0.001 | 0.013 | 0.047 | 0.047 | 0.013 | 0.947 | 0.959 | 0.093 | 0.014 | 0.899
$\rho=0.9$ | 0.001 | 0.012 | 0.050 | 0.043 | 0.013 | 0.918 | 0.934 | 0.086 | 0.013 | 0.767
PCA | $k=0.99$ | 0.000 | 0.015 | 0.053 | 0.050 | 0.016 | 0.946 | 0.882 | 0.099 | 0.016 | 0.924
$k=0.999$ | 0.001 | 0.014 | 0.051 | 0.048 | 0.014 | 0.952 | 0.922 | 0.096 | 0.015 | 0.919
JAM | $\rho=0.6$ | 0.001 | 0.013 | 0.046 | 0.051 | 0.014 | 0.952 | 0.965 | 0.100 | 0.014 | 0.938
$\rho=0.8$ | 0.001 | 0.013 | 0.046 | 0.051 | 0.014 | 0.952 | 0.966 | 0.100 | 0.014 | 0.938
$\rho=0.9$ | 0.001 | 0.013 | 0.046 | 0.051 | 0.014 | 0.953 | 0.966 | 0.100 | 0.014 | 0.938
$\rho=0.95$ | 0.001 | 0.013 | 0.046 | 0.051 | 0.014 | 0.953 | 0.966 | 0.100 | 0.014 | 0.938
F-LIML | —– | 0.001 | 0.014 | 0.061 | 0.051 | 0.015 | 0.944 | 0.936 | 0.100 | 0.016 | 0.942
CLR | —– | —– | —– | 0.047 | —– | —– | 0.951 | 0.928 | —– | —– | 0.955
Three Causal Variants
Top SNP | —– | 0.000 | 0.016 | 0.051 | 0.050 | 0.016 | 0.956 | 0.852 | 0.098 | 0.017 | 0.934
Pruning | $\rho=0.1$ | 0.000 | 0.015 | 0.047 | 0.049 | 0.016 | 0.953 | 0.868 | 0.097 | 0.016 | 0.928
$\rho=0.3$ | 0.000 | 0.014 | 0.046 | 0.049 | 0.014 | 0.942 | 0.919 | 0.096 | 0.015 | 0.930
$\rho=0.5$ | 0.000 | 0.013 | 0.049 | 0.047 | 0.014 | 0.945 | 0.928 | 0.095 | 0.014 | 0.915
$\rho=0.7$ | 0.000 | 0.013 | 0.046 | 0.046 | 0.013 | 0.931 | 0.934 | 0.091 | 0.014 | 0.867
$\rho=0.9$ | 0.000 | 0.012 | 0.042 | 0.041 | 0.013 | 0.891 | 0.904 | 0.083 | 0.013 | 0.724
PCA | $k=0.99$ | 0.000 | 0.015 | 0.051 | 0.049 | 0.015 | 0.941 | 0.884 | 0.096 | 0.016 | 0.932
$k=0.999$ | 0.000 | 0.014 | 0.047 | 0.047 | 0.014 | 0.944 | 0.905 | 0.095 | 0.015 | 0.925
JAM | $\rho=0.6$ | 0.000 | 0.014 | 0.046 | 0.049 | 0.014 | 0.945 | 0.939 | 0.097 | 0.015 | 0.930
$\rho=0.8$ | 0.000 | 0.014 | 0.046 | 0.049 | 0.014 | 0.946 | 0.938 | 0.097 | 0.014 | 0.929
$\rho=0.9$ | 0.000 | 0.014 | 0.046 | 0.049 | 0.014 | 0.942 | 0.943 | 0.097 | 0.014 | 0.923
$\rho=0.95$ | 0.000 | 0.013 | 0.044 | 0.049 | 0.014 | 0.946 | 0.948 | 0.097 | 0.014 | 0.929
F-LIML | —– | 0.000 | 0.014 | 0.060 | 0.050 | 0.014 | 0.932 | 0.933 | 0.100 | 0.016 | 0.949
CLR | —– | —– | —– | 0.049 | —– | —– | 0.940 | 0.919 | —– | —– | 0.959
Table 6: Simulation results for the SHBG region with either 1 or 3 causal
variants.

The results closely resembled the ones for our original simulation scenario
with “strong instruments” and six causal variants. All methods were quite
accurate for $\theta=0$. When $\theta\neq 0$, stepwise pruning was subject to
weak instrument bias to some extent, especially for large correlation
thresholds. The performance of the method depended on the correlation
threshold used. PCA and JAM were more consistent in terms of their tuning
parameters and were generally quite accurate, while factor-based methods were
even more accurate, with a small inflation of Type I error rates under the
null for F-LIML.

Most methods performed slightly better with one than with three causal
variants, but the differences in their performance were small and it seems
that the number of causal signals in the region should have little impact on
the choice of which cis-MR method to use. An exception relates to the use of
the top-SNP approach, which was expectedly quite accurate when only a single
causal variant existed in the region. With three causal variants, top-SNP
analysis was subject to the same issues as in our original simulations (namely
larger standard errors and lower power than other methods), but to a lesser
degree.

### Additional Simulations - Small Reference Dataset

In our third set of additional simulations, we wanted to assess the impact of
the reference dataset on the various MR methods. In particular, we assessed
whether the use of a small reference dataset, which would imply inaccurate
genetic correlation estimates, affects some cis-MR methods worse than others.
We therefore repeated the simulations from the main part of our paper, but
bootstrapped the rows of the UK Biobank matrix and only selected
$N_{ref}=1000$ individuals from which to compute genetic correlations, instead
of using the entire UK Biobank dataset of $N_{bb}=367643$ individuals. In
practice, using small datasets as reference data is not uncommon: for example,
the 1000 genomes dataset is commonly used for this purpose. The use of a small
reference dataset does not systematically bias the estimated genetic
covariance matrix, as would be the case with population stratification.

We used both genetic regions and the “strong instruments” scenario. Again, the
simulation design was identical to those reported in the main part of the
paper, except using a smaller reference dataset. Results are reported in Table
7.

|  | $\theta=0$ | $\theta=0.05$ | $\theta=0.1$
---|---|---|---|---
Method | $\hat{\theta}$ | $se(\hat{\theta})$ | Type I | $\hat{\theta}$ | $se(\hat{\theta})$ | Cov | Power | $\hat{\theta}$ | $se(\hat{\theta})$ | Cov
SHBG Region
Top SNP | —– | -0.001 | 0.019 | 0.039 | 0.048 | 0.019 | 0.941 | 0.676 | 0.095 | 0.020 | 0.917
Pruning | $\rho=0.1$ | -0.001 | 0.016 | 0.046 | 0.048 | 0.016 | 0.925 | 0.799 | 0.095 | 0.017 | 0.892
$\rho=0.3$ | 0.000 | 0.014 | 0.045 | 0.048 | 0.014 | 0.942 | 0.900 | 0.095 | 0.015 | 0.901
$\rho=0.5$ | 0.000 | 0.013 | 0.049 | 0.047 | 0.014 | 0.933 | 0.914 | 0.093 | 0.014 | 0.892
$\rho=0.7$ | 0.000 | 0.013 | 0.048 | 0.046 | 0.013 | 0.917 | 0.916 | 0.090 | 0.014 | 0.863
$\rho=0.9$ | 0.000 | 0.012 | 0.057 | 0.043 | 0.012 | 0.871 | 0.900 | 0.084 | 0.013 | 0.750
PCA | $k=0.99$ | 0.000 | 0.015 | 0.041 | 0.049 | 0.015 | 0.939 | 0.883 | 0.097 | 0.015 | 0.922
$k=0.999$ | 0.000 | 0.014 | 0.052 | 0.048 | 0.014 | 0.925 | 0.904 | 0.094 | 0.015 | 0.919
JAM | $\rho=0.6$ | 0.000 | 0.014 | 0.038 | 0.047 | 0.014 | 0.937 | 0.896 | 0.093 | 0.015 | 0.906
$\rho=0.8$ | 0.000 | 0.014 | 0.039 | 0.047 | 0.014 | 0.922 | 0.896 | 0.092 | 0.014 | 0.887
$\rho=0.9$ | 0.000 | 0.014 | 0.063 | 0.047 | 0.014 | 0.910 | 0.895 | 0.092 | 0.014 | 0.867
$\rho=0.95$ | 0.000 | 0.014 | 0.093 | 0.047 | 0.014 | 0.888 | 0.885 | 0.091 | 0.014 | 0.843
F-LIML | —– | 0.000 | 0.014 | 0.052 | 0.050 | 0.014 | 0.921 | 0.907 | 0.099 | 0.016 | 0.929
CLR | —– | —– | —– | 0.040 | —– | —– | 0.929 | 0.895 | —– | —– | 0.933
HMGCR Region
Top SNP | —– | -0.001 | 0.018 | 0.047 | 0.048 | 0.018 | 0.951 | 0.753 | 0.099 | 0.019 | 0.921
Pruning | $\rho=0.1$ | -0.001 | 0.018 | 0.048 | 0.048 | 0.018 | 0.951 | 0.757 | 0.099 | 0.019 | 0.919
$\rho=0.3$ | -0.001 | 0.017 | 0.056 | 0.047 | 0.017 | 0.951 | 0.770 | 0.097 | 0.018 | 0.928
$\rho=0.5$ | -0.001 | 0.016 | 0.052 | 0.047 | 0.017 | 0.947 | 0.807 | 0.096 | 0.017 | 0.920
$\rho=0.7$ | -0.001 | 0.016 | 0.052 | 0.046 | 0.016 | 0.943 | 0.818 | 0.093 | 0.017 | 0.893
$\rho=0.9$ | 0.000 | 0.014 | 0.113 | 0.042 | 0.015 | 0.836 | 0.785 | 0.084 | 0.015 | 0.735
PCA | $k=0.99$ | -0.001 | 0.017 | 0.047 | 0.048 | 0.017 | 0.950 | 0.816 | 0.098 | 0.018 | 0.921
$k=0.999$ | 0.000 | 0.016 | 0.051 | 0.047 | 0.017 | 0.949 | 0.814 | 0.097 | 0.017 | 0.922
JAM | $\rho=0.6$ | -0.001 | 0.017 | 0.047 | 0.047 | 0.018 | 0.949 | 0.782 | 0.096 | 0.018 | 0.933
$\rho=0.8$ | -0.001 | 0.017 | 0.062 | 0.046 | 0.017 | 0.934 | 0.794 | 0.094 | 0.018 | 0.910
$\rho=0.9$ | 0.000 | 0.017 | 0.093 | 0.045 | 0.017 | 0.895 | 0.790 | 0.093 | 0.018 | 0.859
$\rho=0.95$ | -0.001 | 0.017 | 0.150 | 0.044 | 0.017 | 0.824 | 0.759 | 0.092 | 0.018 | 0.791
F-LIML | —– | -0.001 | 0.016 | 0.052 | 0.049 | 0.017 | 0.941 | 0.825 | 0.101 | 0.019 | 0.945
CLR | —– | —– | —– | 0.049 | —– | —– | 0.945 | 0.810 | —– | —– | 0.942
Table 7: Simulation results for the “strong instruments” scenario, with a
small reference sample ($N_{ref}=1000$).

The JAM algorithm proved to be more sensitive to the reference dataset. The
algorithm adjusts marginal SNP-trait associations using the reference genetic
correlations. By using an inaccurate correlation pattern, the algorithm’s
adjustment was adversely affected and this resulted in selecting more genetic
variants than in runs with a larger reference dataset. This was especially the
case for large values of the correlation threshold, in which case the pre-
pruning step only discards a small number of variants before running JAM. For
example, for the SHBG region with $\theta=0$ and $\rho=0.95$, the algorithm
had a posterior model size of $9.27$ compared to $3.75$ when the entire UK
Biobank was used as a reference dataset. This meant that the algorithm was
more susceptible to weak instrument bias, made JAM causal effect estimates
slightly more variable for larger $\rho$ values and increased Type I error
rates. For smaller values of the correlation threshold the algorithm’s
performance was affected less.

The performance of LD-pruning also attenuated for large values of its
correlation threshold, but this attenuation was in line with what we observed
in our original simulations, using a large reference dataset. Likewise,
principal components analysis and factor-based methods seemed fairly robust to
using a smaller reference dataset.