3 The data: a world with a known answer

To judge an estimator honestly you need a case where you already know the right answer. Real studies never give you that. Simulation does. That makes it the natural setting for learning a method, not a compromise forced on us by data access.

The generator lives in R/simulate.R. It builds each dataset from an explicit causal model, so the true effect is whatever we say it is.

3.1 The causal model, in code

G_j  ~ Binomial(2, maf_j)                       # allele dosages: the instruments
U    ~ N(0, 1)                                   # unmeasured confounder
X    = sum_j bx_j * G_j + conf_x * U + noise     # exposure
Y    = THETA * X + sum_j alpha_j * G_j + conf_y * U + noise   # outcome

Three pieces deserve a comment. THETA is the causal effect we are hunting, fixed at 0.35. conf_x and conf_y are how strongly the confounder $U$ leaks into the exposure and the outcome; set them to zero and the naïve regression from Chapter 1 would be correct. alpha_j is the direct effect of variant $j$ on the outcome: the pleiotropy that breaks the exclusion restriction. When every alpha_j is zero the instruments are valid, which is the world we work in until Chapter 6.

sim <- simulate_mr(pleiotropy = "none")   # defaults: 10,000 people, 60 variants
str(sim$params[c("theta", "n", "n_snps", "conf_x", "conf_y")])

List of 5
 $ theta : num 0.35
 $ n     : num 10000
 $ n_snps: num 60
 $ conf_x: num 0.6
 $ conf_y: num 0.6

3.2 What the variables look like

The instruments are allele dosages: each variant is coded 0, 1, or 2 copies of an allele, with frequencies drawn to resemble common variants. The exposure and outcome are continuous.

g1 <- ggplot(data.frame(g = factor(sim$G[, 1])), aes(g)) +
  geom_bar(fill = mrcol("blue"), width = 0.7, alpha = 0.9) +
  labs(title = "One instrument (rs0001)", x = "Allele dosage", y = "Individuals") +
  theme_mr()
g2 <- ggplot(data.frame(x = sim$x), aes(x)) +
  geom_histogram(bins = 40, fill = mrcol("slate"), alpha = 0.85) +
  labs(title = "The exposure X", x = "Exposure", y = "Individuals") +
  theme_mr()
g1 + g2

Figure 3.1: Left: dosages for one variant take three values, as genotypes do. Right: the exposure is continuous and, thanks to many small genetic contributions plus noise, roughly bell-shaped.

3.3 Are the instruments strong enough?

Relevance (the one assumption we can check) is usually reported as an F-statistic from regressing the exposure on the instruments. The rule of thumb is that $F > 10$ keeps weak-instrument bias small. We can compute it directly.

fit  <- lm(sim$x ~ sim$G)
fs   <- summary(fit)$fstatistic
Fval <- unname(fs["value"])
sprintf("Instrument F-statistic = %.1f  (rule of thumb: comfortable above 10)", Fval)

[1] "Instrument F-statistic = 295.1  (rule of thumb: comfortable above 10)"

r2 <- summary(fit)$r.squared
sprintf("The instruments jointly explain %.1f%% of the variance in the exposure.", 100 * r2)

[1] "The instruments jointly explain 64.0% of the variance in the exposure."

An F well past 10 means these variants move the exposure decisively, so the ratios we build in later chapters rest on solid denominators. Weak instruments are one of the classic ways MR goes wrong, and it is reassuring to rule the problem out before going further.

3.4 The two shapes the data comes in

MR is done in two settings, and the rest of the book follows that split.

One-sample. Individual-level records for exposure, outcome, and genotype, all on the same people. This is Chapter 4, where we run two-stage least squares.
Two-sample. You never see individuals, only summary statistics: for each variant, its association with the exposure (from one study) and with the outcome (from another). This is how most published MR is done, because the summary numbers are freely shared while the raw data are not. Chapters 5 and 6 live here.

We can manufacture the summary-statistic form from our individual data by regressing each trait on each variant, one variant at a time. That is what summarise_snps() does.

ss <- summarise_snps(sim$G, sim$x, sim$y)
head(ss, 5)

	snp	bx	sx	by	sy
beta	rs0001	0.2299716	0.0289360	0.0876550	0.0212258
beta1	rs0002	0.3128183	0.0301328	0.1002602	0.0221489
beta2	rs0003	0.1790455	0.0387234	0.0312247	0.0283686
beta3	rs0004	0.5549845	0.0339791	0.1703367	0.0251398
beta4	rs0005	0.3612507	0.0301488	0.1019630	0.0221991

