Package {MR.RGM}

Posterior probability of a network motif

Description

NetworkMotif quantifies the uncertainty of a given network structure among the response variables: it returns the proportion of posterior draws in which every edge of the motif is present, that is, the posterior probability that all edges of the motif are present simultaneously. The input draws are the GammaPst output of RGM.

Usage

NetworkMotif(Gamma, GammaPst)

Arguments

Value

A number between 0 and 1: the proportion of posterior draws that contain all edges of Gamma.

References

Ni, Y., Ji, Y., & Müller, P. (2018). Reciprocal graphical models for integrative gene regulatory network analysis. Bayesian Analysis, 13(4), 1095-1110. doi:10.1214/17-BA1087.

Examples

# Simulate a cascade Y1 -> Y2 -> Y3 with two own instruments per response
set.seed(9154)
n = 2000; p = 3; k = 6
A = matrix(0, p, p); A[2, 1] = 0.3; A[3, 2] = -0.3     # A[i, j]: effect of Y_j on Y_i
D = matrix(0, p, k)
for (i in 1:p) D[i, (2 * i - 1):(2 * i)] = 1
X = matrix(rnorm(n * k), n, k)
Y = (X %*% t(D) + matrix(rnorm(n * p), n, p)) %*% t(solve(diag(p) - A))
X = scale(X); Y = scale(Y)                              # the model has no intercept

# Fit MR.RGM (a short chain to keep the example fast; the defaults are the recommended length)
fit = RGM(X = X, Y = Y, D = D, nIter = 2000, nBurnin = 1000, Thin = 2)

# Posterior probability of the cascade Y1 -> Y2 -> Y3 (both edges present)
Gamma = matrix(0, p, p); Gamma[2, 1] = 1; Gamma[3, 2] = 1
NetworkMotif(Gamma = Gamma, GammaPst = fit$GammaPst)

# Posterior probability of the feedback loop Y1 -> Y2 -> Y3 -> Y1
Gamma[1, 3] = 1
NetworkMotif(Gamma = Gamma, GammaPst = fit$GammaPst)

Fitting Multivariate Bidirectional Mendelian Randomization Networks Using Bayesian Directed Cyclic Graphical Models

Description

RGM fits the Bayesian directed cyclic graphical model

Y = A Y + B X + C U + E,\qquad E \sim N_p(0, \Sigma^{\ast}),

where Y are p response variables (traits), X are k instrumental variables (for example genetic variants), U are l optional measured covariates, A (p x p, zero diagonal) contains the direct causal effects among the responses (A[i, j] is the effect of Y_j on Y_i; directed cycles are allowed), B (p x k) the instrument effects, C (p x l) the covariate effects and \Sigma^{\ast} the error covariance, whose off-diagonal entries represent unmeasured confounding between responses.

Two versions of the method differ only in how the instruments are mapped to the responses (in the accompanying methods paper they are called IV-RGM and IV-RGM+, since the model applies to instrumental variables in general):

MR.RGM (known instrument map)

D is a binary p x k matrix with D[i, j] = 1 if instrument j is allowed to affect response i (for example a cis-eQTL for its own gene). Entries of B with D = 0 are fixed at zero (the exclusion restriction) and the allowed entries receive a slab prior N(0, eta) without any selection: every allowed instrument stays in the model. This is the default when D contains zeros (selectB = FALSE).

MR.RGM+ (instrument selection)

D is a matrix of ones: every instrument may affect every response, and each entry of B receives a spike-and-slab prior whose inclusion indicator is sampled, so that the instrument–response map is learned from the data (horizontal pleiotropy). This is the default when D has no zeros (selectB = TRUE).

The argument selectB overrides these defaults (see its description).

Inference is by MCMC on the sufficient statistics of the Gaussian model, so the function accepts individual-level data (X, Y, optionally U), the second-moment matrices (Syy, Syx, Sxx, optionally Syu, Sxu, Suu) computed from one sample, or the matrices Sxx, Beta and SigmaHat of marginal instrument–response regressions; responses alone (Y or Syy) fit the model without instruments. The marginal-regression input does not contain the covariances between the responses: they are reconstructed under independent errors (no unmeasured confounding), so this input requires SigmaStarModel = "diagonal", no covariates and, for every response, an instrument that affects only that response, and it is exact only for uncorrelated instruments. Use individual-level data or the second-moment matrices to model confounding. The model has no intercept and the stochastic-search prior on \Sigma^{\ast} assumes responses on roughly unit scale: center and scale the columns of Y, X and U (or compute the summary matrices from centered and scaled data) before fitting.

