Model Assumptions

For the first research question, we introduced a mixture model to account for effect heterogeneity of the instruments. This model has several components:

1. Measurement error: SNP-exposure effects and SNP-outcome effects are observed with (known) measurement error with unknown means.

2. Independence: All the estimated SNP associations are mutually independent. The independence between the SNP-exposure and SNP-outcome effect for a given SNP is guaranteed if they are computed using non-overlapping samples. Independence across SNPs is a reasonable assumption if the selected SNPs are in linkage equilibrium.

3. Effect heterogeneity: We assume each instrument may estimate a different “causal” effect and the instrument-specific effects follow a Gaussian mixture distribution.

Assumptions 1 & 2 are standard in MR analysis, see for example Zhao et al. (2018). Assumption 3 differs from the standard MR assumption and the mixture distribution captures the idea that instruments on the same genetic pathway may correspond to similar causal effect estimates. Small and balanced horizontal pleiotropy is captured by the within-cluster variation.

Statistical methods

1. Visual inspection of effect heterogeneity: We first use a modal plot developed in Wang et al. (forthcoming, implemented in the modal.plot function from the mr.raps package) to visualize the data. This plot shows a robustified log-likelihood function over the causal effect and will exhibit multiple modes when there are clustered effect heterogeneity. Visual inspection of the scatterplot between estimated SNP-exposure and SNP-outcome effects can also suggest heterogeneity among variant-specific estimates.

2. Bayesian inference for the mixture model: We used a Monte-Carlo EM algorithm to estimate the prior parameters in the mixture model and then obtained (empirical Bayes) posterior distribution of the instrument-specific effects. See the Technical Appendix below for more detail.

3. Examine heatmap of lipoprotein subfractions: We ranked the variants based on the posterior mean of their instrument-specific effects and plot a heatmap of the associations of the SNPs with lipoprotein subfractions. This heatmap may provide some insights if groups of SNPs with similar instrument-specific effects correspond to distinct biological mechanisms.

Conclusions

• The 2 SNPs with the highest (posterior) probability of belonging to the 2nd cluster (rs1532085, rs588136) are both in chromosome 15 and are close to the LIPC gene.

• rs1532085, rs588136, and rs174546 (cluster 2) are positively associated with all of the large and extra-large HDL subfraction traits. (p-value \(\leq 10^{-11}\)).

• rs1532085, rs588136, and rs174546 are negatively associated with concentration of small HDL particles (p-value \(\leq 10^{-3}\)). rs1532085 and rs588136 are negatively associated with total lipids in small HDL (p-value \(\leq 10^{-14}\)).

• rs1532085, rs588136, and rs174546 are positively associated with LDL diameter (p-value \(\leq 10^{-7}\)).

• rs14320985 and rs588136 are positively associated with all of the very small VLDL subfraction traits (p-value \(\leq 10^{-14}\)). rs174546 is positively associated with total lipids and phospholipids in very small VLDL (p-value \(\leq 10^{-3}\)).

• rs1532085, rs588136, and rs174546 are positively associated with all of the IDL subtraction traits (p-value \(\leq 10^{-4}\)).

Model fitting Procedure

1. Obtain Empirical Bayes estimates \(\hat{\theta}\) of \(\theta\) using a Monte-Carlo EM algorithm (details provided below).

2. Sample from \(P(\mu_{X_i}, \beta_i, Z_i | X_i, Y_i, \hat{\theta})\), the posterior distribution of the latent variables given \(\hat{\theta}\) (details provided below).

Posterior Sampling

We will use importance sampling (Liu 2004) to sample from \(P(\mu_{X_i}, \beta_i, Z_i | X_i, Y_i, \theta)\). Let \(\phi(\cdot ; \mu, \sigma^2)\) denote the Gaussian density with mean \(\mu\) and variance \(\sigma^2\).

First, we can re-write the latent variable posterior as

\[ \begin{align*} P(\beta_i, Z_i = k, \mu_{X_i} | X_i, Y_i, \theta) & = P(\mu_{X_i}| X_i, Y_i, \theta) P(\beta_i | \mu_{X_i}, Y_i, \theta) P(Z_i = k | \beta_i, \theta) \\ & = P(\mu_{X_i} | X_i, Y_i, \theta) \Big[ \sum_{k=1}^K \pi_k \phi(\beta_i; \tilde{\mu}_{ik}, \tilde{\sigma}_{ik}^2) \Big] \Big[ \frac{\pi_k \phi(\beta_i; \mu_k, \sigma_k^2)}{\sum_{j=1}^K \pi_j \phi(\beta_i; \mu_j, \sigma_j^2)} \Big] \\ \end{align*} \] where

