JAMMR: JAMMR - Bayesian variable selection for Mendelian...

In pjnewcombe/R2BGLiMS: R interface to the java package BGLiMS (Bayesian Generalised Linear Model Selection).

Description

This function implements the JAM-MR algorithm for Mendelian randomization. Given a set of genetic variants and corresponding summary statistics for associations with a risk factor and a disease outcome of interest, JAM-MR performs Bayesian stochastic search through a reversible-jump MCMC algorithm to identify the most suitable variants for inclusion in a subsequent Mendelian randomization analysis. The algorithm prioritizes genetic variants with strong associations with the risk factor while downweighting variants with heterogeneous ratio estimates through the use of the general Bayesian framework and a heterogeneity- penalizing loss function. The stochastic search returns a list of models with associated posterior probabilties. For each of them, JAM-MR obtains a model-specific causal effect estimate by fitting a normal or (if mretn = TRUE) a truncated normal multiplictive random-effects model. Finally, the algorithm averages across the model-specific estimates to return an aggregate causal effect estimate and its standard error.

Usage

JAMMR(
  bx,
  sx,
  by,
  sy,
  N1,
  eafs = NULL,
  G.matrix = NULL,
  G.cor = NULL,
  trait.var = NULL,
  iter = 1e+06,
  w = NULL,
  n.grid = 26,
  grid.limits = NULL,
  initial.model = NULL,
  n.models = 10,
  mretn = TRUE,
  jam.model.priors = c(1, length(bx)),
  loss.function = "steve",
  jam.seed = NULL
)

Arguments

bx	Univariate effect estimates between each genetic variant and the risk factor.
sx	Standard errors for the varant-risk factor genetic effects.
by	Univariate effect estimates between each genetic variant and the outcome.
sy	Standard errors for the variant-outcome genetic effects.
N1	Sample size of the study from which the bx effects were obtained.
eafs	A vector of effect allele frequencies of length same as bx. Must be provided if genetic variants are assumed to be independent.
G.matrix	A reference matrix for the genetic variants used in the analysis. If present, will override eafs and force JAM-MR to model genetic variants as correlated. If absent, JAM-MR will assume that genetic variants are independent. One of eafs ad G.matrix must be specified. Note: the algorithm takes longer to run with coorrelated instruments, therefore when genetic variants are truly independent it is better to specify eafs and not G.matrix.
G.cor	The genetic correlation matrix. If an analysis with correlated SNPs is implemented but the matrix of raw genetic data (G.matrix) is unavailable or is too big to store in memory, JAM-MR can be run by providing both eafs and G.cor instead.
trait.var	An estimate of the variance in risk factor measurements. Can be obtained from the G-X GWAS and will be equal to 1 if the GWAS estimates are reported based on standardized data. If not provided, it is internally estimated by JAM-MR.
iter	The number of reversible-jump MCMC iterations to perform.
w	Tuning parameter(s) for pleiotropy penalization. If scalar, a single JAM-MR implementation will be run. If vector, one implementation will be run for each value and the minimum-causal-standard-error run will be selected. Values should typically be multiples of N1. Larger values encourage stronger penalization and smaller models.
n.grid	If w is not specified, it is estimated by a grid search; n.grid is the number of grid-points to visit during the grid search. Defaults to 26. Redundant if w is specified.
grid.limits	Lower and upper limits for the w values to consider during the grid search. The algorithm will create a grid of size n.grid by logarithmic interpolation between the limits, plus the value w = 0. The default is a grid search between w = 0.01 N1 and w = 100 N1. Redundant if w is specified.
initial.model	The model to be used in the first iteration of the stochastic search. By default the algorithm starts from a model with all genetic variants included, when running on independent variants, and a model with only the smallest p-value genetic variant included, when running on correlated variants.
n.models	The (maximum) number of highest posterior probability models to be reported by JAM-MR. If set to zero, no models will be reported.
mretn	Logical. If TRUE, the algorithm will fit a multiplicative random-effects truncated normal model in order to compute model-specific standard errors when averaging across models. If FALSE, random-effects IVW estimates will be used instead. Currently this is only implemented for independent genetic variants; with correlated variants the algorithm will always use IVW as model-specific estimates.
jam.model.priors	Prior parameters for JAM's beta-binomial model space prior.
loss.function	The loss function to use. Weighted if "steve", unweighted if "variance".
jam.seed	The seed to use for initializing the RJMCMC algorithm, if any. If a grid search for w is implemented, the implementations for different grid points will be seeded using the values jam.seed, jam.seed + 1, jam.seed + 2, ...