The output contains posterior means of all parameters, posterior inclusion probabilities of the causal edges (GammaEst), of the instrument effects (PhiEst, MR.RGM+ only) and of the confounding links (ZEst, SSSL only), the graph estimated by thresholding the inclusion probabilities at 0.5, the posterior draws of the causal edge indicators (GammaPst, used by NetworkMotif), the log-likelihood trace, acceptance rates and convergence/identification diagnostics. Warnings are issued for signs of weak identification or poor convergence (see diagnostics in the Value section and the Note).

MCMC algorithm. Every iteration updates, in turn, the instrument effects B, the covariate effects C, the error covariance \Sigma^{\ast} (with its confounding indicators) and the causal effects A, each with its prior hyperparameters. Given the other parameters, the likelihood is quadratic in B and its (spike-and-slab or slab) prior is normal given the inclusion indicators, so each row of B is drawn exactly from its multivariate normal full conditional (Gibbs; gibbsB = TRUE). The entries of A enter the likelihood also through \log|\det(I - A)| and are updated one at a time by random-walk Metropolis–Hastings with adaptive proposals: each entry has its own proposal variance, started at 2.38^2 times an approximation of its posterior variance and tuned during the burn-in towards an acceptance rate of 0.44 (adaptive Metropolis-within-Gibbs of Roberts and Rosenthal, 2009, with larger adaptation steps early in the burn-in); the variances are frozen after the burn-in, so the kept draws come from a fixed Markov kernel. With more than 100 allowed instruments for a response, its row of B is drawn entry by entry (each entry exactly from its scalar full conditional) instead of jointly. The same adaptive scheme updates B when it cannot be drawn by Gibbs (Threshold prior with selection, or gibbsB = FALSE). The chain starts at A = 0, at the least-squares effects of each response on its allowed instruments (known map; B = 0 with selection) and at the corresponding residual covariance.

Chain length and convergence. One call runs one chain. The defaults (50,000 iterations, 10,000 burn-in, thinning 10) are the settings validated in the simulation study of the accompanying paper for networks of up to about 20 responses; larger networks or weak instruments may need longer chains. For a reported analysis, run several chains (for example four) with different seeds (set.seed() before each call), check their agreement (for example the rank-normalized split-\hat R of LLPst and of the edge indicators in GammaPst, computed with the package posterior, and the posterior inclusion probabilities of the chains), and pool the chains only if they agree. The function also reports diagnostics and warnings of its own (see diagnostics in the value).

Usage

RGM(
  X = NULL,
  Y = NULL,
  U = NULL,
  Syy = NULL,
  Syx = NULL,
  Sxx = NULL,
  Syu = NULL,
  Sxu = NULL,
  Suu = NULL,
  Beta = NULL,
  SigmaHat = NULL,
  D,
  n,
  nIter = 50000,
  nBurnin = 10000,
  Thin = 10,
  prior = c("Spike and Slab", "Threshold"),
  SigmaStarModel = c("SSSL", "IW", "diagonal"),
  aRho = 3,
  bRho = 1,
  nu1 = 0.001,
  aPsi = 0.5,
  bPsi = 0.5,
  nu2 = 1e-04,
  aSigma = 0.01,
  bSigma = 0.01,
  PropVarA = NULL,
  PropVarB = NULL,
  gibbsB = TRUE,
  selectB = NULL,
  tauC = 2,
  omega1 = 5,
  omega2 = 0.005,
  varpi = 0.5,
  lambda = 1
)

Arguments

Value

A list with the following components (components that do not apply to the chosen prior or model are omitted):

Note

Weak identification. With D = 1 (MR.RGM+), the instrument effects are identified only through sparsity and an instrument-diversity condition, which require instruments that explain a non-negligible share of the variance of several responses. With weak, response-specific instruments (for example cis-eQTLs) the posterior can stay close to the prior; in that case RGM warns (error variances far above the second moments of the responses) and the instrument map D should encode which instrument affects which response (MR.RGM). Compare the posterior inclusion probabilities with their prior values (0.75 for causal edges under the default aRho = 3, bRho = 1; 0.5 for confounding links) and the causal inclusion probabilities computed from the two halves of GammaPst.

Identification of MR.RGM. With a known map, the causal effects are identified when, for every response, the instruments excluded from its equation satisfy the classical rank condition; a sufficient condition is that every response has at least one instrument that affects only that response. RGM warns when this sufficient condition fails and no selection is performed.

References

Roberts, G. O., & Rosenthal, J. S. (2009). Examples of adaptive MCMC. Journal of Computational and Graphical Statistics, 18(2), 349-367. doi:10.1198/jcgs.2009.06134.

Wang, H. (2015). Scaling it up: stochastic search structure learning in graphical models. Bayesian Analysis, 10(2), 351-377. doi:10.1214/14-BA916.

