Learn R Programming

RoBMA (version 1.2.0)

RoBMA: Estimate a Robust Bayesian Meta-Analysis

Description

RoBMA is used to estimate a Robust Bayesian Meta-Analysis. Either t-statistics (t) and sample sizes of the original studies (n or n1 and n2), or effect sizes (d) and standard errors (se) can be used to estimate the model.

Usage

RoBMA(
  t = NULL,
  d = NULL,
  r = NULL,
  y = NULL,
  OR = NULL,
  se = NULL,
  n = NULL,
  n1 = NULL,
  n2 = NULL,
  lCI = NULL,
  uCI = NULL,
  test_type = "two.sample",
  study_names = NULL,
  mu_transform = if (!is.null(r)) "cohens_d" else if (!is.null(OR)) "log_OR" else NULL,
  effect_direction = "positive",
  priors_mu = prior(distribution = "normal", parameters = list(mean = 0, sd = 1)),
  priors_tau = prior(distribution = "invgamma", parameters = list(shape = 1, scale =
    0.15)),
  priors_omega = list(prior(distribution = "two.sided", parameters = list(alpha = c(1,
    1), steps = c(0.05)), prior_odds = 1/2), prior(distribution = "two.sided", parameters
    = list(alpha = c(1, 1, 1), steps = c(0.05, 0.1)), prior_odds = 1/2)),
  priors_mu_null = prior(distribution = "point", parameters = list(location = 0)),
  priors_tau_null = prior(distribution = "point", parameters = list(location = 0)),
  priors_omega_null = prior(distribution = "point", parameters = list(location = 1)),
  chains = 3,
  iter = 10000,
  burnin = 5000,
  thin = 1,
  parallel = FALSE,
  control = NULL,
  save = "all",
  seed = NULL
)

Arguments

t

a vector of t-statistics.

d

a vector of effect sizes measured as Cohen's d.

r

a vector of effect sizes measured as correlations.

y

a vector of unspecified effect sizes.

OR

a vector of odds ratios.

se

a vector of standard errors of the effect sizes.

n

a vector of overall sample sizes.

n1

a vector of sample sizes for the first group.

n2

a vector of sample sizes for the second group.

lCI

a vector of lower bounds of confidence intervals.

uCI

a vector of upper bounds of confidence intervals.

test_type

a type of test used in the original studies. Options are "two.sample" (default) and "one.sample". Only available if d is supplied.

study_names

an optional argument with the names of the studies.

mu_transform

transformation to be applied to the supplied effect sizes before fitting the individual models. Only available if correlations or odds ratios are supplied as input. Defaults to "cohens_d" for correlations (another options is "fishers_z") and "log_OR" for odds ratios (another options is "cohens_d"). Note that priors are specified on the transformed scale and estimates are transformed back (apart from tau).

effect_direction

the expected direction of the effect. The one-sided selection sets the weights omega to 1 to significant results in the expted direction. Defaults to "positive" (another oprion is "negative").

priors_mu

list of prior distributions for the mu parameter that will be treated as belonging to the alternative hypothesis. Defaults to prior(distribution = "normal", parameters = list(mean = 0, sd = 1)).

priors_tau

list of prior distributions for the tau parameter that will be treated as belonging to the alternative hypothesis. Defaults to prior(distribution = "invgamma", parameters = list(shape = 1, scale = .15)).

priors_omega

list of prior weight functions for the omega parameter that will be treated as belonging to the alternative hypothesis. Defaults to list( prior(distribution = "two.sided", parameters = list(alpha = c(1, 1), steps = c(.05)), prior_odds = 1/2), prior(distribution = "two.sided", parameters = list(alpha = c(1, 1, 1), steps = c(.05, .10)), prior_odds = 1/2) ).

priors_mu_null

list of prior distributions for the mu parameter that will be treated as belonging to the null hypothesis. Defaults to point distribution with location at 0 ( prior(distribution = "point", parameters = list(location = 0))).

priors_tau_null

list of prior distributions for the tau parameter that will be treated as belonging to the null hypothesis. Defaults to point distribution with location at 0 ( prior(distribution = "point", parameters = list(location = 0))).

priors_omega_null

list of prior weight functions for the omega parameter that will be treated as belonging to the null hypothesis. Defaults to point distribution with location at 1 ( prior(distribution = "point", parameters = list(location = 0))).

chains

a number of chains of the MCMC algorithm.

iter

a number of sampling iterations of the MCMC algorithm. Defaults to 10000, with a minimum of 4000.

burnin

a number of burnin iterations of the MCMC algorithm. Defaults to 5000.

thin

a thinning of the chains of the MCMC algorithm. Defaults to 1.

parallel

whether the individual models should be fitted in parallel. Defaults to FALSE. The cores argument within the control list will overwrite the setting if specified to a number higher than 1.

control

a list of additional arguments for the MCMC algorithm. Possible options are:

autofit

Whether the models should be refitted until convergence. Defaults to FALSE

max_error

The target MCMC error for the autofit function. The argument is passed to raftery.diag as 'r'. Defaults to .01.

max_rhat

The target Rhat error for the autofit function. The argument is passed to add.summary as 'psrf.target'. Defaults to 1.05.

max_time

A string specifying the maximum fitting time in case of autofit. Defaults to Inf. Can be specified as a number and a unit (Acceptable units include <U+2019>seconds<U+2019>, <U+2019>minutes<U+2019>, <U+2019>hours<U+2019>, <U+2019>days<U+2019>, <U+2019>weeks<U+2019>, or the first letter(s) of each), i.e. "1hr".

