run_sim = function() {
N = 1000
P = 1
beta = rep(0.5, P)
alpha = 0.5
U_sd = 0.3
G = rnorm(N * P, mean = 0, sd = 1)
U = rnorm(N, mean = 0, sd = U_sd)
X = G * beta + rnorm(N, mean = 0, sd = sqrt(1.0 - beta ^ 2 - U_sd ^ 2)) + U
Y = X * alpha + rnorm(N, mean = 0, sd = sqrt(1.0 - alpha ^ 2 * var(X) - U_sd ^ 2)) + U
return(coef(lm(Y ~ X))[2])
}
sims = map_dbl(1:1000, ~run_sim())GWAS IV: Mendelian Randomization
BIOS 770
Josh Weinstock, PhD
Assistant Professor, Department of Human Genetics, Emory University
Learning objectives
- Discuss the assumptions of Mendelian randomization (MR)
- Discuss a few notable MR methods
- Close with discussion on the MR literature broadly
Inferring causal effects of X on Y
Consider two phenotypes, measured in humans:
- X, some “exposure”, e.g., low-density lipoprotein, C-reactive protein, etc.
- Y, some “outcome”, e.g., ischemic heart disease, systolic blood pressure, etc.
Inferring the causal effect of X on Y
We are interested in inferring the causal effect of X on Y. Perhaps ideally, we could randomly assign values to X, i.e., do(X=x), and then measure how Y changes as we change do(X=x).
This works (e.g, RCTs). But obviously, for numerous pairs of (X,Y), RCT isn’t tractable.
Transmission of alleles as a pseudo-RCT
Perhaps we can instead use natural variation, with the right kind of randomness, as a proxy for an RCT, a kind of “natural experiment.” Consider genetic inheritance:
- In the autosomes, offspring randomly inherit one of mother’s two alleles, and one of the father’s two alleles.
- These alleles are transmitted indepenently on one another (putting aside linkage disequilibrium for now)
Perhaps with these principles in mind, we can view each trait relevant variant as a kind of RCT.
Genetic variation as a parallel to an RCT
MR study design, Sanderson et al., 2022, Nature reviews methods primers
Mendelian Randomization
To do this, we can use Mendelian Randomization (MR), which leverages genetic variants as instrumental variables (IVs) to estimate the causal effect of an exposure (X) on an outcome (Y).
In the case of no confounding between X and Y:
Xi=μx+Giβ+ϵxi
Yi=μy+Xiα+ϵyi
Assuming E(ϵx)=E(ϵy)=0 and Cov(ϵx,ϵy)=0. Gi are genotypes for the ith individual.
Implications of this two stage model: Expectation
E(Yi|Gi)=E(E(Yi|Xi,Gi))=μy+αE(Xi|Gi)=μy+α(Giβ+μx)
Implications of this two stage model: Covariance
Cov(Y,X∣G)=Cov(μy+α∗X+ϵy, X∣G)=α∗Var(X∣G)
Cov(Y,X)=E(Cov(Y,X∣G))+Cov(E(Y∣G),E(X∣G))=α∗E(Var(X∣G))+α∗Cov(E(X∣G),E(X∣G))=α∗Var(ϵx)+α∗β2Var(G)
Note: Only true if Cov(ϵy,X)=0 !
Now, a bit more realistic, with a confounder
However, if there was no confounding between X and Y, why do we need MR anyway?
Now, a bit more realistic, with a confounder
Consider the following scenario where there is a confounder U that affects both X and Y:
Xi=μx+Giβ+Uiγ+ϵxiYi=μy+Xiα+Uiδ+ϵyi
Here, Ui represents the confounder, γ is the effect of the confounder on X, and δ is the effect of the confounder on Y.
Implications of this model with a confounder
The presence of the confounder U introduces bias in the estimation of the causal effect α.
Covariance
Cov(Y,X∣G)=Cov(μy+αX+Uδ+ϵy, X∣G)=α∗Var(X∣G)+δ∗Cov(U,X∣G)
Note: The term δ∗Cov(U,X∣G) introduces bias due to the confounder U.
Addressing confounding in MR
To address confounding, we can use genetic variants as instrumental variables (IVs) that are not associated with the confounder U. This allows us to estimate the causal effect α without bias from U.
Example
Suppose we have identified a set of genetic variants that are associated with body mass index (BMI) and we want to estimate the causal effect of BMI on blood pressure. We can use these genetic variants as instruments in an MR analysis to infer the causal relationship between BMI and blood pressure.
Example
Sanderson et al., 2022, Nature reviews methods primers
Example confounders
- Age
- Smoking
- Diet
- Alcohol consumption
- Comorbid diseases
- Many others…
Assumptions of MR
- Relevance: (G is associated with X)
- Independence: (G is not associated with U)
- Exclusion restriction: The total effect of G on Y is mediated through X
How to estimate α
Consider the least squares estimate of the effect of X on Y
ˆα=(XtX)−1XtY
What if we didn’t use the observed X, but a proxy for it using G ?
Estimation
ˆβ=(GtG)−1GtXˆX=G(GtG)−1GtX
We can estimate α by plugging in our estimate of X
ˆα=((G(GtG)−1GtX)t(G(GtG)−1GtX))−1(G(GtG)−1GtX)tYPG=G(GtG)−1Gˆα=((PGX)t(PGX))−1(PGX)tYˆα=(XtPGX)−1XtPGY
Estimation
E((XtPGX)−1XtPGY)
=E(E((XtPGX)−1XtPGY∣X))=E((XtPGX)−1(XtPGX)α+(XtPGX)−1E(ϵy|X))=α∗I+E((XtPGX)−1E(ϵy|X))=αif A1
A1: E(Y|X,G)=αX
In general A1 is not true in finite samples
The asymptotic bias
To understand the asymptotic bias of the IV estimator, consider the bias:
E((XtPGX)−1XtPGϵy)
- Expectation is difficult to evaluate, because we can’t move the expectation operator in the inverse (a non-linear function)
- However, we can approximate this with some asymptotic arguments.
1nXtPGX→QXG
1nXtPGϵy→0by assumption
Then: E((XtPGX)−1XtPGϵy)→Q−1XG∗0=0
This is an estimator of α if one has access to G, X, and Y all at the same time. Often, however, this is not the case.
Two-sample MR
It’s quite natural in MR to ask whether some molecular phenotype mediates the assocation between genotype and phenotype
- Does gene expression mediate the association between G and Y? (e.g., TWAS)
- Does protein abundance mediate the association? (e.g. PWAS)
- Do metabolites mediate the association?
Often, interesting molecular phenotypes are measured in distinct cohorts (e.g., GTEx, Geuvadis) from the GWAS.
Inverse-variance weighted estimator (IVW)
Assume that summary statistics for the following have been generated:
- Association between G and X : ˆβ, ^σ2x,
- Association between G and Y : ˆψ, ^σ2y
ˆαIVW=∑Pi=1Wiˆψiˆβi∑Pi=1WiwithWi=(ˆβ^σy)2
How would you estimate this the lm command in R?
What if some of our instruments are invalid?
- Perhaps some of the SNPs don’t predict the phenotype
- Perhaps some of the SNPs have direct effects on the outcome phenotype
Examining this behavior through simulation
Now, with instrumental variables
run_sim_iv = function() {
N = 1000
P = 1
beta = 0.5
alpha = 0.5
U_sd = 0.3
G = rnorm(N * P, mean = 0, sd = 1)
U = rnorm(N, mean = 0, sd = U_sd)
X = G * beta + rnorm(N, mean = 0, sd = sqrt(1.0 - beta ^ 2 - U_sd ^ 2)) + U
Y = X * alpha + rnorm(N, mean = 0, sd = sqrt(1.0 - alpha ^ 2 * var(X) - U_sd ^ 2)) + U
YG = coef(lm(Y ~ G))[2]
XG = coef(lm(X ~ G))[2]
return(YG / XG)
}
sims_iv = map_dbl(1:1000, ~run_sim_iv())
Weak instrument simulation
run_sim_iv_weak = function() {
N = 1000
P = 1
beta = 0.03
alpha = 0.5
U_sd = 0.3
G = rnorm(N * P, mean = 0, sd = 1)
U = rnorm(N, mean = 0, sd = U_sd)
X = G * beta + rnorm(N, mean = 0, sd = sqrt(1.0 - beta ^ 2 - U_sd ^ 2)) + U
Y = X * alpha + rnorm(N, mean = 0, sd = sqrt(1.0 - alpha ^ 2 * var(X) - U_sd ^ 2)) + U
YG = coef(lm(Y ~ G))[2]
XG = coef(lm(X ~ G))[2]
return(YG / XG)
}
set.seed(3)
sims_iv_weak = map_dbl(1:1000, ~run_sim_iv_weak())Bias when instruments are weak
Median estimator
- Let ˆαj denote the jth ordered estimate.
- Simply take the median of these values.
- Can be improved by taking the inverse of the variance of the estimates as weights.
- Consistent even if 50% of the instruments are invalid.
- (See Bowden et al., 2016, Genetic Epidemiology).
MR-RAPS
Consider the likelihood
l(α,β1,…,βP)=−0.5[P∑i=1(ˆβi−βi)2σ2xi+P∑i=1(ˆψi−βiα)2σ2yi]
After profiling out the β terms:
l(α)=−0.5P∑i=1(ˆψi−αˆβi)2σ2Xiα2+σ2Yi
Zhao et al., 2020, Annals of Statistics.
Robust likelihood
l(α)=−0.5P∑i=1ρ(ˆψi−αˆβi)√σ2Xiα2+σ2Yi
with ρ defined as some “robust” loss function, e.g., the Huber loss:
ρ(r;k)={12r2,if |r|≤k,k|r|−12k2,if |r|>k.
Returning to pleiotropy
Zhao et al., 2020, AOS
Complex traits are highly pleiotropic
Table 1, Watanabe et al., Nature Genetics, 2019
Complex traits are highly polygenic
Figure 1, Zhang et al., Nature Communications, 2022
Correlated and Uncorrelated pleiotropy
Figure 1, Morrison et al., Nature Genetics, 2020
MR-CAUSE
Methods, Morrison et al., Nature Genetics, 2020
Pleiotropy in MR-RAPS
Assume: ψ−α∗β∼N(0,τ2)
Then:
l(α)=−0.5∗P∑i=1(ˆψi−αˆβ)2σ2Xiα2+σ2Yi+ˆτ2
Sparsity assumptions on uncorrelated pleiotropy terms
Berzuini et al., 2020, Biostatistics
Bayesian inference of pleiotropy terms
Consider the following structural equation models: Xi=μx+Giβ+Uiγ+ϵxiYi=μy+Xiα+Uiδ+Giθ+ϵyi
What probabilistic assumptions can we place on θ to make these terms identifiable?
Sparse prior on θ
θj∼N(0,λ2jτ2)λj∼C+(0,1)τ∼C+(0,1)
What is this doing?
Carvalho et al., 2009, JMLR
Posterior computation performed with Hamiltonian Monte-Carlo implemented in the Stan probabilistic programming language.
MR-PRESSO
Verbanck et al., 2018, Nature Genetics
MR-PRESSO
Figure 1a, Verbanck et al.
MR-PRESSO
Figure 1b, Verbanck et al.
MR-PRESSO
Figure 1c, Verbanck et al.
Network MR: bidirectional MR on many phenotypes
Brown and Knowles, 2020, biorxiv
- Estimate all pairwise total effects between all possible X and Y
- Define all total effects as the matrix T
- Estimate a sparse inverse G=T−1, which estimates direct effects
MR “in action”
Voight et al., 2012, The Lancet
- Single SNP MR to assess HDL→MI
- Missense variant in LIPG, reported to affect HDL but not other lipids
- Examine association between this SNP with MI in many cohorts
Figure 2, Voight et al.
MR literature crisis
MR literature crisis
Stender et al., Lipids in Health and Disease, 2024
Call to editors:
“First, the public availability of summary results from genome-wide association studies has led to an explosion of low-quality two-sample mendelian randomization (2SMR) studies. These studies add minimal – if any – value and overwhelm reviewers and journals. Second, the availability of large datasets with individual-level genotype data, like UK Biobank, has spurred the development and use of novel MR methods. However, some methods are being applied without proper testing, leading to misleading results, as exemplified by recent spurious findings that are being retracted and/or corrected relating to vitamin D. What can editors and peer reviewers do to handle the deluge of 2SMR studies and the premature application of highly complex MR methods? We advise editors to simply reject papers that only report 2SMR findings, with no additional supporting evidence.”
Stender et al., Lipids in Health and Disease, 2024
In summary
- MR is an IV approach for estimating the causal effect of X on Y
- Challenges include
- weak instruments
- extensive pleiotropy (both correlated and uncorrelated)
- how to model several phenotypes
- The “convenience” of 2SMR studies has led to a deluge of weak papers
GWAS IV: Mendelian Randomization BIOS 770 Josh Weinstock, PhD Assistant Professor, Department of Human Genetics, Emory University