\[ \begin{align*} \tilde{\sigma}_{ik}^2 &= \Big(\frac{1}{\sigma_k^2} + \frac{\mu_{X_i}^2}{\sigma_{Y_i}^2}\Big)^{-1} \\ \tilde{\mu}_{ik} & = \tilde{\sigma}_{ik}^2 \Big( \frac{Y_i \mu_{X_i}}{\sigma_{Y_i}^2} + \frac{\mu_k}{\sigma_k^2} \Big) \end{align*} \]

Thus, we can sample from \(P(\mu_{X_i} | X_i, Y_i, \theta)\) using a (importance) proposal distribution, then sequentially sample from \(P(\beta_i | \mu_{X_i}, Y_i, \theta)\) and \(P(Z_i = k | \beta_i, \theta)\) directly, and weight the samples according to the importance weights.

\[ \begin{align*} P(\mu_{X_i} | X_i, Y_i, \theta) & \propto P(\mu_{X_i} | X_i, \theta) P(Y_i | \mu_{X_i}, \theta) \\ & = \phi(\mu_{X_i}; \tilde{m}_{X_i},\tilde{\lambda}_{X_i}^2) \Big[ \sum_{k=1}^K \pi_k \phi(Y_i; \mu_{X_i}\mu_k, \mu_{X_i}^2 \sigma_k^2 + \sigma_{Y_i}^2) \Big] =: \pi(\mu_{X_i}) \\ \end{align*} \]

where

\[ \begin{align*} \tilde{\lambda}_{X_i}^2 &= \Big(\frac{1}{\sigma_{X_i}^2} + \frac{1}{\lambda_x^2}\Big)^{-1} \\ \tilde{m}_{X_i} & = \tilde{\lambda}_{X_i}^2 \Big(\frac{X_i}{\sigma_{X_i}^2} + \frac{m_x}{\lambda_x^2}\Big) \end{align*} \]

We chose to use \(g(\mu_{X_i}) = P(\mu_{X_i} | X_i, \theta)\) as the proposal distribution with associated importance weights

\[ w_i = \frac{\pi(\mu_{X_i})}{g(\mu_{X_i})} = \sum_{k=1}^K \pi_k \phi(Y_i; \mu_{X_i}\mu_k, \mu_{X_i}^2 \sigma_k^2 + \sigma_{Y_i}^2) \] Since \(\phi(Y_i; \mu_{X_i}\mu_k, \mu_{X_i}^2 \sigma_k^2 + \sigma_{Y_i}^2) \leq \frac{1}{\sqrt{2\pi \sigma_{Y_i}^2}}\), \[ 0 \leq w_i \leq \frac{1}{\sqrt{2\pi \sigma_{Y_i}^2}} \] Therefore, the importance weights are bounded.

Monte-Carlo EM Algorithm

a Briefly, the EM algorithm (Dempster, Laird, and Rubin 1977) is an iterative algorithm to compute maximum likelihood estimates in latent variable models that consists of two steps:

• E Step: Compute the expectation of the complete-data log-likelihood with respect to the posterior distribution of the latent variables, given the parameter estimates from the previous iteration.

• M Step: Maximize the expectation computed in the E step with respect to the parameters.

In our model, the expectation in the E step does not have an analytical form so we resort to Monte-Carlo methods to approximate it (Levine and Casella 2001).

The complete-data log-likelihood is given by

\[ \begin{align*} l(\theta | \mathbf{X}, \mathbf{Y}, \{\beta_i\}, \{Z_i\}, \{\mu_{X_i}\}) & := \sum_{i=1}^p l(\theta | X_i, Y_i, \beta_i, Z_i, \mu_{X_i}) \\ & = \sum_{i=1}^p \Big\{ \text{log} \phi(\mu_{X_i}; m_x, \lambda_x^2) + \sum_{k=1}^K Z_{ik} \big[\text{log} \pi_k + \text{log} \phi(\beta_i; \mu_k, \sigma_k^2) \big] \Big\} \end{align*} \] where \(Z_{ik} = 1\) if \(Z_i = k\) and \(Z_{ik} = 0\) otherwise.

To approximate the E-step expectation, we approximate \(E(\mu_{X_i} | X_i, Y_i, \theta)\), \(E(\mu_{X_i}^2 | X_i, Y_i, \theta)\), \(E(Z_{ik} \beta_i | X_i, Y_i, \theta)\), \(E(Z_{ik} \beta_i^2 | X_i, Y_i, \theta)\), and \(E(Z_{ik} | X_i, Y_i, \theta)\) using the importance sampling procedure described in the previous section. After obtaining the importance sampling estimates, the M step is straightforward so we omit the details here.