Ni, Y., Ji, Y., & Müller, P. (2018). Reciprocal graphical models for integrative gene regulatory network analysis. Bayesian Analysis, 13(4), 1095-1110. doi:10.1214/17-BA1087.

Sarkar, B. & Ni, Y. (2025). MR.RGM: an R package for fitting Bayesian multivariate bidirectional Mendelian randomization networks. Bioinformatics, 41(4), btaf130. doi:10.1093/bioinformatics/btaf130.

Examples

# ------------------------------------------------------------------------------------------
# Simulated data: Y1 -> Y2 -> Y3 and Y3 -> Y1 (a feedback loop), two own instruments per
# response, one covariate, and an unmeasured confounder shared by Y1 and Y3
set.seed(9154)
n = 2000; p = 3; k = 6; l = 1
A = matrix(0, p, p); A[2, 1] = 0.3; A[3, 2] = -0.3; A[1, 3] = 0.2   # A[i, j]: effect of Y_j on Y_i
D = matrix(0, p, k)                                                  # known instrument map
for (i in 1:p) D[i, (2 * i - 1):(2 * i)] = 1
X = matrix(rnorm(n * k), n, k); U = matrix(rnorm(n * l), n, l); W = rnorm(n)
E = cbind(W, 0, W) + matrix(rnorm(n * p), n, p)                      # confounding of Y1 and Y3
Y = (X %*% t(D) + U %*% matrix(0.5, l, p) + E) %*% t(solve(diag(p) - A))
colnames(Y) = paste0("Y", 1:p); colnames(X) = paste0("SNP", 1:k); colnames(U) = "age"
X = scale(X); Y = scale(Y); U = scale(U)                            # the model has no intercept

# Example 1: MR.RGM (known map D) from individual-level data, spike-and-slab prior and
# SSSL error covariance (the defaults); a short chain to keep the example fast, whereas the
# default chain length is the recommended one
fit = RGM(X = X, Y = Y, U = U, D = D, nIter = 2000, nBurnin = 1000, Thin = 2)
round(fit$GammaEst, 2)   # posterior inclusion probabilities of the causal edges
round(fit$AEst, 2)       # posterior means of the causal effects
round(fit$ZEst, 2)       # posterior inclusion probabilities of the confounding links
fit$diagnostics          # weak-instrument, stability and convergence diagnostics
plot(fit$Graph, main = "Estimated causal network")
plot(fit$LLPst, type = "l", xlab = "Kept draw", ylab = "Log-likelihood")

# Example 2: the same model from the second-moment matrices of the same sample
fit_s = RGM(Syy = crossprod(Y) / n, Syx = crossprod(Y, X) / n, Sxx = crossprod(X) / n,
            Syu = crossprod(Y, U) / n, Sxu = crossprod(X, U) / n, Suu = crossprod(U) / n,
            D = D, n = n, nIter = 2000, nBurnin = 1000, Thin = 2)
round(fit_s$GammaEst, 2)

# Example 3: MR.RGM+ (D = 1): every instrument may affect every response and the instrument
# effects are selected by a spike-and-slab prior (horizontal pleiotropy)
fit_plus = RGM(X = X, Y = Y, U = U, D = matrix(1, p, k), nIter = 2000, nBurnin = 1000, Thin = 2)
round(fit_plus$PhiEst, 2)   # posterior inclusion probabilities of the instrument effects
fit_plus$settings$selectB   # TRUE

# Example 4: threshold prior on A and a diagonal error covariance
fit_t = RGM(X = X, Y = Y, D = D, prior = "Threshold", SigmaStarModel = "diagonal",
            nIter = 2000, nBurnin = 1000, Thin = 2)
fit_t$zAEst

# Example 5: marginal-regression input (Sxx, Beta, SigmaHat), with uncorrelated instruments,
# no covariates and independent errors (SigmaStarModel = "diagonal" is required)
Y0 = scale((X %*% t(D) + matrix(rnorm(n * p), n, p)) %*% t(solve(diag(p) - A)))
Beta = SigmaHat = matrix(0, p, k)
for (i in 1:p) for (j in 1:k) {
  b = sum(X[, j] * Y0[, i]) / sum(X[, j]^2)                          # centered data, no intercept
  Beta[i, j] = b; SigmaHat[i, j] = sum((Y0[, i] - b * X[, j])^2) / n  # residual sum of squares / n
}
fit_b = RGM(Sxx = crossprod(X) / n, Beta = Beta, SigmaHat = SigmaHat, D = D, n = n,
            SigmaStarModel = "diagonal", nIter = 2000, nBurnin = 1000, Thin = 2)
round(fit_b$GammaEst, 2)