IB-PRESSO: instrument borrowing for outlier-robust Mendelian randomization
Source:R/IBPRESSO.R
IBPRESSO.RdExtends MR-PRESSO (Verbanck et al., 2018) by scoring each instrument on a primary outcome AND an auxiliary outcome simultaneously. Every other element of the procedure – leave-one-out IVW fitting through the origin, the global heterogeneity bootstrap, the Bonferroni correction across instruments, and the outlier-corrected IVW estimate – is MR-PRESSO's, unchanged, so that any difference in performance is attributable to the auxiliary trait and to nothing else.
Usage
IBPRESSO(
BetaOutcome,
BetaExposure,
BetaAux,
SdOutcome,
SdExposure,
SdAux,
data,
OUTLIERtest = TRUE,
DISTORTIONtest = FALSE,
SignifThreshold = 0.05,
NbDistribution = 10000,
seed = NULL,
n_cores = 1,
CorOutcomeAux = 0,
OutlierTestMethod = c("chisq", "bootstrap")
)Arguments
- BetaOutcome
Character, column name for primary outcome SNP effects.
- BetaExposure
Character, column name for exposure SNP effects (one only).
- BetaAux
Character, column name for auxiliary trait SNP effects.
- SdOutcome
Character, column name for standard errors of outcome effects.
- SdExposure
Character, column name for standard errors of exposure effects.
- SdAux
Character, column name for standard errors of auxiliary effects.
- data
Data frame containing all the above columns, one row per instrument.
- OUTLIERtest
Logical; perform the per-SNP outlier test and report the corrected estimate. Default
TRUE.- DISTORTIONtest
Logical; test whether the removed instruments distort the estimate more than a random set of the same size. Default
FALSE.- SignifThreshold
Numeric significance threshold, Bonferroni-corrected across instruments for the per-SNP test. Default 0.05.
- NbDistribution
Integer, number of bootstrap null datasets for the GLOBAL test. Default 10000.
- seed
Integer random seed, or
NULL.- n_cores
Retained for backward compatibility and ignored.
- CorOutcomeAux
Cross-outcome sampling correlation \(\rho_j\): the cross-trait LDSC intercept on the CORRELATION scale, \(I_{12}/\sqrt{I_{11} I_{22}}\). Scalar, or a vector of length
nrow(data)or of the number of complete instruments. Must lie strictly inside \((-1, 1)\); a raw LDSC intercept generally does not, and is not a correlation. Default 0, which reproduces independent outcome GWAS.- OutlierTestMethod
"chisq"(default) for the analytic per-SNP null, or"bootstrap"for the previous Monte Carlo route.
Value
An object of class mrpresso_ib, a list with components:
- raw_beta, raw_se
IVW estimate and standard error using all instruments.
- corrected_beta, corrected_se
Outlier-corrected estimate, or
NAwhen no outlier was removed. Callers should fall back to the raw estimate.- p_value
Global heterogeneity test p-value (bootstrap).
- p_distort
Distortion test p-value, or
NA.- outlier_idx
Row indices of the removed instruments.
- outlyingness
Observed per-SNP statistic \(T_j\).
- outlier_pvals
Bonferroni-corrected per-SNP p-values that were used for flagging; a copy of whichever of the next two
outlier_test_methodnames.- outlier_test_method
"chisq"or"bootstrap".- outlier_pvals_chisq
Analytic per-SNP p-values. Always computed.
- outlier_pvals_boot
Bootstrap per-SNP p-values, or all
NAwhen the analytic route was used.- outcome_correlation
The \(\rho_j\) actually used, after alignment.
- covariance_model
Free-text description of the covariance assumption.
- n_instruments, n_outliers
Counts.
Per-SNP statistic
Let \(r_{1j}\) and \(r_{2j}\) be the leave-one-out residuals of the primary and auxiliary outcomes at instrument \(j\), on the raw scale of the reported effect estimates. Conditional on the leave-one-out slopes \(b_{1,-j}\) and \(b_{2,-j}\) their sampling covariance is known: $$V_j = \begin{pmatrix} s_{Y_1j}^2 + b_{1,-j}^2 s_{Xj}^2 & \rho_j s_{Y_1j} s_{Y_2j} + b_{1,-j} b_{2,-j} s_{Xj}^2 \\ \rho_j s_{Y_1j} s_{Y_2j} + b_{1,-j} b_{2,-j} s_{Xj}^2 & s_{Y_2j}^2 + b_{2,-j}^2 s_{Xj}^2 \end{pmatrix},$$ and the statistic is the generalized squared residual $$T_j = (r_{1j}, r_{2j}) V_j^{-1} (r_{1j}, r_{2j})^\top = \frac{v_{2j} r_{1j}^2 - 2 c_{12,j} r_{1j} r_{2j} + v_{1j} r_{2j}^2} {v_{1j} v_{2j} - c_{12,j}^2}.$$ MR-PRESSO's statistic is the primary term alone. Two dependence sources enter \(V_j\), each exactly once: outcome-GWAS sample overlap through \(\rho_j\), and shared exposure noise through \(b_{1,-j} b_{2,-j} s_{Xj}^2\), which is present even when \(\rho_j = 0\).
Why the covariance is not estimated
\(V_j\) is built from the reported GWAS standard errors and a cross-trait LDSC intercept, never from the observed residuals. Estimating a \(2 \times 2\) residual covariance (robustly or otherwise) fails in three ways: the invalid instruments inflate the estimate along precisely the direction that identifies them; the resulting Mahalanobis distance is bounded by \(1/\pi\) where \(\pi\) is the fraction of invalid instruments, so at \(\pi = 0.5\) no instrument can exceed 2 against a cutoff of 5.99; and the minimum covariance determinant has a breakdown point of 50 percent.
Reference distributions
\(T_j\) is pivotal – a quadratic form in a bivariate normal with known
covariance – so conditional on the leave-one-out slopes its null is exactly
\(\chi^2_2\) and the per-SNP p-value is analytic by default. This removes
the bootstrap resolution floor of \(K/\code{NbDistribution}\), which at
\(K = 1270\) and 10000 draws is 0.127 and makes flagging impossible at any
sane threshold. The GLOBAL test keeps the bootstrap, because the
leave-one-out slopes are refitted on every draw and \(\sum_j T_j\) has no
closed form. Setting OutlierTestMethod = "bootstrap" restores the previous
per-SNP behaviour; both sets of p-values are always returned.
Neither route propagates the uncertainty in the leave-one-out slopes, so on self-fitted null data \(T_j\) is inflated by roughly \(1 + 2/K\) (measured mean 2.03 at \(K = 100\), 2.02 at \(K = 200\), against 2).
Weights
Two weights, with a deliberate division of labour. \(W_j = 1/(s_{Y_1j} s_{Y_2j})\) fits the leave-one-out slopes, a single weight applied to both outcomes so that neither trait is privileged; it does not enter the statistic, since \(V_j\) carries all the scaling. The causal estimate is an IVW fit weighted by \(1/s_{Y_1j}^2\), which is MR-PRESSO's own weight; the final estimator is therefore identical to MR-PRESSO's and the two methods differ only in which instruments they retain.
References
Verbanck M, Chen CY, Neale B, Do R (2018). Detection of widespread horizontal pleiotropy in causal relationships inferred from Mendelian randomization between complex traits and diseases. Nature Genetics 50, 693-698.
Examples
set.seed(1)
K <- 40
bx <- runif(K, 0.02, 0.05)
g <- c(rnorm(12, 0.012, 0.004), rep(0, K - 12))
d <- data.frame(
BetaExposure = rnorm(K, bx, 1/sqrt(1e5)), SdExposure = 1/sqrt(1e5),
BetaOutcome = rnorm(K, 0.1 * bx + g, 1/sqrt(5e4)), SdOutcome = 1/sqrt(5e4),
BetaAux = rnorm(K, 0.3 * bx + g, 1/sqrt(1e5)), SdAux = 1/sqrt(1e5)
)
IBPRESSO("BetaOutcome", "BetaExposure", "BetaAux",
"SdOutcome", "SdExposure", "SdAux",
data = d, NbDistribution = 2000, seed = 1, CorOutcomeAux = 0.3)
#> $raw_beta
#> [1] 0.1971852
#>
#> $raw_se
#> [1] 0.02764942
#>
#> $corrected_beta
#> [1] 0.1763317
#>
#> $corrected_se
#> [1] 0.02019307
#>
#> $p_value
#> [1] 0
#>
#> $p_distort
#> [1] NA
#>
#> $outlier_idx
#> [1] 1 5 10 11 29
#>
#> $outlyingness
#> [1] 33.25822072 12.18123801 11.17121753 0.05786788 15.14654017 7.18946394
#> [7] 7.94597913 0.69915812 3.82600252 26.53512227 29.18893312 7.36666237
#> [13] 0.36642495 0.27543522 1.18449415 1.76120563 1.93712145 2.71332315
#> [19] 1.73129527 1.95178069 0.31050107 3.03632990 1.68699309 7.07799187
#> [25] 3.39113823 1.39326998 3.46675986 1.36459230 13.91376996 1.07178933
#> [31] 7.18038110 1.44174514 2.89995586 4.39098916 3.73374830 1.60738905
#> [37] 4.50461617 11.01249254 2.71819784 1.15068627
#>
#> $outlier_pvals
#> [1] 2.399546e-06 9.056028e-02 1.500586e-01 1.000000e+00 2.056035e-02
#> [6] 1.000000e+00 7.526838e-01 1.000000e+00 1.000000e+00 6.918810e-05
#> [11] 1.835539e-05 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [16] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [21] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [26] 1.000000e+00 1.000000e+00 1.000000e+00 3.808231e-02 1.000000e+00
#> [31] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [36] 1.000000e+00 1.000000e+00 1.624530e-01 1.000000e+00 1.000000e+00
#>
#> $n_instruments
#> [1] 40
#>
#> $n_outliers
#> [1] 5
#>
#> $outcome_correlation
#> [1] 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3
#> [20] 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3
#> [39] 0.3 0.3
#>
#> $outlier_test_method
#> [1] "chisq"
#>
#> $outlier_pvals_chisq
#> [1] 2.399546e-06 9.056028e-02 1.500586e-01 1.000000e+00 2.056035e-02
#> [6] 1.000000e+00 7.526838e-01 1.000000e+00 1.000000e+00 6.918810e-05
#> [11] 1.835539e-05 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [16] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [21] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [26] 1.000000e+00 1.000000e+00 1.000000e+00 3.808231e-02 1.000000e+00
#> [31] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
#> [36] 1.000000e+00 1.000000e+00 1.624530e-01 1.000000e+00 1.000000e+00
#>
#> $outlier_pvals_boot
#> [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#> [26] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#>
#> $covariance_model
#> [1] "analytic residual covariance: outcome LDSC intercept + shared exposure error"
#>
#> attr(,"class")
#> [1] "mrpresso_ib" "list"