Each row is one variant: bx is its effect on the exposure, by its effect on the outcome, with standard errors. Everything from here on is built out of these two columns. On to the estimators, starting with the individual-level workhorse: two-stage least squares.

# The data: a world with a known answer ```{r} #| label: setup #| include: false source("R/theme_mr.R") source("R/simulate.R") source("R/mr_estimators.R") library(ggplot2) library(patchwork) ``` To judge an estimator honestly you need a case where you already know the right answer. Real studies never give you that. Simulation does. That makes it the natural setting for learning a method, not a compromise forced on us by data access. The generator lives in [`R/simulate.R`](https://github.com/Lexiyao/mendelian-randomisation-from-scratch/blob/main/R/simulate.R). It builds each dataset from an explicit causal model, so the true effect is whatever we say it is. ## The causal model, in code ```{r} #| label: model #| eval: false G_j ~ Binomial(2, maf_j) # allele dosages: the instruments U ~ N(0, 1) # unmeasured confounder X = sum_j bx_j * G_j + conf_x * U + noise # exposure Y = THETA * X + sum_j alpha_j * G_j + conf_y * U + noise # outcome ``` Three pieces deserve a comment. `THETA` is the causal effect we are hunting, fixed at `r TRUE_THETA`. `conf_x` and `conf_y` are how strongly the confounder $U$ leaks into the exposure and the outcome; set them to zero and the naïve regression from Chapter 1 would be correct. `alpha_j` is the *direct* effect of variant $j$ on the outcome: the pleiotropy that breaks the exclusion restriction. When every `alpha_j` is zero the instruments are valid, which is the world we work in until Chapter 6. ```{r} #| label: simulate sim <- simulate_mr(pleiotropy = "none") # defaults: 10,000 people, 60 variants str(sim$params[c("theta", "n", "n_snps", "conf_x", "conf_y")]) ``` ## What the variables look like The instruments are allele dosages: each variant is coded 0, 1, or 2 copies of an allele, with frequencies drawn to resemble common variants. The exposure and outcome are continuous. ```{r} #| label: fig-marginals #| fig-cap: "Left: dosages for one variant take three values, as genotypes do. Right: the exposure is continuous and, thanks to many small genetic contributions plus noise, roughly bell-shaped." g1 <- ggplot(data.frame(g = factor(sim$G[, 1])), aes(g)) + geom_bar(fill = mrcol("blue"), width = 0.7, alpha = 0.9) + labs(title = "One instrument (rs0001)", x = "Allele dosage", y = "Individuals") + theme_mr() g2 <- ggplot(data.frame(x = sim$x), aes(x)) + geom_histogram(bins = 40, fill = mrcol("slate"), alpha = 0.85) + labs(title = "The exposure X", x = "Exposure", y = "Individuals") + theme_mr() g1 + g2 ``` ## Are the instruments strong enough? Relevance (the one assumption we *can* check) is usually reported as an **F-statistic** from regressing the exposure on the instruments. The rule of thumb is that $F > 10$ keeps weak-instrument bias small. We can compute it directly. ```{r} #| label: fstat fit <- lm(sim$x ~ sim$G) fs <- summary(fit)$fstatistic Fval <- unname(fs["value"]) sprintf("Instrument F-statistic = %.1f (rule of thumb: comfortable above 10)", Fval) ``` ```{r} #| label: r2 r2 <- summary(fit)$r.squared sprintf("The instruments jointly explain %.1f%% of the variance in the exposure.", 100 * r2) ``` An F well past 10 means these variants move the exposure decisively, so the ratios we build in later chapters rest on solid denominators. Weak instruments are one of the classic ways MR goes wrong, and it is reassuring to rule the problem out before going further. ## The two shapes the data comes in MR is done in two settings, and the rest of the book follows that split. - **One-sample.** Individual-level records for exposure, outcome, and genotype, all on the same people. This is Chapter 4, where we run two-stage least squares. - **Two-sample.** You never see individuals, only *summary statistics*: for each variant, its association with the exposure (from one study) and with the outcome (from another). This is how most published MR is done, because the summary numbers are freely shared while the raw data are not. Chapters 5 and 6 live here. We can manufacture the summary-statistic form from our individual data by regressing each trait on each variant, one variant at a time. That is what `summarise_snps()` does. ```{r} #| label: summary-stats ss <- summarise_snps(sim$G, sim$x, sim$y) head(ss, 5) ``` Each row is one variant: `bx` is its effect on the exposure, `by` its effect on the outcome, with standard errors. Everything from here on is built out of these two columns. On to the estimators, starting with the individual-level workhorse: [two-stage least squares](04-one-sample.qmd).