R/mr_radialegger_rjags.R
mr_radialegger_rjags.Rd
Bayesian radial MR-Egger model with a choice of prior distributions fitted using JAGS.
mr_radialegger_rjags( object, prior = "default", betaprior = "", sigmaprior = "", n.chains = 3, n.burn = 1000, n.iter = 5000, seed = NULL, rho = 0.5, ... )
object | A data object of class |
---|---|
prior | A character string for selecting the prior distributions;
|
betaprior | A character string in JAGS syntax to allow a user defined prior for the causal effect. |
sigmaprior | A character string in JAGS syntax to allow a user defined prior for the residual standard deviation. |
n.chains | Numeric indicating the number of chains used in the MCMC estimation, the default is |
n.burn | Numeric indicating the burn-in period of the Bayesian MCMC estimation. The default is |
n.iter | Numeric indicating the number of iterations in the Bayesian MCMC estimation. The default is |
seed | Numeric indicating the random number seed. The default is the rjags default. |
rho | Numeric indicating the correlation coefficient input into the joint prior distribution. The default is |
... | Additional arguments passed through to |
An object of class radialeggerjags
containing the following components:
The mean of the simulated pleiotropic effect
The mean of the simulated causal effect
Standard deviation of the simulated causal effect
The mean of the simaulted residual standard deviation
The credible interval for the causal effect, which includes the lower (2.5%), median (50%) and upper intervals (97.5%)
Output of the Bayesian MCMC samples
The specified priors
Bowden, J., et al., Improving the visualization, interpretation and analysis of two-sample summary data Mendelian randomization via the Radial plot and Radial regression. International Journal of Epidemiology, 2018. 47(4): p. 1264-1278. doi: 10.1093/ije/dyy101 .
if (requireNamespace("rjags", quietly = TRUE)) { fit <- mr_radialegger_rjags(bmi_insulin) summary(fit) plot(fit$samples) # 90% credible interval fitdf <- do.call(rbind.data.frame, fit$samples) cri90 <- quantile(fitdf$Estimate, probs = c(0.05,0.95)) print(cri90) }#> Warning: The mean of the sigma parameter, the residual standard deviation, is less than 1, we recommend refitting the model with sigma constrained to be >= 1.#> Prior : #> #> Pleiotropy ~ dnorm(0, 1E-3) #> Estimate ~ dnorm(0, 1E-3) #> sigma ~ dunif(.0001, 10) #> #> Estimation results: #> #> MCMC iterations = 6000 #> Burn in = 1000 #> Sample size by chain = 5000 #> Number of Chains = 3 #> Number of SNPs = 14 #> #> Inflating Parameter: 0.02679302 #> #> Estimate SD 2.5% 50% 97.5% #> Avg Pleio -27.288326 11.597682 -49.444400 -27.306644 -2.723596 #> Causal Effect 5.747213 2.222423 1.094817 5.757201 9.996314#> 5% 95% #> 1.974631 9.355917