adapt

A number of iterations used for MCMC adaptation. Defaults to 1000.

bridge_max_iter

Maximum number of iterations for the bridge_sampler function. Defaults to 10000

allow_max_error

Maximum allowed MCMC error for a model to be taken into consideration. The model will be removed from the ensemble if it fails to achieve the set MCMC error. Defaults to NULL - no model will be removed based on MCMC error.

allow_max_rhat

Maximum allowed Rhat for a model to be taken into consideration. Model will be removed from the ensemble if it fails to achieve the set Rhat. Defaults to NULL - no model will be removed based on Rhat.

allow_min_ESS

Minimum allowed ESS for a model to be taken into consideration. Model will be removed from the ensemble if it fails to achieve the set ESS. Defaults to NULL - no model will be removed based on ESS.

balance_prob

Whether the prior probability of the removed model should be redistributed to other models with the same type if possible (crossing of effect / heterogeneity / publication bias). Defaults to TRUE.

silent

Whether all fitting messages should be suppressed. Defaults to FALSE. Ideal for getting rid of the "full precision may not have been achieved in pntfinal'" warning that cannot be suppressed in any other way.

boost

Whether the likelihood functions implemented using the boost C++ library should be used as the first option. The higher precision of boost allows to estimate models in difficult cases. Defaults to FALSE. The R distributions are used as default and boost is used only if they fail. Warning: the estimation using boost takes about twice as long.

cores

Number of cores to bu used fo parallel computation. If parallel == TRUE, the default number is equal to number of cores - 1, and 1 (no parallel processing otherwise).

save

whether all models posterior distributions should be kept after obtaining a model-averaged result. Defaults to "all" which does not remove anything. Set to "min" to significantly reduce the size of final object, however, some model diagnostics check() will not be available.

seed

a seed to be set before model fitting, marginal likelihood computation, and posterior mixing for exact results reproducibility. Defaults to NULL - no seed is set.

Value

RoBMA returns an object of class "RoBMA".

Details

The default settings with either t-statistics / Cohen's d effect sizes and sample sizes / standard errors correspond to the ensemble proposed by maier2020RoBMA. The vignette("CustomEnsembles") and vignette("ReproducingBMA") vignettes describe how to use RoBMA() to fit custom meta-analytic ensembles (see prior() for more information about prior distributions). To get help with the error and warning messages, see vignette("WarningsAndErrors").

The RoBMA function first generates models from a combination of the provided priors for each of the model parameters. Then, the individual models are fitted using autorun.jags function. A marginal likelihood is computed using bridge_sampler function. The individual models are then combined into an ensemble using the posterior model probabilities.

Generic summary.RoBMA(), print.RoBMA(), and plot.RoBMA() functions are provided to facilitate manipulation with the ensemble. A visual check of the individual model diagnostics can be obtained using the diagnostics() function. The fitted model can be further updated or modified by update.RoBMA() function.

References

See Also

summary.RoBMA(), update.RoBMA(), prior(), check_setup()

Examples

Run this code
# NOT RUN {
# using the example data from Anderson et al. 2010 and fitting the default model
# (note that the model can take a while to fit)
fit <- RoBMA(r = Anderson2010$r, n = Anderson2010$n, study_names = Anderson2010$labels)

# in order to speed up the process, we can reduce the default number of chains, iteration,
# and disable the autofit functionality (see ?RoBMA for all possible settings)
fit_faster <- RoBMA(r = Anderson2010$r, n = Anderson2010$n, study_names = Anderson2010$labels,
chains = 2, iter = 5000, control = list(autofit = FALSE))

# RoBMA function allows to use different prior specifications
# for example, change the prior for tau to be half normal and specify one-sided selection only
# on significant p-values (see '?.prior' for all options regarding prior distributions)
fit1 <- RoBMA(r = Anderson2010$r, n = Anderson2010$n, study_names = Anderson2010$labels,
              priors_tau = prior("normal",
                                 parameters = list(mean = 0, sd = 1),
                                 truncation = list(lower = 0, upper = Inf)),
              priors_omega = prior("one-sided",
                                   parameters = list(cuts = c(.05), alpha = c(1, 1))))

# the priors for the null models can be modified or even omited in a similar manner,
# allowing to test different (non-nill-null) hypotheses
fit2 <- RoBMA(r = Anderson2010$r, n = Anderson2010$n, study_names = Anderson2010$labels,
              priors_mu_null  = prior("normal",
                                 parameters = list(mean = 0, sd = .1),
                                 truncation = list(lower = -0.1, upper = 0.1)))

# an already fitted RoBMA model can be further updated or modified by using the update function
# for example, the prior model probabilities can be changed after the fitting by
# (but see '?update.RoBMA' for other posibilities including refitting or adding more models)
fit3 <- update(fit2, prior_odds = c(10,1,1,1,1,1,1,1,1,1,1,1))

# we can get a quick overview of the model coefficients just by printing the model
fit

# a more detailed overview using the summary function (see '?summary.RoBMA' for all options)
summary(fit)

# results of the models can be visualized using the plot function (see ?plot.RoBMA for all options)
# for example, the model-averaged mean estimate
plot(fit, parameter = "mu")

# diagnostics for the individual parameters in individual models can be obtained using diagnostics
# function (see 'diagnostics' for all options)
diagnostics(fit, parameter = "mu", type = "chains")
# }
# NOT RUN {
# }

Run the code above in your browser using DataLab