---
title: "Two-Sample MR Practical"
author: "Caroline Brito Nunes"
date: "2026-07-01"
output:
  html_document:
    toc: true
    toc_float:
      collapsed: false
      smooth_scroll: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

#### With thanks to Kaitlin Wade and Chris Flatley.

```{r}
rm(list=ls())
```

# Preliminary steps

If you haven't already installed the necessary packages with libraries, please do so!

```{r, eval=FALSE}
install.packages("metafor")
install.packages("plyr")
install.packages("meta")
install.packages("rmeta")
```

```{r, results='hide',  message=FALSE}
library(metafor)
library(plyr)
library(meta)
library(rmeta)
```

Set working directory to your local directory.

For example:

```{r}
setwd("/Users/uqcbrito/GGWS26/module4/Practice_2_TSMR")
```

# Practical overview<br>

In this practical we will be testing whether an exposure (here, **BMI**) has a *causal* effect on an outcome (here, **coronary heart disease, CHD**).

The basic logic of a **two-sample MR** (the design used in this practical) is:

1. Find SNPs that are strongly and independently associated with the exposure (BMI) in a large GWAS (here, GIANT/Locke et al. 2015).
2. Look up the association of those *same* SNPs with the outcome (CHD) in an independent GWAS (here, CARDIOGRAM).
3. Make sure the effect estimates in both datasets refer to the *same* effect allele ("harmonization").
4. For each SNP, calculate a **Wald ratio**: the SNP-outcome effect divided by the SNP-exposure effect.
   This estimates the causal effect of the exposure on the outcome.
5. Combine the SNP-specific ratios into a single overall causal estimate (meta-analysis), and run sensitivity
   analyses (e.g. MR-Egger, weighted median, weighted mode) to test the robustness of the causal inference.

## Objectives<br>

+ **Identify independent SNPs from BMI GWAS for use as the instrument variable of MR.**<br><br>

+ **Merge and harmonize with SNPs from the CHD GWAS.**<br><br>

+ **Check for palindromic SNPs and for SNPs in opposing directions.**<br><br>

+ **Estimate Wald Ratio and meta analyze results.**<br><br>

+ **Calculate heterogeneity statistics.**<br><br>

+ **Run sensitivity analyses.**<br><br>


# PART 1: Understanding the data and defining our BMI genetic instruments.<br>