Value

A list with the following arguments:

• causal - Estimated causal effect.

• se - Causal standard error.

• model.matrix - If n.models > 0, a matrix containing the highest posterior probability models visited by JAM-MR during the variable selection procedure. See the function TopModels for more details.

• snp.probs - Posterior inclusion probabilities for each genetic variant.

• w - The value of the tuning parameter selected.

In addition, if a grid search for w is implemented:

• all.causal - Causal effect estimates for each w value.

• all.se - Causal standard errors for each w value.

• all.probs - Posterior inclusion probabilities per genetic variant for each w value.

• all.w - All w values used as part of the grid search.

Author(s)

Apostolos Gkatzionis and Paul Newcombe

See Also

See JAMMR_ManhattanPlot for a visual summary of covariate selection probabilities and JAMMR_SensitivityPlot for a visual comparison of results obtained using different w values. Also see JAM for an implementation of the fine-mapping algorithm on which JAM-MR relies.

Examples

### Using an artificial dataset.
set.seed(21)
P <- 50
bx.sim <- rnorm(P, 0.15, 0.05)
sx.sim <- abs(rnorm(P, 0.02, 0.008))
which.pl <- rbinom(P, 1, 0.3)   ## Which SNPs are pleiotropic.
by.sim <- bx.sim * 0.5 + rnorm(P, 0, 0.05)   ## True causal effect is 0.5.
by.sim[which.pl == 1] <- by.sim[which.pl == 1] + rep(0.25, sum(which.pl))   ## Pleiotropic effects.
sy.sim <- abs(rnorm(P, 0.04, 0.01))

### Visualize the data.
plot(bx.sim, by.sim, type = "p", pch = 19, xlim = c(0.00, 0.30), ylim = c(-0.10, 0.45), main = "Bx-By Plot")
for (i in 1:50) lines(c(bx.sim[i], bx.sim[i]), c(by.sim[i] - 1.96 * sy.sim[i], by.sim[i] + 1.96 * sy.sim[i]))
for (i in 1:50) lines(c(bx.sim[i] - 1.96 * sx.sim[i], bx.sim[i] + 1.96 * sx.sim[i]), c(by.sim[i], by.sim[i]))
abline(a = 0, b = median(by.sim / bx.sim), col = "brown", lty = 2)

### Run JAM-MR with a single value of w and independent data.
jammr.results1 <- JAMMR(bx.sim, sx.sim, by.sim, sy.sim, N1 = 10000, eafs = rep(0.3, 50),
                         iter = 1e6, w = 10000, jam.seed = 22)
jammr.results1$causal
jammr.results1$se

jammr.results1$snp.probs   ## Probabilities of inclusion per SNP.
jammr.results1$snp.probs[which.pl == 0]   ## Probabilities for valid SNPs.
jammr.results1$snp.probs[which.pl == 1]   ## Probabilities for pleiotropic SNPs.

### Run with multiple values of w.
jammr.results1.1 <- JAMMR(bx.sim, sx.sim, by.sim, sy.sim, N1 = 10000, eafs = rep(0.3, 50),
                           iter = 1e6, n.grid = 6, jam.seed = 22)
jammr.results1.1$causal
jammr.results1.1$se
jammr.results1.1$all.causal
jammr.results1.1$all.se
jammr.results1.1$snp.probs



### Using a real dataset from the MendelianRandomization package.
library(MendelianRandomization)

### Run existing MR methods.
mr_ivw( mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse) )
mr_egger( mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse) )
mr_median( mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse) )

### Run JAM-MR.
jammr.results <- JAMMR(ldlc, ldlcse, chdlodds, chdloddsse, N1 = 17723, eafs = lipid_eaf, iter = 1e6,
                        n.grid = 6, grid.limits = c(100, 20000), n.models = 0, jam.seed = 4189)
jammr.results$causal
jammr.results$se
jammr.results$snp.probs

pjnewcombe/R2BGLiMS documentation built on Feb. 10, 2020, 8:52 p.m.