**1. We will be using results from the Locke et al. 2015 paper using data from the GIANT Consortium (downloaded from https://portals.broadinstitute.org/collaboration/giant/index.php/GIANT_consortium_data_files)** <br><br>

Note: If you want to have a look at the full GIANT data, then download, save to your working directory and load in the full results (this is a fairly large file) using the following code:
<br>

giant_full <- read.table("./GIANT_raw.uniq", header = T)
<br>

In the Locke et al. paper, the authors describe a certain number of SNPs that are "approximately independently associated with BMI" across all ancestries.<br>

Let's have a look at these SNPs.

Note: This file was generated by taking the first rows from the Supplementary Table 8 of the Locke et al. paper.

```{r}
snps <- read.csv("./Data/giant_snps_all.csv", header = T)
```
<br>
a. How many SNPs do authors describe as being independently associated with BMI in all ancestries?

```{r}
dim(snps)
head(snps)
```
<br>
b. What information does this file contain that are needed for Mendelian randomization analyses?
```{r}
colnames(snps)
```
<br>
c. How many are associated within only Europeans?
```{r}
table(snps$Sig_Analysis)
```
<br>
d. Why might it be best to use the SNPs that have been identified as being associated with BMI in Europeans only?

<br><br>
**2. Read in the second sheet of this file to get the estimates of the SNPs associated with BMI in Europeans.**
```{r}
euro_snps <- read.csv("./Data/giant_snps_euro.csv", header = T)
dim(euro_snps)
head(euro_snps)
```
Note on units: `Beta` and `SE` are the change in **BMI, in standard deviation (SD) units**, per copy of the
effect allele — Locke et al. inverse-normal transformed BMI before running the GWAS, which is what puts these
estimates on an SD scale rather than raw kg/m². Later, in Part 4, we'll compare these genetic estimates against
an observational BMI-CHD association reported per 4.56 kg/m² (i.e., that study's SD of BMI); we use the same
4.56 kg/m² figure to translate the SD units here into kg/m², as an approximation for interpretability.
<br><br>
a. Check that these are all associated with BMI at a conventional level of genome-wide significance.

```{r}
sort(euro_snps$P)
length(which(euro_snps$P<=5E-8))
```
<br>
b. Are all of these SNPs "good instruments"? What else might we want to check to see if they are strongly and independently associated with BMI?

<br><br>
**3.	We're going to make sure the effect allele is the allele that increases BMI using the effect allele and beta column.
Browse the data, are all SNPs coded so that the effect allele increases BMI?**
<br><br>
a. Are all SNP effects in the same direction?
```{r}
euro_snps[,c("SNP","Effect_Allele","Beta","SE","P")]
summary(euro_snps$Beta)
```
*E.g., a Beta of 0.05 means each copy of that allele raises BMI by 0.05 SD, roughly 0.05 × 4.56 ≈ 0.23 kg/m².*
<br>

# PART 2: SNP lookup in the GWAS of CHD.<br>


**1. These were downloaded from the CARDIOGRAM website http://www.cardiogramplusc4d.org/ and provided for you (however, this is quite a large file so we have truncated for you using the following code):**
<br><br>
CARDIOGRAM <- read.table("./CARDIoGRAM_GWAS_RESULTS.txt", header = T)

BMI_SNPS_in_CARDIOGRAM <- CARDIOGRAM[CARDIOGRAM\$SNP %in% euro_snps$SNP,]

top_1000 <- CARDIOGRAM[c(1:1000),]

bottom_1000 <- CARDIOGRAM[c((nrow(CARDIOGRAM)-1000):(nrow(CARDIOGRAM))),]

CARDIOGRAM_TRUNC <- rbind(BMI_SNPS_in_CARDIOGRAM, top_1000, bottom_1000)

write.table(CARDIOGRAM_TRUNC, "./CARDIOGRAM_CLEANED.txt", quote = F, col.names = T, row.names = F, sep = "\t")
<br><br>
a. How many SNPs are in this truncated CARDIOGRAM dataset?
```{r}
CARDIOGRAM <- read.table("./Data/CARDIOGRAM_CLEANED.txt", sep="\t", header=T, colClasses = "character")
dim(CARDIOGRAM)
head(CARDIOGRAM)
```
Note on units: `log_odds` is the **log odds ratio (log OR) for CHD**, per copy of the reference allele — it's on
the *log* scale, not the OR scale, until it's exponentiated (we'll do that once we get to results in Part 4).
<br><br>
b. Does this file contain everything that is required to perform a two-sample Mendelian randomization analysis?

*You need, at minimum: the SNP identifier, effect (reference) allele, other allele, effect estimate (log OR) and
its standard error, and ideally the effect allele frequency (useful for checking strand alignment and for
calculating instrument strength/F-statistics).*
```{r}
colnames(CARDIOGRAM)
```
<br>
**2. How many of the BMI SNPs are included in the CARDIOGRAM dataset?**
```{r}
BMI_SNPs <- euro_snps$SNP
BMI_SNPs <- as.vector(BMI_SNPs)
matches <- unique(grep(paste(BMI_SNPs, collapse="|"), CARDIOGRAM$SNP, value=TRUE))
matches
```
<br>
View the data from CARDIOGRAM for our BMI SNPs.
```{r}
CARDIOGRAM_BMI <- CARDIOGRAM[grepl(paste(BMI_SNPs, collapse="|"), CARDIOGRAM$SNP),]
CARDIOGRAM_BMI
```
<br>
**3. Merge the GIANT and CARDIOGRAM SNP summary associations.**
<br><br>
First, make sure the column headings are easy to understand (i.e., add "BMI" and "CHD" onto the respective datasets).
```{r}
colnames(euro_snps) <- paste("BMI", colnames(euro_snps), sep = "_")
head(euro_snps)
colnames(CARDIOGRAM) <- paste("CHD", colnames(CARDIOGRAM), sep = "_")
head(CARDIOGRAM)
merged <- merge(euro_snps,CARDIOGRAM, by.x="BMI_SNP", by.y="CHD_SNP")
dim(merged)
head(merged)
```
<br>

# PART 3: Harmonizing the effect alleles in the BMI and CHD datasets.<br>

**1. Make sure that the effect alleles in the CARDIOGRAM and GIANT datasets are the same. We want the CARDIOGRAM effect allele to be the allele that increases BMI.**
<br><br>
But be careful of palindromic SNPs or SNPs on different strands.<br>
First we need to see whether the effect alleles are the same.<br>
Browse the data.
```{r}
merged[,c("BMI_SNP", "BMI_Effect_Allele","CHD_reference_allele","BMI_Other_Allele","CHD_other_allele","BMI_EAF","CHD_ref_allele_frequency")]
```
<br>
a. How can we tell if the CARDIOGRAM and GIANT SNPs are coded using the same reference strand?

*Compare the effect and other allele pairs at each SNP across the two datasets. If, for a non-palindromic SNP,
the two "effect alleles" don't match and neither do the complementary bases, one dataset may be reporting
alleles on the opposite DNA strand.*
<br><br>
b. Are CARDIOGRAM and the GIANT SNPs coded using the same reference strand?
<br><br>
c. Are there any palindromic SNPs?

*Palindromic SNPs (A/T or G/C) are ambiguous because the same two letters are used to describe the alleles
regardless of which DNA strand a dataset happens to report them on — so you cannot tell just from the letters
whether two datasets agree on which allele is which.*
```{r}
palindromic_at<-subset(merged,BMI_Effect_Allele %in% "A" & BMI_Other_Allele %in% "T")
palindromic_ta<-subset(merged,BMI_Effect_Allele %in% "T" & BMI_Other_Allele %in% "A")
palindromic_gc<-subset(merged,BMI_Effect_Allele %in% "G" & BMI_Other_Allele %in% "C")
palindromic_cg<-subset(merged,BMI_Effect_Allele %in% "C" & BMI_Other_Allele %in% "G")
dim(palindromic_at)
dim(palindromic_ta)
dim(palindromic_gc)
dim(palindromic_cg)
```
<br>
d. How can we tell whether the effect alleles are the same in both datasets for palindromic SNPs (i.e., the allele that increases BMI is the same as the reference allele in CARDIOGRAM)?

*For palindromic SNPs, compare the effect allele frequency (EAF) between the two datasets. If a SNP's EAF is
close to 0.5 in either dataset it remains ambiguous and is usually best excluded; otherwise, a similar EAF in
both datasets (e.g., both close to 0.2, or both close to 0.8) implies the alleles are aligned on the same
strand, while a "mirror image" EAF (e.g., 0.2 vs 0.8) implies they are on opposite strands and need flipping.*
<br><br><br>
**2. Make sure the CARDIOGRAM log odds ratio reflects the allele that increases BMI in the GIANT data.**
<br><br>
First, find the positions of SNPs with different effect alleles.
```{r}
head(merged)
effect_diff <- which(merged$BMI_Effect_Allele != merged$CHD_reference_allele) # The position of SNPs where effect alleles are different
merged[effect_diff,c("BMI_SNP", "BMI_Effect_Allele","CHD_reference_allele","BMI_Other_Allele","CHD_other_allele","BMI_EAF","CHD_ref_allele_frequency")]
```
<br>
a. How many SNPs have effect alleles that are coded in the opposite direction?
<br><br>
b. Where the effect alleles are different, flip the direction of the log odds ratio by multiplying it by -1.

*This works because for a log OR, flipping which allele is treated as the "effect" allele simply reverses the
sign of the effect (a doubling in risk for allele X is equivalent to a halving of risk for the other allele).
After this step, `CHD_flip_log_odds` is the log OR for CHD **per copy of the BMI-increasing allele**, matching
the direction of `BMI_Beta`.*
<br><br>
If you want, you can also generate new columns that reflect the effect allele change but this isn't used in the causal estimate.
```{r}
merged$CHD_flip_log_odds <- as.numeric(merged$CHD_log_odds) # Make log odds ratio numeric
merged$CHD_log_odds_se <- as.numeric(merged$CHD_log_odds_se) # Make standard error numeric
merged$CHD_flip_log_odds[effect_diff] <- merged$CHD_flip_log_odds[effect_diff]*(-1)
head(merged)
dim(merged)
```
<br>
c. Check that all of the effect estimates have been flipped appropriately.
```{r}
merged[effect_diff,c("BMI_SNP", "BMI_Effect_Allele","CHD_reference_allele", "CHD_log_odds", "CHD_flip_log_odds")]
merged[-effect_diff,c("BMI_SNP", "BMI_Effect_Allele","CHD_reference_allele", "CHD_log_odds", "CHD_flip_log_odds")]
```
<br>
d. Check that the effect allele frequencies are correlated.

*If the two datasets are describing the same allele at each SNP, their allele frequencies should be strongly,
positively correlated across SNPs. A weak or even negative correlation before harmonization (which improves
after flipping) is a useful confirmation that the harmonization step has worked correctly.*
<br><br>
Check correlation of effect allele frequency between BMI and CARDIOGRAM datasets before harmonising alleles.
```{r}
merged$BMI_EAF <- as.numeric(merged$BMI_EAF)
merged$CHD_ref_allele_frequency <- as.numeric(merged$CHD_ref_allele_frequency)
cor(merged$BMI_EAF,merged$CHD_ref_allele_frequency)
```
<br>
Check correlation of effect allele frequency between BMI and CARDIOGRAM datasets after harmonising alleles
```{r}
merged$CHD_ref_allele_frequency[effect_diff] <- 1-merged$CHD_ref_allele_frequency[effect_diff]
cor(merged$BMI_EAF,merged$CHD_ref_allele_frequency)
```
<br>
e. What happened and why?

*The correlation should increase (move closer to +1) after harmonization. Before flipping, some SNPs had their
frequency reported for the opposite allele, which weakens or reverses the correlation with `BMI_EAF`. Once we
flip the CARDIOGRAM allele frequency for the same SNPs whose log OR we flipped (frequency of the opposite allele
= 1 − original frequency), both datasets are describing the same, BMI-increasing allele, and the frequencies
line up as expected.*
<br><br>

# PART 4: Estimate the effect of BMI on CHD.<br>


**1. Estimate the Wald ratios for each SNP and their delta approximated standard errors.**

*Reminder on units: `gp`/`segp` are in SD-BMI units (1 SD ≈ 4.56 kg/m²); `gd`/`segd` are in log OR units for
CHD. The resulting `wald_ratio` is therefore the **log OR for CHD per 1-SD (≈4.56 kg/m²) increase in BMI** for
each SNP.*
```{r}
gp <- merged$BMI_Beta # The effect of the SNP on BMI (SD-BMI units per effect allele; 1 SD ~ 4.56 kg/m^2)
segp <- merged$BMI_SE # The standard error of the SNP effect on BMI (SD-BMI units)
gd <- merged$CHD_flip_log_odds # The log odds ratio for CHD (that were harmonized to reflect an increase in BMI), log OR per effect allele
segd <- merged$CHD_log_odds_se # Standard error of the log odds ratio (log OR scale)
wald_ratio <- gd/gp # The log odds ratio of CHD per unit (1-SD) change in BMI
Cov <- 0 # Only required when the SNP-BMI and SNP-CHD associations are estimated in the same participants (therefore for two-sample MR with non-overlapping samples, this is set to 0)
wald_ratio_se <- sqrt((segd^2/gp^2) + (gd^2/gp^4)*segp^2 - 2*(gd/gp^3)*Cov) # Delta approximated standard error of the wald ratio; see Thomas, D. C., Lawlor, D. a, & Thompson, J. R. (2007). Re: Estimation of bias in nongenetic observational studies using Mendelian triangulation by Bautista et al. Annals of Epidemiology, 17(7), 5113. doi:10.1016/j.annepidem.2006.12.005
z <- wald_ratio/wald_ratio_se # Z statistic for the wald ratio
p <- 2*pnorm(abs(z) ,lower.tail=F) # P value for the z statistics under the null hypothesis that there is not effect
wald_ratio_var = wald_ratio_se^2 # Variance
weight <- 1/wald_ratio_var # Inverse variance weight
snps <- merged$BMI_SNP # SNPs that we will use in the estimates
```
<br>
**2. Combine the Wald ratios by fixed effects meta-analysis.**

*This produces a single, precision-weighted, overall causal estimate across all instrument SNPs, still on the
log OR per 1-SD BMI scale.*
```{r}
meta_results <- metagen(wald_ratio,wald_ratio_se,comb.fixed=T,sm="OR") #combine the SNPs by fixed effects meta-analysis
meta_results
```
<br>
**3. Create a table of the results, which you could export to other programs e.g. excel, STATA etc.**
```{r}
mr_results <- data.frame(matrix(c(as.character(merged$BMI_SNP),round(wald_ratio,2),round(wald_ratio_se,2),round(p,3)),nrow=length(merged$BMI_SNP),ncol=4))
names(mr_results) <- c("SNP","log_odds","se","p")
mr_results_order <- mr_results[order(mr_results$log_odds),]
overall_genetic_effect <- data.frame(matrix(c("Overall genetic effect",meta_results$TE.fixed,meta_results$seTE.fixed,meta_results$pval.fixed),nrow=1,ncol=4))
names(overall_genetic_effect) <- c("SNP","log_odds","se","p")
overall_genetic_effect
twosampleResults <- rbind.fill(mr_results_order,overall_genetic_effect)
twosampleResults
```
<br>
All `log_odds` values in this table are **log odds ratios for CHD per 1-SD (≈4.56 kg/m²) genetically-predicted
increase in BMI**.

<br>
**4. Calculate heterogeneity p value**

*Heterogeneity here means: do the individual SNP-specific causal estimates (Wald ratios) agree with each other
more than expected by chance alone? High heterogeneity (low p-value) can indicate that one or more SNPs are
affecting CHD through a pathway other than BMI (horizontal pleiotropy), rather than every SNP giving a
consistent estimate of the same true causal effect.*
```{r}
p_het <- pchisq(meta_results$Q,meta_results$df.Q,lower.tail=F)
twosampleResults$p_chi<-NA
```
<br>
Add the p value for heterogeneity to the results table.
```{r}
twosampleResults$p_chi[twosampleResults$SNP=="Overall genetic effect"] <- p_het
write.table(twosampleResults,"./twosample_results_BMI_CHD.txt",sep="\t",col.names=T,row.names=F,quote=F)
twosampleResults
```
<br><br>
**5. Create a forest plot of the results and compare the genetic and observational associations.**
<br><br>
The observational effect is 1.23 (95% CI: 1.17, 1.29) per 4.56 kg/m2 (i.e., per SD) increase in BMI. This is an
**odds ratio for CHD per 1-SD higher (observed, not genetically-predicted) BMI**, estimated in a standard
observational (e.g., cohort/case-control) study.<br>

Formula for SE from 95% confidence interval: (log(uci)-log(lci))/(1.96*2)
```{r, fig.height=10, fig.width=7.5}
effect <- c(wald_ratio[order(wald_ratio)],meta_results$TE.fixed,log(1.23))
se <- c(wald_ratio_se[order(wald_ratio)],meta_results$seTE.fixed,(log(1.29)-log(1.17))/(1.96*2))
snps <- c(as.character(merged$BMI_SNP)[order(wald_ratio)],"Overall genetic effect","Observational effect")
metaplot(effect,se,labels=snps,conf.level=0.95,logeffect=T,nn=0.1,boxsize=0.8,
         xlab="Odds ratio for CHD per 1-SD (~4.56 kg/m^2) higher BMI, with 95% CI",ylab="SNP",cex=0.7)
```
<br>
The plotted boxes and lines are odds ratios and 95% confidence intervals (the `logeffect=T` argument tells
`metaplot` that `effect`/`se` were supplied on the log scale and to exponentiate them for the x-axis); an OR of
1 (vertical reference line) means no effect of BMI on CHD risk, OR > 1 means higher BMI genetically predicts
higher CHD risk, and OR < 1 means the opposite.
<br><br>
**6. Interpret the results.**
<br><br>

**a.** Is the MR-derived effect similar to the observational association?

*Compare the "Overall genetic effect" OR (from the MR analysis, on the per-1-SD-BMI scale) with the
observational OR of 1.23 (95% CI 1.17–1.29) per SD of BMI. Are the point estimates and confidence intervals
similar, or does the MR estimate suggest a smaller/larger/absent causal effect than the observational
association implies? Bear in mind both are expressed in the same units (OR per 1-SD, i.e., per ~4.56 kg/m²,
increase in BMI), so they are directly comparable.*

**b.** Is there evidence of heterogeneity in the genetic effects? How do you interpret this?

*Look at the heterogeneity p-value (`p_chi`) calculated above. A small p-value (e.g., < 0.05) suggests the
SNP-specific Wald ratios are more different from one another than chance alone would predict — a signal that
should prompt closer inspection of individual SNPs (e.g., via the forest plot) and motivates the sensitivity
analyses in Part 5, several of which are specifically designed to still give a valid causal estimate even when
some instruments are pleiotropic.*

**c.** Can you think of reasons for caution?

<br><br>

# PART 5: Sensitivity analyses.<br>


All that is required is summary level results for each SNP (remember gp, segp, gd, segd from PART 4).
<br>

gp <- merged$BMI_Beta **(The effect of the SNP on BMI, in SD-BMI units; 1 SD ≈ 4.56 kg/m²)**

segp <- merged$BMI_SE **(The standard error of the SNP effect on BMI, SD-BMI units)**

gd <- merged$CHD_flip_log_odd **(The log odds ratio for CHD (that were harmonized to reflect an increase in BMI), log OR units)**

segd <- merged$CHD_log_odds_se **(Standard error of the log odds ratio, log OR units)**

<br>
These sensitivity analyses make different assumptions about how pleiotropy (if present) behaves across the
instrument SNPs, and are used together as a "triangulation" exercise: if the IVW, MR-Egger, weighted median and
weighted mode estimates broadly agree in direction and magnitude, that gives more confidence that the overall
causal estimate is not being driven by a small number of pleiotropic SNPs.<br>

**1. These functions define the IVW, MR-Egger, weighted median and weighted mode estimators, respectively, and a function that wraps up the results.**<br>
<br>
set seed for replication purposes

```{r}
set.seed(50)
two.sample.iv.ivw <- function(x, y, sigmax, sigmay) {
  beta.ivw.fit  = summary(lm(y~x-1, weights=sigmay^-2))
  beta.ivw.fit.only  = lm(y~x-1, weights=sigmay^-2)
  beta.ivw      = beta.ivw.fit$coef[1,1]
  beta.se.ivw   = beta.ivw.fit$coef[1,2]/min(beta.ivw.fit$sigma,1)
  beta.df.ivw   = length(y) - 1
  beta.p.ivw    = 2*(1-pt(abs(beta.ivw/beta.se.ivw),beta.df.ivw))
  beta.lower.ivw   = beta.ivw + (-1*qt(df=beta.df.ivw, 0.975)*beta.se.ivw)
  beta.upper.ivw   = beta.ivw + (1*qt(df=beta.df.ivw, 0.975)*beta.se.ivw)
  return(list(beta.ivw=beta.ivw,beta.se.ivw=beta.se.ivw,beta.lower.ivw=beta.lower.ivw,beta.upper.ivw=beta.upper.ivw,beta.t.ivw=beta.ivw/beta.se.ivw,beta.p.ivw=beta.p.ivw,
              beta.ivw.fit.only=beta.ivw.fit.only,beta.df.ivw=beta.df.ivw,beta.ivw.fit=beta.ivw.fit))
}
weighted.median <- function(x, w) {
  N    = length(x)
  ord  = order(x);
  x    = x[ord];
  w    = w[ord];
  Sn   = cumsum(w)
  S_N  = Sn[N]
  Pn   = (100/S_N)*(Sn-w/2)
  if(sort(abs(Pn-50))[1] == 0){M = which(Pn==50); return(x[M])}
  Q  = length(Pn[sign(Pn-50)==-1])
  V1 = Q;  V2 = Q+1
  M  = x[V1] + (50 - Pn[V1])*(x[V2]-x[V1])/(Pn[V2]-Pn[V1])
  return(list(beta.median=M,CumSum.median=Sn,ordX.median=x))
}
weighted.median.boot <- function(x, y, sigmax, sigmay, Nsim, alpha, W) {
  med = NULL
  for (i in 1:Nsim){
    y_boot  = rnorm(length(y), mean=y, sd=sigmay)
    x_boot  = rnorm(length(x), mean=x, sd=sigmax)
    iv_boot = y_boot/x_boot
    run = weighted.median(iv_boot,W)
    med[i]  = run$beta.median
  }
  lower = Nsim*alpha/2
  upper = Nsim*(1-alpha/2)
  Sort  = sort(med)
  lowerCI = Sort[lower]
  upperCI = Sort[upper]
  se    = sd(med)
  t     = mean(med)/se
  p     = 2*(1-pt(abs(t),length(y)-1))

  return(list(beta.se.median=se,beta.lower.median=lowerCI,beta.upper.median=upperCI,beta.t.median=t,beta.p.median=p))
}
two.sample.iv.egger <- function(x, y, sigmax, sigmay) {
  egger.fit       = summary(lm(y~x, weights=sigmay^-2))
  df.egger        = length(y) - 2
  beta.egger      = egger.fit$coef[2,1]
  beta.se.egger   = egger.fit$coef[2,2] / min(egger.fit$sigma, 1)
  beta.p.egger    = 2*(1-pt(abs(beta.egger/beta.se.egger),df.egger))
  beta.lower.egger   = beta.egger + (-1*qt(df=df.egger, 0.975)*beta.se.egger)
  beta.upper.egger   = beta.egger + (1*qt(df=df.egger, 0.975)*beta.se.egger)
  alpha.egger     = egger.fit$coef[1,1]
  alpha.se.egger  = egger.fit$coef[1,2] / min(egger.fit$sigma, 1)
  alpha.p.egger   = 2*(1-pt(abs(alpha.egger/alpha.se.egger),df.egger))
  alpha.lower.egger   = alpha.egger + (-1*qt(df=df.egger, 0.975)*alpha.se.egger)
  alpha.upper.egger   = alpha.egger + (1*qt(df=df.egger, 0.975)*alpha.se.egger)
  return(list(beta.egger=beta.egger,beta.se.egger=beta.se.egger,beta.lower.egger=beta.lower.egger,beta.upper.egger=beta.upper.egger,beta.t.egger=beta.egger/beta.se.egger,beta.p.egger=beta.p.egger,
              alpha.egger=alpha.egger,alpha.se.egger=alpha.se.egger,alpha.lower.egger=alpha.lower.egger,alpha.upper.egger=alpha.upper.egger,alpha.t.egger=alpha.egger/alpha.se.egger,alpha.p.egger=alpha.p.egger))
}
ModeEstimator <- function(x, y, sigmax, sigmay, phi=c(1,0.5,0.25), n_boot=1e4, alpha=0.05) {
  beta <- function(BetaIV.in, seBetaIV.in) {
    s <- 0.9*(min(sd(BetaIV.in), mad(BetaIV.in)))/length(BetaIV.in)^(1/5)
    weights <- seBetaIV.in^-2/sum(seBetaIV.in^-2)
    beta <- NULL
    for(cur_phi in phi) {
      h <- s*cur_phi
      densityIV <- density(BetaIV.in, weights=weights, bw=h)
      beta[length(beta)+1] <- densityIV$x[densityIV$y==max(densityIV$y)]
    }
    return(beta)
  }
  boot <- function(BetaIV.in, seBetaIV.in, beta_Mode.in) {
    beta.boot <- matrix(nrow=n_boot, ncol=length(beta_Mode.in))
    for(i in 1:n_boot) {
      BetaIV.boot      <- rnorm(length(BetaIV.in), mean=BetaIV.in, sd=seBetaIV.in[,1])
      BetaIV.boot_NOME <- rnorm(length(BetaIV.in), mean=BetaIV.in, sd=seBetaIV.in[,2])
      beta.boot[i,1:length(phi)]                     <- beta(BetaIV.in=BetaIV.boot, seBetaIV.in=rep(1, length(BetaIV)))
      beta.boot[i,(length(phi)+1):(2*length(phi))]   <- beta(BetaIV.in=BetaIV.boot, seBetaIV.in=seBetaIV.in[,1])
      beta.boot[i,(2*length(phi)+1):(3*length(phi))] <- beta(BetaIV.in=BetaIV.boot_NOME, seBetaIV.in=rep(1, length(BetaIV)))
      beta.boot[i,(3*length(phi)+1):(4*length(phi))] <- beta(BetaIV.in=BetaIV.boot_NOME, seBetaIV.in=seBetaIV.in[,2])
    }
    return(beta.boot)
  }
  BetaIV   <- y/x
  seBetaIV <- cbind(sqrt((sigmay^2)/(x^2) + ((y^2)*(sigmax^2))/(x^4)), sigmay/abs(x))
  beta_SimpleMode        <- beta(BetaIV.in=BetaIV, seBetaIV.in=rep(1, length(BetaIV)))
  beta_WeightedMode      <- beta(BetaIV.in=BetaIV, seBetaIV.in=seBetaIV[,1])
  beta_WeightedMode_NOME <- beta(BetaIV.in=BetaIV, seBetaIV.in=seBetaIV[,2])
  beta_Mode <- rep(c(beta_SimpleMode, beta_WeightedMode,
                     beta_SimpleMode, beta_WeightedMode_NOME))
  beta_Mode.boot <- boot(BetaIV.in=BetaIV, seBetaIV.in=seBetaIV, beta_Mode.in=beta_Mode)
  se_Mode <- apply(beta_Mode.boot, 2, mad)
  CIlow_Mode <- beta_Mode-qnorm(1-alpha/2)*se_Mode
  CIupp_Mode <- beta_Mode+qnorm(1-alpha/2)*se_Mode
  P_Mode <- pt(abs(beta_Mode/se_Mode), df=length(x)-1, lower.tail=F)*2
  Method <- rep(c('Simple', 'Weighted', 'Simple (NOME)', 'Weighted (NOME)'), each=length(phi))
  Results <- data.frame(Method, phi, beta_Mode, se_Mode, CIlow_Mode, CIupp_Mode, P_Mode)
  colnames(Results) <- c('Method', 'phi', 'Estimate', 'SE', 'CI_low', 'CI_upp', 'P')
  return(Results)
}
MR_output <- function(ivw,egger,median, mode) {
  output = data.frame(matrix(NA, nrow=5, ncol=7))
  names(output) = c("test", "parameter", "estimate", "se", "lower_CI", "upper_CI","p_value")
  output[1:5,1] = c("IVW","MR-Egger","MR-Egger","Weighted_median","Weighted_mode")
  output[1:5,2] = c("beta","beta","alpha","beta","beta")
  output[1,3:7] = c(IVW$beta.ivw,IVW$beta.se.ivw,IVW$beta.lower.ivw,IVW$beta.upper.ivw,IVW$beta.p.ivw)
  output[2,3:7] = c(Egger$beta.egger,Egger$beta.se.egger,Egger$beta.lower.egger,Egger$beta.upper.egger,Egger$beta.p.egger)
  output[3,3:7] = c(Egger$alpha.egger,Egger$alpha.se.egger,Egger$alpha.lower.egger,Egger$alpha.upper.egger,Egger$alpha.p.egger)
  output[4,3:7] = c(Median$beta.median,MedianBoot$beta.se.median,MedianBoot$beta.lower.median,MedianBoot$beta.upper.median,MedianBoot$beta.p.median)
  output[5,3:7] = c(Mode$Estimate[Mode$Method=="Weighted (NOME)" & Mode$phi==1.00],Mode$SE[Mode$Method=="Weighted (NOME)" & Mode$phi==1.00],Mode$CI_low[Mode$Method=="Weighted (NOME)" & Mode$phi==1.00],Mode$CI_upp[Mode$Method=="Weighted (NOME)" & Mode$phi==1.00],Mode$P[Mode$Method=="Weighted (NOME)" & Mode$phi==1.00])
  return(output)
}
```
<br>

**2. Use the functions to estimate the results.**

*As with Part 4, every "beta" below is a **log OR for CHD per 1-SD (~4.56 kg/m²) increase in BMI**, except for
the MR-Egger **alpha** (intercept), which is on the log OR scale too, but represents average pleiotropic bias
across the SNPs — i.e., how far the SNP-specific ratios' regression line sits away from zero when the SNP-BMI
effect is zero. An alpha significantly different from zero suggests directional horizontal pleiotropy.*
```{r}
IVW <- two.sample.iv.ivw(gp,gd,segp,segd)
Egger <- two.sample.iv.egger(gp,gd,segp,segd)
Median <- weighted.median(wald_ratio,weight)
MedianBoot <- weighted.median.boot(gp,gd,segp,segd,1000,0.05,weight)
Mode <- ModeEstimator(gp,gd,segp,segd)
sensitivity <- MR_output(IVW,Egger,Median,Mode)
sensitivity
write.table(sensitivity,"./twosample_sensitivity_BMI_CHD.txt",sep="\t",col.names=T,row.names=F,quote=F)
```

**IVW** - Inverse Variance Weighted - Combines the Wald ratios using an inverse variance weighted meta-analysis, where the weight of each ratio is the inverse of the variance of the association between the SNP and the outcome. It assumes that either every SNP is a valid instrument, or that any pleiotropy "balances out" across SNPs (balanced pleiotropy). It's the most efficient (precise) of these methods when its assumptions hold, so it's usually treated as the primary/main MR estimate, with the others used as sensitivity checks.
<br><br>
**MR-Egger** - combines the Wald ratio's together into a meta-regression to estimate the causal effect adjusted for any directional pleiotropy. This approach is less powered than the IVW, so its confidence interval is typically much wider — a non-significant MR-Egger result doesn't necessarily contradict a significant IVW result, it may simply reflect lower statistical power.
<br><br>
**MR-Egger (alpha)** - the intercept of the MR-egger meta-regression. Provides an indication of horizontal pleiotropy when it is not null (i.e., when its 95% CI excludes zero).
<br><br>
**Weighted_median** -  assigns a weight to each SNP derived from the inverse variance of each SNP's effect on the outcome. This robust method requires only 50% of the variants to be valid and not exhibit horizontal pleiotropy, unmeasured confounding, etc. — so, unlike IVW, it can still give a consistent estimate even if up to half of the SNPs are invalid instruments.
<br><br>
**Weighted_mode** - also assigns a weight to each SNP derived from the inverse variance of each SNP's effect on the outcome. It assumes that the largest sub-set of SNPs giving similar causal estimates (the "mode") are valid instruments, even if this isn't a strict majority — making it robust as long as no single subset of invalid instruments giving a similar (but wrong) answer is larger than the subset of valid instruments.
<br><br>
**Bringing it together:**

If IVW, weighted median, and weighted mode all point in a similar direction and magnitude, and the MR-Egger
intercept (alpha) is close to zero, this triangulation of evidence supports a genuine causal effect of BMI
on CHD of roughly the size given by the IVW/overall genetic effect estimate.

If the methods disagree substantially, or the MR-Egger intercept is significantly different from zero, this
points towards horizontal pleiotropy biasing at least some of the SNP-specific estimates, and the causal
conclusion should be treated with more caution — the weighted median/mode estimates (which are more robust
to a minority of invalid instruments) may then be more trustworthy than the IVW estimate.