diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml index 6cce515..46f531d 100644 --- a/.github/workflows/pkgdown.yaml +++ b/.github/workflows/pkgdown.yaml @@ -6,7 +6,7 @@ on: - main - master tags: - -'*' + - '*' jobs: build: @@ -47,12 +47,16 @@ jobs: - name: Setup renv uses: r-lib/actions/setup-renv@v2 - # Install the package + # Install Pandoc + - name: Setup Pandoc + uses: r-lib/actions/setup-pandoc@v2 + + # Install the package and its dependencies - name: Install devtools and the RoBMA Package run: | install.packages('devtools') install.packages('pkgdown') - install.packages(c('metaBMA', 'metafor', 'weightr', 'lme4', 'fixest', 'emmeans', 'metadat', 'vdiffr', 'testthat', 'covr', 'pandoc')) + install.packages(c('metaBMA', 'metafor', 'weightr', 'lme4', 'fixest', 'emmeans', 'metadat', 'vdiffr', 'testthat', 'covr')) devtools::install() shell: Rscript {0} diff --git a/_pkgdown.yml b/_pkgdown.yml new file mode 100644 index 0000000..15655d6 --- /dev/null +++ b/_pkgdown.yml @@ -0,0 +1,3 @@ +url: 'https://https://fbartos.github.io/RoBMA/' +template: + bootstrap: 5 \ No newline at end of file diff --git a/docs/404.html b/docs/404.html new file mode 100644 index 0000000..ff9091b --- /dev/null +++ b/docs/404.html @@ -0,0 +1,87 @@ + + +
+ + + + +vignettes/CustomEnsembles.Rmd
+ CustomEnsembles.Rmd
This vignette provides a step-by-step guide to fitting custom +meta-analytic ensembles using the RoBMA R package. By the end of this +guide, you will be able to construct and evaluate custom meta-analytic +models.
+By default, the RoBMA()
function specifies models as a
+combination of all supplied prior distributions (across null and
+alternative specification), with their prior model weights being equal
+to the product of prior distributions’ weights. This results in the 36
+meta-analytic models using the default settings (Bartoš et al.,
+2023).
+In another vignette, we illustrated
+that RoBMA can be also utilized for reproducing Bayesian Model-Averaged
+Meta-Analysis (BMA) (Bartoš et al., 2021; Gronau
+et al., 2017, 2021). However, the package was built as a
+framework for estimating highly customized meta-analytic model
+ensembles. Here, we are going to illustrate how to do exactly that (see
+Bartoš et al. (2022) for a tutorial paper
+on customizing the model ensemble with JASP).
Please keep in mind that all models should be justified by theory. +Furthermore, the models should be tested to make sure that the ensemble +can perform as intended a priori to drawing inference from it. The +following sections are only for illustrating the functionality of the +package. We provide a complete discussion with the relevant sources in +the Example section of Bartoš et al. +(2023).
+To illustrate the custom model building procedure, we use data from +the infamous Bem (2011) “Feeling the +future” precognition study. We use coding of the results as summarized +by Bem in one of his later replies (Bem et al., +2011).
+
+library(RoBMA)
+#> Loading required namespace: runjags
+#> Loading required namespace: mvtnorm
+
+data("Bem2011", package = "RoBMA")
+Bem2011
+#> d se study
+#> 1 0.25 0.10155048 Detection of Erotic Stimuli
+#> 2 0.20 0.08246211 Avoidance of Negative Stimuli
+#> 3 0.26 0.10323629 Retroactive Priming I
+#> 4 0.23 0.10182427 Retroactive Priming II
+#> 5 0.22 0.10120277 Retroactive Habituation I - Negative trials
+#> 6 0.15 0.08210765 Retroactive Habituation II - Negative trials
+#> 7 0.09 0.07085372 Retroactive Induction of Boredom
+#> 8 0.19 0.10089846 Facilitation of Recall I
+#> 9 0.42 0.14752627 Facilitation of Recall II
We consider the following scenarios as plausible explanations for the +data, and decide to include only those models into the meta-analytic +ensemble:
+If we were to fit the ensemble using the RoBMA()
+function and specifying all of the priors, we would have ended with 2
+(effect or no effect) * 2 (heterogeneity or no heterogeneity) * 5 (no
+publication bias or 4 ways of adjusting for publication bias) = 20
+models. That is 13 models more than requested. Furthermore, we could not
+specify different parameters for the prior distributions for each model.
+The following process allows this, though we do not utilize it here.
We start with fitting only the first model using the
+RoBMA()
function and we will continuously update the fitted
+object to include all of the models.
We initiate the model ensemble by specifying only the first model
+with the RoBMA()
function. We explicitly specify prior
+distributions for all components and set the prior distributions to
+correspond to the null hypotheses and set the seed to ensure
+reproducibility of the results.
+fit <- RoBMA(d = Bem2011$d, se = Bem2011$se, study_names = Bem2011$study,
+ priors_effect = NULL, priors_heterogeneity = NULL, priors_bias = NULL,
+ priors_effect_null = prior("spike", parameters = list(location = 0)),
+ priors_heterogeneity_null = prior("spike", parameters = list(location = 0)),
+ priors_bias_null = prior_none(),
+ seed = 1)
We verify that the ensemble contains only the single specified model
+with the summary()
function by setting
+type = "models"
.
+summary(fit, type = "models")
+#> Call:
+#> RoBMA(d = Bem2011$d, se = Bem2011$se, study_names = Bem2011$study,
+#> priors_effect = NULL, priors_heterogeneity = NULL, priors_bias = NULL,
+#> priors_effect_null = prior("spike", parameters = list(location = 0)),
+#> priors_heterogeneity_null = prior("spike", parameters = list(location = 0)),
+#> priors_bias_null = prior_none(), seed = 1)
+#>
+#> Robust Bayesian meta-analysis
+#> Models overview:
+#> Model Prior Effect Prior Heterogeneity Prior prob. log(marglik) Post. prob.
+#> 1 Spike(0) Spike(0) 1.000 -3.28 1.000
+#> Inclusion BF
+#> Inf
Before we add the second model to the ensemble, we need to decide on
+the prior distribution for the mean parameter. If precognition were to
+exist, the effect would be small since all casinos would be bankrupted
+otherwise. The effect would also be positive, since any deviation from
+randomness could be characterized as an effect. Therefore, we decide to
+use a normal distribution with mean = 0.15, standard deviation 0.10, and
+truncated to the positive range. This sets the prior density around
+small effect sizes. To get a better grasp of the prior distribution, we
+visualize it using the plot())
function (the figure can
+also be created using the ggplot2 package by adding
+plot_type = "ggplot"
argument).
We add the second model to the ensemble using the
+update.RoBMA()
function. The function can also be used for
+many other purposes - updating settings, prior model weights, and
+refitting failed models. Here, we supply the fitted ensemble object and
+add an argument specifying the prior distributions of each component for
+the additional model. Since we want to add Model 2 - we set the prior
+for the
+
+parameter to be treated as a prior belonging to the alternative
+hypothesis of the effect size component and the remaining priors treated
+as belonging to the null hypotheses. If we wanted, we could also specify
+prior_weights
argument, to change the prior probability of
+the fitted model but we do not utilize this option here and keep the
+default value, which sets the prior weights for the new model to
+1
. (Note that the arguments for specifying prior
+distributions in update.RoBMA()
function are
+prior_X
- in singular, in comparison to
+RoBMA()
function that uses priors_X
in
+plural.)
+fit <- update(fit,
+ prior_effect = prior("normal", parameters = list(mean = .15, sd = .10), truncation = list(lower = 0)),
+ prior_heterogeneity_null = prior("spike", parameters = list(location = 0)),
+ prior_bias_null = prior_none())
We can again inspect the updated ensemble to verify that it contains +both models. We see that Model 2 notably outperformed the first model +and attained all of the posterior model probability.
+
+summary(fit, type = "models")
+#> Call:
+#> RoBMA(d = Bem2011$d, se = Bem2011$se, study_names = Bem2011$study,
+#> priors_effect = NULL, priors_heterogeneity = NULL, priors_bias = NULL,
+#> priors_effect_null = prior("spike", parameters = list(location = 0)),
+#> priors_heterogeneity_null = prior("spike", parameters = list(location = 0)),
+#> priors_bias_null = prior_none(), seed = 1)
+#>
+#> Robust Bayesian meta-analysis
+#> Models overview:
+#> Model Prior Effect Prior Heterogeneity Prior prob. log(marglik)
+#> 1 Spike(0) Spike(0) 0.500 -3.28
+#> 2 Normal(0.15, 0.1)[0, Inf] Spike(0) 0.500 14.91
+#> Post. prob. Inclusion BF
+#> 0.000 0.000
+#> 1.000 79422247.251
Before we add the remaining models to the ensemble using the
+update()
function, we need to decide on the remaining prior
+distributions. Specifically, on the prior distribution for the
+heterogeneity parameter
+,
+and the publication bias adjustment parameters
+
+(for the selection models’ weightfunctions) and PET and PEESE for the
+PET and PEESE adjustment.
For Model 3, we use the usual inverse-gamma(1, .15) prior +distribution based on empirical heterogeneity estimates (Erp et al., 2017) for the heterogeneity +parameter +. +For Models 4.1-4.4 we use the default settings for the publication bias +adjustments as outlined in the Appendix B of (Bartoš et al., 2023).
+Now, we just need to add the remaining models to the ensemble using
+the update()
function as already illustrated.
+### adding Model 3
+fit <- update(fit,
+ prior_effect = prior("normal", parameters = list(mean = .15, sd = .10), truncation = list(lower = 0)),
+ prior_heterogeneity = prior("invgamma", parameters = list(shape = 1, scale = .15)),
+ prior_bias_null = prior_none())
+
+### adding Model 4.1
+fit <- update(fit,
+ prior_effect_null = prior("spike", parameters = list(location = 0)),
+ prior_heterogeneity_null = prior("spike", parameters = list(location = 0)),
+ prior_bias = prior_weightfunction("one.sided", parameters = list(alpha = c(1, 1), steps = c(0.05))))
+
+### adding Model 4.2
+fit <- update(fit,
+ prior_effect_null = prior("spike", parameters = list(location = 0)),
+ prior_heterogeneity_null = prior("spike", parameters = list(location = 0)),
+ prior_bias = prior_weightfunction("one.sided", parameters = list(alpha = c(1, 1, 1), steps = c(0.05, 0.10))))
+
+### adding Model 4.3
+fit <- update(fit,
+ prior_effect_null = prior("spike", parameters = list(location = 0)),
+ prior_heterogeneity_null = prior("spike", parameters = list(location = 0)),
+ prior_bias = prior_PET("Cauchy", parameters = list(0, 1), truncation = list(lower = 0)))
+
+### adding Model 4.4
+fit <- update(fit,
+ prior_effect_null = prior("spike", parameters = list(location = 0)),
+ prior_heterogeneity_null = prior("spike", parameters = list(location = 0)),
+ prior_bias = prior_PEESE("Cauchy", parameters = list(0, 5), truncation = list(lower = 0)))
We again verify that all of the requested models are included in the
+ensemble using the summary())
function with
+type = "models"
argument.
+summary(fit, type = "models")
+#> Call:
+#> RoBMA(d = Bem2011$d, se = Bem2011$se, study_names = Bem2011$study,
+#> priors_effect = NULL, priors_heterogeneity = NULL, priors_bias = NULL,
+#> priors_effect_null = prior("spike", parameters = list(location = 0)),
+#> priors_heterogeneity_null = prior("spike", parameters = list(location = 0)),
+#> priors_bias_null = prior_none(), seed = 1)
+#>
+#> Robust Bayesian meta-analysis
+#> Models overview:
+#> Model Prior Effect Prior Heterogeneity
+#> 1 Spike(0) Spike(0)
+#> 2 Normal(0.15, 0.1)[0, Inf] Spike(0)
+#> 3 Normal(0.15, 0.1)[0, Inf] InvGamma(1, 0.15)
+#> 4 Spike(0) Spike(0)
+#> 5 Spike(0) Spike(0)
+#> 6 Spike(0) Spike(0)
+#> 7 Spike(0) Spike(0)
+#> Prior Bias Prior prob. log(marglik)
+#> 0.143 -3.28
+#> 0.143 14.91
+#> 0.143 12.85
+#> omega[one-sided: .05] ~ CumDirichlet(1, 1) 0.143 13.70
+#> omega[one-sided: .1, .05] ~ CumDirichlet(1, 1, 1) 0.143 12.58
+#> PET ~ Cauchy(0, 1)[0, Inf] 0.143 15.75
+#> PEESE ~ Cauchy(0, 5)[0, Inf] 0.143 15.65
+#> Post. prob. Inclusion BF
+#> 0.000 0.000
+#> 0.168 1.210
+#> 0.021 0.132
+#> 0.050 0.318
+#> 0.016 0.100
+#> 0.391 3.845
+#> 0.353 3.278
Finally, we use the summary()
function to inspect the
+model results. The results from our custom ensemble indicate weak
+evidence for the absence of the precognition effect,
+
+->
+,
+moderate evidence for the absence of heterogeneity,
+
+->
+,
+and moderate evidence for the presence of the publication bias,
+.
+summary(fit)
+#> Call:
+#> RoBMA(d = Bem2011$d, se = Bem2011$se, study_names = Bem2011$study,
+#> priors_effect = NULL, priors_heterogeneity = NULL, priors_bias = NULL,
+#> priors_effect_null = prior("spike", parameters = list(location = 0)),
+#> priors_heterogeneity_null = prior("spike", parameters = list(location = 0)),
+#> priors_bias_null = prior_none(), seed = 1)
+#>
+#> Robust Bayesian meta-analysis
+#> Components summary:
+#> Models Prior prob. Post. prob. Inclusion BF
+#> Effect 2/7 0.286 0.189 0.584
+#> Heterogeneity 1/7 0.143 0.021 0.132
+#> Bias 4/7 0.571 0.811 3.212
+#>
+#> Model-averaged estimates:
+#> Mean Median 0.025 0.975
+#> mu 0.036 0.000 0.000 0.226
+#> tau 0.002 0.000 0.000 0.000
+#> omega[0,0.05] 1.000 1.000 1.000 1.000
+#> omega[0.05,0.1] 0.938 1.000 0.014 1.000
+#> omega[0.1,1] 0.935 1.000 0.012 1.000
+#> PET 0.820 0.000 0.000 2.601
+#> PEESE 7.284 0.000 0.000 25.508
+#> The estimates are summarized on the Cohen's d scale (priors were specified on the Cohen's d scale).
The finalized ensemble can be treated as any other RoBMA
+ensemble using the summary()
, plot()
,
+plot_models()
, forest()
, and
+diagnostics()
functions. For example, we can use the
+plot.RoBMA()
with the
+parameter = "mu", prior = TRUE
arguments to plot the prior
+(grey) and posterior distribution (black) for the effect size. The
+function visualizes the model-averaged estimates across all models by
+default. The arrows represent the probability mass at the value 0 (a
+spike at 0). The secondary y-axis (right) shows the probability mass at
+the zero effect size, which increased from the prior probability of 0.71
+to the posterior the posterior probability of 0.81.
+plot(fit, parameter = "mu", prior = TRUE)
We can also inspect the posterior distributions of the publication
+bias adjustments. To visualize the model-averaged weightfunction, we set
+parameter = weightfunction
argument. The resulting figure
+shows the light gray prior distribution and the dark gray the posterior
+distribution.
+plot(fit, parameter = "weightfunction", prior = TRUE)
We can also inspect the posterior estimate of the regression
+relationship between the standard errors and effect sizes by setting
+parameter = "PET-PEESE"
.
+plot(fit, parameter = "PET-PEESE", prior = TRUE)
+- The default setting used to produce 12 models in RoBMA versions < +2, which corresponded to an earlier an article by Maier et al. (2023) in which we applied Bayesian +model-averaging only across selection models.
+vignettes/HierarchicalBMA.Rmd
+ HierarchicalBMA.Rmd
Hierarchical (or multilevel/3-level) meta-analysis adjusts for the +dependency of effect sizes due to clustering in the data. For example, +effect size estimates from multiple experiments reported in the same +manuscript might be expected to be more similar than effect sizes from a +different paper (Konstantopoulos, 2011). +This vignette illustrates how to deal with such dependencies among +effect size estimates (in cases with simple nested structure) using the +Bayesian model-averaged meta-analysis (BMA) (Bartoš et al., 2021; Gronau et al., 2017, +2021). (See other vignettes for more details on BMA: Reproducing BMA or Informed BMA in medicine.)
+First, we introduce the example data set. Second, we illustrate the
+frequentist hierarchical meta-analysis with the metafor
R
+package and discuss the results. Third, we outline the hierarchical
+meta-analysis parameterization. Fourth, we estimate the Bayesian
+model-averaged hierarchical meta-analysis. Finally, we conclude by
+discussing further extensions and publication bias adjustment.
We use the dat.konstantopoulos2011
data set from the
+metadat
R package (Thomas et al.,
+2019) that is used for the same functionality in the metafor
+(Wolfgang, 2010) R package. We roughly
+follow the example in the data set’s help file,
+?dat.konstantopoulos2011
. The data set consists of 56
+studies estimating the effects of modified school calendars on students’
+achievement. The 56 studies were run in individual schools, which can be
+grouped into 11 districts. We might expect more similar effect size
+estimates from schools in the same district – in other words, the effect
+size estimates from the same district might not be completely
+independent. Consequently, we might want to adjust for this dependency
+(clustering) between the effect size estimates to draw a more
+appropriate inference.
First, we load the data set, assign it to the dat
+object, and inspect the first few rows.
+data("dat.konstantopoulos2011", package = "metadat")
+dat <- dat.konstantopoulos2011
+
+head(dat)
+#> district school study year yi vi
+#> 1 11 1 1 1976 -0.18 0.118
+#> 2 11 2 2 1976 -0.22 0.118
+#> 3 11 3 3 1976 0.23 0.144
+#> 4 11 4 4 1976 -0.30 0.144
+#> 5 12 1 5 1989 0.13 0.014
+#> 6 12 2 6 1989 -0.26 0.014
In the following analyses, we use the following variables:
+yi
, standardized mean differences,vi
, sampling variances of the standardized mean
+differences,district
, district id which distinguishes among the
+districts,school
, that distinguishes among different schools
+within the same district.metafor
+We follow the data set’s help file and fit a simple random effects
+meta-analysis using the rma()
function from
+metafor
package. This model ignores the dependency between
+effect size estimates. We use this simple model as our starting point
+and as a comparison with the later models.
+fit_metafor.0 <- metafor::rma(yi = yi, vi = vi, data = dat)
+fit_metafor.0
+#>
+#> Random-Effects Model (k = 56; tau^2 estimator: REML)
+#>
+#> tau^2 (estimated amount of total heterogeneity): 0.0884 (SE = 0.0202)
+#> tau (square root of estimated tau^2 value): 0.2974
+#> I^2 (total heterogeneity / total variability): 94.70%
+#> H^2 (total variability / sampling variability): 18.89
+#>
+#> Test for Heterogeneity:
+#> Q(df = 55) = 578.8640, p-val < .0001
+#>
+#> Model Results:
+#>
+#> estimate se zval pval ci.lb ci.ub
+#> 0.1279 0.0439 2.9161 0.0035 0.0419 0.2139 **
+#>
+#> ---
+#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The model summary returns a small but statistically significant +effect size estimate + +() +and a considerable heterogeneity estimate +.
+We extend the model to account for the hierarchical structure of the
+data, i.e., schools within districts, by using the rma.mv()
+function from the metafor
package and extending it with the
+random = ~ school | district
argument.
+fit_metafor <- metafor::rma.mv(yi, vi, random = ~ school | district, data = dat)
+fit_metafor
+#>
+#> Multivariate Meta-Analysis Model (k = 56; method: REML)
+#>
+#> Variance Components:
+#>
+#> outer factor: district (nlvls = 11)
+#> inner factor: school (nlvls = 11)
+#>
+#> estim sqrt fixed
+#> tau^2 0.0978 0.3127 no
+#> rho 0.6653 no
+#>
+#> Test for Heterogeneity:
+#> Q(df = 55) = 578.8640, p-val < .0001
+#>
+#> Model Results:
+#>
+#> estimate se zval pval ci.lb ci.ub
+#> 0.1847 0.0846 2.1845 0.0289 0.0190 0.3504 *
+#>
+#> ---
+#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
We find that accounting for the hierarchical structure of the data
+results in (1) a slightly larger effect size estimate
+()
+and (2) larger standard error of the effect size estimate
+().
+The larger standard error is a natural consequence of accounting for the
+dependency between the effect sizes. Because the effect sizes are
+dependent, they contribute less additional information than independent
+effect sizes would. Specifying the hierarchical model then accounts for
+the dependency by estimating similarity between the estimates from the
+same cluster (school) and discounting the information borrowed from each
+estimate. The estimate of the similarity among estimates from the same
+cluster is summarized in the \rho = 0.666
estimate.
We specify a simple hierarchical meta-analytic model (see Konstantopoulos (2011) for an example). Using +distributional notation, we can describe the data generating process as +a multi-stage sampling procedure. In a nutshell, we assume the existence +of an overall mean effect +. +Next, we assume that the effect sizes in each district +, +, +systematically differ from the mean effect, with the variance of the +district-level effects summarized with heterogeneity + +(as between). Furthermore, we assume that the true effects + +of each study + +systematically differ from the district-level effect, with the variance +of the study effects from the district-level effect summarized with +heterogeneity + +(as within). Finally, the observed effect sizes + +that differ from the true effects + +due to random errors +.
+Mathematically, we can describe such a model as: + Where N() denotes a normal distribution +with mean and variance.
+Conveniently, and with a bit of algebra, we do not need to estimate +the district-level and true study effects. Instead, we marginalize them +out, and we sample the observed effect sizes from each district + +directly from a multivariate normal distributions, MN(), with a common +mean + +and covariance matrix S: + The random effects marginalization is +helpful as it allows us to sample fewer parameters from the posterior +distribution (which significantly simplifies marginal likelihood +estimation via bridge sampling). Furthermore, the marginalization allows +us to properly specify selection model publication bias adjustment +models – the marginalization propagates the selection process up through +all the sampling steps at once (we cannot proceed with the sequential +sampling as the selection procedure on the observed effect sizes +modifies the sampling distributions of all the preceding levels).
+We can further re-parameterize the model by performing the following
+substitution,
+ and specifying the covariance matrix
+using the inter-study correlation
+,
+total heterogeneity
+,
+and the standard errors
+:
+ This specification corresponds to the
+compound symmetry covariance matrix of random effects, the default
+settings in the metafor::rma.mv()
function. More
+importantly, it allows us to easily specify prior distributions on the
+correlation coefficient
+
+and the total heterogeneity
+.
RoBMA
+Before we estimate the complete Hierarchical Bayesian Model-Averaged
+Meta-Analysis (hBMA) with the RoBMA
package, we briefly
+reproduce the simpler models we estimated with the metafor
+package in the previous section.
First, we estimate a simple Bayesian random effects meta-analysis
+(corresponding to fit_metafor.0
). We use
+the RoBMA()
function and specify the effect sizes and
+sampling variances via the d = dat$yi
and
+v = dat$vi
arguments. We set the
+priors_effect_null
, priors_heterogeneity_null
,
+and priors_bias
arguments to null to omit models assuming
+the absence of the effect, heterogeneity, and the publication bias
+adjustment components.
+fit.0 <- RoBMA(d = dat$yi, v = dat$vi,
+ priors_effect_null = NULL,
+ priors_heterogeneity_null = NULL,
+ priors_bias = NULL,
+ parallel = TRUE, seed = 1)
We generate a complete summary for the only estimated model by adding
+the type = "individual"
argument to the
+summary()
function.
+summary(fit.0, type = "individual")
+#> Call:
+#> RoBMA(d = dat$yi, v = dat$vi, priors_bias = NULL, priors_effect_null = NULL,
+#> priors_heterogeneity_null = NULL, parallel = TRUE, seed = 1)
+#>
+#> Robust Bayesian meta-analysis
+#> Model 1 Parameter prior distributions
+#> Prior prob. 1.000 mu ~ Normal(0, 1)
+#> log(marglik) 17.67 tau ~ InvGamma(1, 0.15)
+#> Post. prob. 1.000
+#> Inclusion BF Inf
+#>
+#> Parameter estimates:
+#> Mean SD lCI Median uCI error(MCMC) error(MCMC)/SD ESS R-hat
+#> mu 0.126 0.043 0.041 0.127 0.211 0.00044 0.010 9757 1.000
+#> tau 0.292 0.033 0.233 0.290 0.364 0.00034 0.010 9678 1.000
+#> The estimates are summarized on the Cohen's d scale (priors were specified on the Cohen's d scale).
We verify that the effect size, + +(), +and heterogeneity, + +(), +estimates closely correspond to the frequentist results (as we would +expect from parameter estimates under weakly informative priors).
+Second, we account for the clustered effect size estimates within
+districts by extending the previous function call with the
+study_ids = dat$district
argument. This allows us to
+estimate the hierarchical Bayesian random effects meta-analysis
+(corresponding to fit_metafor
). We use the default prior
+distribution for the correlation parameter
+\rho \sim \text{Beta}(1, 1)
, set via the
+priors_hierarchical
argument, which restricts the
+correlation to be positive and uniformly distributed on the interval
+.
+fit <- RoBMA(d = dat$yi, v = dat$vi, study_ids = dat$district,
+ priors_effect_null = NULL,
+ priors_heterogeneity_null = NULL,
+ priors_bias = NULL,
+ parallel = TRUE, seed = 1)
Again, we generate the complete summary for the only estimated +model,
+
+summary(fit, type = "individual")
+#> Call:
+#> RoBMA(d = dat$yi, v = dat$vi, study_ids = dat$district, priors_bias = NULL,
+#> priors_effect_null = NULL, priors_heterogeneity_null = NULL,
+#> parallel = TRUE, seed = 1)
+#>
+#> Robust Bayesian meta-analysis
+#> Model 1 Parameter prior distributions
+#> Prior prob. 1.000 mu ~ Normal(0, 1)
+#> log(marglik) 25.70 tau ~ InvGamma(1, 0.15)
+#> Post. prob. 1.000 rho ~ Beta(1, 1)
+#> Inclusion BF Inf
+#>
+#> Parameter estimates:
+#> Mean SD lCI Median uCI error(MCMC) error(MCMC)/SD ESS R-hat
+#> mu 0.181 0.083 0.017 0.180 0.346 0.00088 0.011 9041 1.000
+#> tau 0.308 0.056 0.223 0.299 0.442 0.00090 0.016 3859 1.000
+#> rho 0.627 0.142 0.320 0.641 0.864 0.00219 0.015 4202 1.000
+#> The estimates are summarized on the Cohen's d scale (priors were specified on the Cohen's d scale).
and verify that our estimates, again, correspond to the frequentist +counterparts, with the estimated effect size, + +(), +heterogeneity, + +(), +and correlation, + +().
+We can further visualize the prior and posterior distribution of the
+
+parameter using the plot()
function.
Third, we extend the previous model into a model ensemble that also
+includes models assuming the absence of the effect and/or heterogeneity
+(we do not incorporate models assuming presence of publication bias due
+to computational complexity explained in the summary). Including those
+additional models allows us to evaluate evidence in favor of the effect
+and heterogeneity. Furthermore, specifying all those additional models
+allows us to incorporate the uncertainty about the specified models and
+weight the posterior distribution according to how well the models
+predicted the data. We estimate the remaining models by removing the
+priors_effect_null
and
+priors_heterogeneity_null
arguments from the previous
+function calls, which include the previously omitted models of no effect
+and/or no heterogeneity.
+fit_BMA <- RoBMA(d = dat$yi, v = dat$vi, study_ids = dat$district,
+ priors_bias = NULL,
+ parallel = TRUE, seed = 1)
Now we generate a summary for the complete model-averaged ensemble by
+not specifying any additional arguments in the summary()
+function.
+summary(fit_BMA)
+#> Call:
+#> RoBMA(d = dat$yi, v = dat$vi, study_ids = dat$district, priors_bias = NULL,
+#> parallel = TRUE, seed = 1)
+#>
+#> Robust Bayesian meta-analysis
+#> Components summary:
+#> Models Prior prob. Post. prob. Inclusion BF
+#> Effect 2/4 0.500 0.478 9.170000e-01
+#> Heterogeneity 2/4 0.500 1.000 9.326943e+92
+#> Hierarchical 2/4 0.500 1.000 9.326943e+92
+#>
+#> Model-averaged estimates:
+#> Mean Median 0.025 0.975
+#> mu 0.087 0.000 0.000 0.314
+#> tau 0.326 0.317 0.231 0.472
+#> rho 0.659 0.675 0.354 0.879
+#> The estimates are summarized on the Cohen's d scale (priors were specified on the Cohen's d scale).
We find the ensemble contains four models, the combination of models
+assuming the presence/absence of the effect/heterogeneity, each with
+equal prior model probabilities. Importantly, the models assuming
+heterogeneity are also specified with the hierarchical structure and
+account for the clustering. A comparison of the specified models reveals
+weak evidence against the effect,
+,
+and extreme evidence for the presence of heterogeneity,
+.
+Moreover, we find that the Hierarchical
component summary
+has the same values as the Heterogeneity
component summary
+because the default settings specify that all models assuming the
+presence of heterogeneity also include the hierarchical structure.
We also obtain the model-averaged posterior estimates that combine +the posterior estimates from all models according to the posterior model +probabilities, the effect size, + +(), +heterogeneity, + +(), +and correlation, + +().
+In the previous analyses, we assumed that the effect sizes are indeed +clustered within the districts, and we only adjusted for the clustering. +However, the effect sizes within the same cluster may not be more +similar than effect sizes from different clusters. Now, we specify a +model ensemble that allows us to test this assumption by specifying two +sets of random effect meta-analytic models. The first set of models +assumes that there is indeed clustering and that the correlation of +random effects is uniformly distributed on the + +interval (as in the previous analyses). The second set of models assumes +that there is no clustering, i.e., the correlation of random effects +, +which simplifies the structured covariance matrix to a diagonal matrix. +Again, we model average across models assuming the presence and absence +of the effect to account for the model uncertainty.
+To specify this ‘special’ model ensemble with the
+RoBMA()
function, we need to modify the previous model call
+in the following ways. We removed the fixed effect models by specifying
+the priors_heterogeneity_null = NULL
+argument.
+Furthermore, we specify the prior distribution for models assuming the
+absence of the hierarchical structure by adding the
+priors_hierarchical_null = prior(distribution = "spike", parameters = list("location" = 0))
+argument.
+hierarchical_test <- RoBMA(d = dat$yi, v = dat$vi, study_ids = dat$district,
+ priors_heterogeneity_null = NULL,
+ priors_hierarchical_null = prior(distribution = "spike", parameters = list("location" = 0)),
+ priors_bias = NULL,
+ parallel = TRUE, seed = 1)
+summary(hierarchical_test)
+#> Call:
+#> RoBMA(d = dat$yi, v = dat$vi, study_ids = dat$district, priors_bias = NULL,
+#> priors_heterogeneity_null = NULL, priors_hierarchical_null = prior(distribution = "spike",
+#> parameters = list(location = 0)), parallel = TRUE, seed = 1)
+#>
+#> Robust Bayesian meta-analysis
+#> Components summary:
+#> Models Prior prob. Post. prob. Inclusion BF
+#> Effect 2/4 0.500 0.478 0.917
+#> Heterogeneity 4/4 1.000 1.000 Inf
+#> Hierarchical 2/4 0.500 1.000 4624.794
+#>
+#> Model-averaged estimates:
+#> Mean Median 0.025 0.975
+#> mu 0.087 0.000 0.000 0.314
+#> tau 0.326 0.317 0.231 0.472
+#> rho 0.659 0.675 0.354 0.879
+#> The estimates are summarized on the Cohen's d scale (priors were specified on the Cohen's d scale).
We summarize the resulting model ensemble and find out that the
+Hierarchical
component is no longer equivalent to the
+Heterogeneity
component – the new model specification
+allowed us to compare random effect models assuming the presence of the
+hierarchical structure to random effect models assuming the absence of
+the hierarchical structure. The resulting inclusion Bayes factor of the
+hierarchical structure shows extreme evidence in favor of clustering of
+the effect sizes,
+,
+i.e., there is extreme evidence that the intervention results in more
+similar effects within the districts.
We illustrated how to estimate a hierarchical Bayesian model-averaged
+meta-analysis using the RoBMA
package. The hBMA model
+allows us to test for the presence vs absence of the effect and
+heterogeneity while simultaneously adjusting for clustered effect size
+estimates. While the current implementation allows us to draw a fully
+Bayesian inference, incorporate prior information, and acknowledge model
+uncertainty, it has a few limitations in contrast to the
+metafor
package. E.g., the RoBMA
package only
+allows a simple nested random effects (i.e., estimates within studies,
+schools within districts etc). The simple nesting allows us to break the
+full covariance matrix into per cluster block matrices which speeds up
+the already demanding computation. Furthermore, the computational
+complexity significantly increases when considering selection models as
+we need to compute an exponentially increasing number of multivariate
+normal probabilities with the increasing cluster size (existence of
+clusters with more than four studies makes the current implementation
+impractical due to the computational demands). However, these current
+limitations are not the end of the road, as we are exploring other
+approaches (e.g., only specifying PET-PEESE style publication bias
+adjustment and other dependency adjustments) in a future vignette.
+We could also model-average across the hierarchical structure assuming +fixed effect models, i.e., + +and +. +However specifying such a model ensemble is a beyond the scope of this +vignette, see Custom ensembles +vignette for some hints.
+vignettes/MedicineBMA.Rmd
+ MedicineBMA.Rmd
Bayesian model-averaged meta-analysis allows researchers to +seamlessly incorporate available prior information into the analysis +(Bartoš et al., 2021; Gronau et al., 2017, +2021). This vignette illustrates how to do this with an example +from Bartoš et al. (2021), who developed +informed prior distributions for meta-analyses of continuous outcomes +based on the Cochrane database of systematic reviews. Then, we extend +the example by incorporating publication bias adjustment with robust +Bayesian meta-analysis (Bartoš et al., 2023; +Maier et al., 2023).
+We illustrate how to fit the informed BMA (not adjusting for
+publication bias) using the RoBMA
R package. For this
+purpose, we reproduce the dentine hypersensitivity example from Bartoš et al. (2021), who reanalyzed five
+studies with a tactile outcome assessment that were subjected to a
+meta-analysis by Poulsen et al.
+(2006).
We load the dentine hypersensitivity data included in the +package.
+
+library(RoBMA)
+
+data("Poulsen2006", package = "RoBMA")
+Poulsen2006
+#> d se study
+#> 1 0.9073050 0.2720456 STD-Schiff-1994
+#> 2 0.7207589 0.1977769 STD-Silverman-1996
+#> 3 1.3305829 0.2721648 STD-Sowinski-2000
+#> 4 1.7688872 0.2656483 STD-Schiff-2000-(2)
+#> 5 1.3286828 0.3582617 STD-Schiff-1998
To reproduce the analysis from the example, we need to set informed +empirical prior distributions for the effect sizes +() +and heterogeneity +() +parameters that Bartoš et al. (2021) +obtained from the Cochrane database of systematic reviews. We can either +set them manually,
+
+fit_BMA <- RoBMA(d = Poulsen2006$d, se = Poulsen2006$se, study_names = Poulsen2006$study,
+ priors_effect = prior(distribution = "t", parameters = list(location = 0, scale = 0.51, df = 5)),
+ priors_heterogeneity = prior(distribution = "invgamma", parameters = list(shape = 1.79, scale = 0.28)),
+ priors_bias = NULL,
+ transformation = "cohens_d", seed = 1, parallel = TRUE)
with priors_effect
and priors_heterogeneity
+corresponding to the
+
+and
+
+informed prior distributions for the “oral health” subfield and removing
+the publication bias adjustment models by setting
+priors_bias = NULL
.
+Note that the package contains function NoBMA()
from
+version 3.1 which skips publication bias adjustment directly.
Alternatively, we can utilize the prior_informed
+function that prepares informed prior distributions for the individual
+medical subfields automatically.
+fit_BMA <- RoBMA(d = Poulsen2006$d, se = Poulsen2006$se, study_names = Poulsen2006$study,
+ priors_effect = prior_informed(name = "oral health", parameter = "effect", type = "smd"),
+ priors_heterogeneity = prior_informed(name = "oral health", parameter = "heterogeneity", type = "smd"),
+ priors_bias = NULL,
+ transformation = "cohens_d", seed = 1, parallel = TRUE)
The name
argument specifies the medical subfield name
+(use print(BayesTools::prior_informed_medicine_names)
to
+check names of all available subfields). The parameter
+argument specifies whether we want prior distribution for the effect
+size or heterogeneity. Finally, the type
argument specifies
+what type of measure we use for the meta-analysis (see
+?prior_informed
for more details regarding the informed
+prior distributions).
We obtain the output with the summary
function. Adding
+the conditional = TRUE
argument allows us to inspect the
+conditional estimates, i.e., the effect size estimate assuming that the
+models specifying the presence of the effect are true and the
+heterogeneity estimates assuming that the models specifying the presence
+of heterogeneity are
+true.
+summary(fit_BMA, conditional = TRUE)
+#> Call:
+#> RoBMA(d = Poulsen2006$d, se = Poulsen2006$se, study_names = Poulsen2006$study,
+#> transformation = "cohens_d", priors_effect = prior_informed(name = "oral health",
+#> parameter = "effect", type = "smd"), priors_heterogeneity = prior_informed(name = "oral health",
+#> parameter = "heterogeneity", type = "smd"), priors_bias = NULL,
+#> parallel = TRUE, seed = 1)
+#>
+#> Robust Bayesian meta-analysis
+#> Components summary:
+#> Models Prior prob. Post. prob. Inclusion BF
+#> Effect 2/4 0.500 0.995 217.517
+#> Heterogeneity 2/4 0.500 0.778 3.511
+#>
+#> Model-averaged estimates:
+#> Mean Median 0.025 0.975
+#> mu 1.076 1.088 0.664 1.422
+#> tau 0.231 0.208 0.000 0.726
+#> The estimates are summarized on the Cohen's d scale (priors were specified on the Cohen's d scale).
+#>
+#> Conditional estimates:
+#> Mean Median 0.025 0.975
+#> mu 1.082 1.090 0.701 1.422
+#> tau 0.297 0.255 0.075 0.779
+#> The estimates are summarized on the Cohen's d scale (priors were specified on the Cohen's d scale).
The output from the summary.RoBMA()
function has 3
+parts. The first one under the ‘Robust Bayesian Meta-Analysis’ heading
+provides a basic summary of the fitted models by component types
+(presence of the Effect and Heterogeneity). The table summarizes the
+prior and posterior probabilities and the inclusion Bayes factors of the
+individual components. The results show that the inclusion Bayes factor
+for the effect corresponds to the one reported in Bartoš et al. (2021),
+
+and
+
+(up to an MCMC error).
The second part under the ‘Model-averaged estimates’ heading displays +the parameter estimates model-averaged across all specified models +(i.e., including models specifying the effect size to be zero). We +ignore this section and move to the last part.
+The third part under the ‘Conditional estimates’ heading displays the +conditional effect size estimate corresponding to the one reported in +Bartoš et al. (2021), +, +95% CI [0.686,1.412], and a heterogeneity estimate (not reported +previously).
+The RoBMA
package provides extensive options for
+visualizing the results. Here, we visualize the prior (grey) and
+posterior (black) distribution for the mean parameter.
+plot(fit_BMA, parameter = "mu", prior = TRUE)
By default, the function plots the model-averaged estimates across
+all models; the arrows represent the probability of a spike, and the
+lines represent the posterior density under models assuming non-zero
+effect. The secondary y-axis (right) shows the probability of the spike
+(at value 0) decreasing from 0.50, to 0.005 (also obtainable from the
+‘Robust Bayesian Meta-Analysis’ field in summary.RoBMA()
+function).
To visualize the conditional effect size estimate, we can add the
+conditional = TRUE
argument,
+plot(fit_BMA, parameter = "mu", prior = TRUE, conditional = TRUE)
+which displays only the model-averaged posterior distribution of the +effect size parameter for models assuming the presence of the +effect.
+We can also visualize the estimates from the individual models used
+in the ensemble. We do that with the plot_models()
+function, which visualizes the effect size estimates and 95% CI of each
+specified model included in the ensemble. Model 1 corresponds to the
+fixed effect model assuming the absence of the effect,
+,
+Model 2 corresponds to the random effect model assuming the absence of
+the effect,
+,
+Model 3 corresponds to the fixed effect model assuming the presence of
+the effect,
+,
+and Model 4 corresponds to the random effect model assuming the presence
+of the effect,
+).
+The size of the square representing the mean estimate reflects the
+posterior model probability of the model, which is also displayed in the
+right-hand side panel. The bottom part of the figure shows the model
+averaged-estimate that is a combination of the individual model
+posterior distributions weighted by the posterior model
+probabilities.
+plot_models(fit_BMA)
+We see that the posterior model probability of the first two models +decreased to essentially zero (when rounding to two decimals), +completely omitting their estimates from the figure. Furthermore, the +much larger box of Model 4 (the random effect model assuming the +presence of the effect) shows that Model 4 received the largest share of +the posterior probability, +)
+The last type of visualization that we show here is the forest plot.
+It displays the original studies’ effects and the meta-analytic estimate
+within one figure. It can be requested by using the
+forest()
function. Here, we again set the
+conditional = TRUE
argument to display the conditional
+model-averaged effect size estimate at the bottom.
+forest(fit_BMA, conditional = TRUE)
For more options provided by the plotting function, see its
+documentation by using ?plot.RoBMA()
,
+?plot_models()
, and ?forest()
.
Finally, we illustrate how to adjust our informed BMA for publication +bias with robust Bayesian meta-analysis Maier et +al. (2023). In short, we specify additional models assuming the +presence of the publication bias and correcting for it by either +specifying a selection model operating on +-values +(Vevea & Hedges, 1995) or by +specifying a publication bias adjustment method correcting for the +relationship between effect sizes and standard errors – PET-PEESE (Stanley, 2017; Stanley & Doucouliagos, +2014). See Bartoš et al. (2022) for +a tutorial.
+To obtain a proper before and after publication bias adjustment +comparison, we fit the informed BMA model but using the default effect +size transformation (Fisher’s +).
+
+fit_BMAb <- RoBMA(d = Poulsen2006$d, se = Poulsen2006$se, study_names = Poulsen2006$study,
+ priors_effect = prior_informed(name = "oral health", parameter = "effect", type = "smd"),
+ priors_heterogeneity = prior_informed(name = "oral health", parameter = "heterogeneity", type = "smd"),
+ priors_bias = NULL,
+ seed = 1, parallel = TRUE)
+summary(fit_BMAb, conditional = TRUE)
+#> Call:
+#> RoBMA(d = Poulsen2006$d, se = Poulsen2006$se, study_names = Poulsen2006$study,
+#> priors_effect = prior_informed(name = "oral health", parameter = "effect",
+#> type = "smd"), priors_heterogeneity = prior_informed(name = "oral health",
+#> parameter = "heterogeneity", type = "smd"), priors_bias = NULL,
+#> parallel = TRUE, seed = 1)
+#>
+#> Robust Bayesian meta-analysis
+#> Components summary:
+#> Models Prior prob. Post. prob. Inclusion BF
+#> Effect 2/4 0.500 0.997 347.932
+#> Heterogeneity 2/4 0.500 0.723 2.608
+#>
+#> Model-averaged estimates:
+#> Mean Median 0.025 0.975
+#> mu 1.045 1.052 0.705 1.344
+#> tau 0.186 0.163 0.000 0.623
+#> The estimates are summarized on the Cohen's d scale (priors were specified on the Cohen's d scale).
+#>
+#> Conditional estimates:
+#> Mean Median 0.025 0.975
+#> mu 1.048 1.053 0.720 1.344
+#> tau 0.256 0.220 0.064 0.681
+#> The estimates are summarized on the Cohen's d scale (priors were specified on the Cohen's d scale).
We obtain noticeably stronger evidence for the presence of the +effect. This is a result of placing more weights on the fixed-effect +models, especially the fixed-effect model assuming the presence of the +effect +. +In our case, the increase in the posterior model probability of + +occurred because this model predicted the data slightly better after +removing the correlation between effect sizes and standard errors (a +consequence of using Fisher’s + +transformation). Nevertheless, the conditional effect size estimate +stayed almost the same.
+Now, we fit the publication bias-adjusted model by simply removing
+the priors_bias = NULL
argument, which allows us to obtain
+the default 36 models ensemble called RoBMA-PSMA (Bartoš et al., 2023).
+fit_RoBMA <- RoBMA(d = Poulsen2006$d, se = Poulsen2006$se, study_names = Poulsen2006$study,
+ priors_effect = prior_informed(name = "oral health", parameter = "effect", type = "smd"),
+ priors_heterogeneity = prior_informed(name = "oral health", parameter = "heterogeneity", type = "smd"),
+ seed = 1, parallel = TRUE)
+summary(fit_RoBMA, conditional = TRUE)
+#> Call:
+#> RoBMA(d = Poulsen2006$d, se = Poulsen2006$se, study_names = Poulsen2006$study,
+#> priors_effect = prior_informed(name = "oral health", parameter = "effect",
+#> type = "smd"), priors_heterogeneity = prior_informed(name = "oral health",
+#> parameter = "heterogeneity", type = "smd"), parallel = TRUE,
+#> seed = 1)
+#>
+#> Robust Bayesian meta-analysis
+#> Components summary:
+#> Models Prior prob. Post. prob. Inclusion BF
+#> Effect 18/36 0.500 0.858 6.022
+#> Heterogeneity 18/36 0.500 0.714 2.502
+#> Bias 32/36 0.500 0.697 2.304
+#>
+#> Model-averaged estimates:
+#> Mean Median 0.025 0.975
+#> mu 0.722 0.880 0.000 1.283
+#> tau 0.202 0.161 0.000 0.799
+#> omega[0,0.025] 1.000 1.000 1.000 1.000
+#> omega[0.025,0.05] 0.943 1.000 0.329 1.000
+#> omega[0.05,0.5] 0.874 1.000 0.071 1.000
+#> omega[0.5,0.95] 0.855 1.000 0.042 1.000
+#> omega[0.95,0.975] 0.866 1.000 0.050 1.000
+#> omega[0.975,1] 0.897 1.000 0.057 1.000
+#> PET 0.931 0.000 0.000 4.927
+#> PEESE 1.131 0.000 0.000 12.261
+#> The estimates are summarized on the Cohen's d scale (priors were specified on the Cohen's d scale).
+#> (Estimated publication weights omega correspond to one-sided p-values.)
+#>
+#> Conditional estimates:
+#> Mean Median 0.025 0.975
+#> mu 0.838 0.938 -0.035 1.297
+#> tau 0.285 0.227 0.064 0.906
+#> omega[0,0.025] 1.000 1.000 1.000 1.000
+#> omega[0.025,0.05] 0.736 0.829 0.092 1.000
+#> omega[0.05,0.5] 0.411 0.373 0.014 0.951
+#> omega[0.5,0.95] 0.320 0.249 0.008 0.919
+#> omega[0.95,0.975] 0.376 0.311 0.009 0.958
+#> omega[0.975,1] 0.518 0.425 0.010 1.000
+#> PET 2.909 3.136 0.171 5.614
+#> PEESE 7.048 6.034 0.375 18.162
+#> The estimates are summarized on the Cohen's d scale (priors were specified on the Cohen's d scale).
+#> (Estimated publication weights omega correspond to one-sided p-values.)
We notice the additional values in the ‘Components summary’ table in +the ‘Bias’ row. The model is now extended with 32 publication bias +adjustment models that account for 50% of the prior model probability. +When comparing the RoBMA to the second BMA fit, we notice a large +decrease in the inclusion Bayes factor for the presence of the effect + +vs. , +which still, however, presents moderate evidence for the presence of the +effect. We can further quantify the evidence in favor of the publication +bias with the inclusion Bayes factor for publication bias +, +which can be interpreted as weak evidence in favor of publication +bias.
+We can also compare the publication bias unadjusted and publication +bias-adjusted conditional effect size estimates. Including models +assuming publication bias into our model-averaged estimate (assuming the +presence of the effect) slightly decreases the estimated effect to +, +95% CI [-0.035, 1.297] with a much wider confidence interval, as +visualized in the prior and posterior conditional effect size estimate +plot.
+
+plot(fit_RoBMA, parameter = "mu", prior = TRUE, conditional = TRUE)
+The additional setting transformation = "cohens_d"
allows
+us to get more comparable results with the metaBMA
R
+package since RoBMA otherwise internally transforms the effect sizes to
+Fisher’s
+
+for the fitting purposes. The seed = 1
and
+parallel = TRUE
options grant us exact reproducibility of
+the results and parallelization of the fitting process.
+The model-averaged estimates that RoBMA
returns by default
+model-averaged across all specified models – a different behavior from
+the metaBMA
package that by default returns what we call
+“conditional” estimates in RoBMA
.
vignettes/MedicineBiBMA.Rmd
+ MedicineBiBMA.Rmd
Bayesian model-averaged meta-analysis can be specified using the +binomial likelihood and applied to data with dichotomous outcomes. This +vignette illustrates how to do this with an example from Bartoš et al. (2023), who implemented a +binomial-normal Bayesian model-averaged meta-analytic model and +developed informed prior distributions for meta-analyses of binary and +time-to-event outcomes based on the Cochrane database of systematic +reviews (see Bartoš et al. (2021) for +informed prior distributions for meta-analyses of continuous outcomes +highlighted in Informed Bayesian +Model-Averaged Meta-Analysis in Medicine vignette.
+We illustrate how to fit the binomial-normal Bayesian model-averaged
+meta-analysis using the RoBMA
R package. For this purpose,
+we reproduce the example of adverse effects of honey in treating acute
+cough in children from Bartoš et al.
+(2023), who reanalyzed two studies with adverse events of
+nervousness, insomnia, or hyperactivity in the honey vs. no treatment
+condition that were subjected to a meta-analysis by Oduwole et al. (2018).
We load the RoBMA package and specify the number of adverse events +and sample sizes in each arm as described on p. 73 (Oduwole et al., 2018).
+
+library(RoBMA)
+
+events_experimental <- c(5, 2)
+events_control <- c(0, 0)
+observations_experimental <- c(35, 40)
+observations_control <- c(39, 40)
+study_names <- c("Paul 2007", "Shadkam 2010")
Notice that both studies reported no adverse events in the control +group. Using a normal-normal meta-analytic model with log odds ratios +would require a continuity correction, which might result in bias. +Binomial-normal models allow us to circumvent the issue by modeling the +observed proportions directly (see Bartoš et al. +(2023) for more details).
+First, we fit the binomial-normal Bayesian model-averaged
+meta-analysis using informed prior distributions based on the
+Acute Respiratory Infections
subfield. We use the
+BiBMA
function and specify the observed events
+(x1
and x2
) and sample size (n1
+and n2
) of adverse events and sample sizes in each arm. We
+use the prior_informed
function to specify the informed
+prior distributions for the individual medical subfields
+automatically.
+fit <- BiBMA(
+ x1 = events_experimental,
+ x2 = events_control,
+ n1 = observations_experimental,
+ n2 = observations_control,
+ study_names = study_names,
+ priors_effect = prior_informed("Acute Respiratory Infections", type = "logOR", parameter = "effect"),
+ priors_heterogeneity = prior_informed("Acute Respiratory Infections", type = "logOR", parameter = "heterogeneity"),
+ seed = 1
+)
with priors_effect
and priors_heterogeneity
+corresponding to the
+
+and
+
+prior distributions (see ?prior_informed
for more details
+regarding the informed prior distributions).
We obtain the output with the summary
function. Adding
+the conditional = TRUE
argument allows us to inspect the
+conditional estimates, i.e., the effect size estimate assuming that the
+models specifying the presence of the effect are true, and the
+heterogeneity estimates assuming that the models specifying the presence
+of heterogeneity are true. We also set the
+output_scale = "OR"
argument to display the effect size
+estimates on the odds ratio scale.
+summary(fit, conditional = TRUE, output_scale = "OR")
+#> Call:
+#> BiBMA(x1 = events_experimental, x2 = events_control, n1 = observations_experimental,
+#> n2 = observations_control, study_names = study_names, priors_effect = prior_informed("Acute Respiratory Infections",
+#> type = "logOR", parameter = "effect"), priors_heterogeneity = prior_informed("Acute Respiratory Infections",
+#> type = "logOR", parameter = "heterogeneity"), seed = 1)
+#>
+#> Bayesian model-averaged meta-analysis (binomial-normal model)
+#> Components summary:
+#> Models Prior prob. Post. prob. Inclusion BF
+#> Effect 2/4 0.500 0.725 2.630
+#> Heterogeneity 2/4 0.500 0.564 1.296
+#>
+#> Model-averaged estimates:
+#> Mean Median 0.025 0.975
+#> mu 3.389 1.642 0.842 15.143
+#> tau 0.420 0.158 0.000 2.594
+#> The effect size estimates are summarized on the OR scale and heterogeneity is summarized on the logOR scale (priors were specified on the log(OR) scale).
+#>
+#> Conditional estimates:
+#> Mean Median 0.025 0.975
+#> mu 4.242 2.261 0.781 17.613
+#> tau 0.747 0.426 0.097 3.233
+#> The effect size estimates are summarized on the OR scale and heterogeneity is summarized on the logOR scale (priors were specified on the log(OR) scale).
The output from the summary.RoBMA()
function has three
+parts. The first part, under the ‘Robust Bayesian Meta-Analysis’ heading
+provides a basic summary of the fitted models by component types
+(presence of the Effect and Heterogeneity). The results show that the
+inclusion Bayes factor for the effect corresponds to the one reported in
+Bartoš et al. (2023),
+
+and
+
+(up to an MCMC error)—weak/undecided evidence for the presence of the
+effect and heterogeneity.
The second part, under the ‘Model-averaged estimates’ heading +displays the parameter estimates model-averaged across all specified +models (i.e., including models specifying the effect size to be zero). +These estimates are shrunk towards the null hypotheses of null effect or +no heterogeneity in accordance with the posterior uncertainty about the +presence of the effect or heterogeneity. We find the model-averaged mean +effect OR = 3.39, 95% CI [0.84, 15.14], and a heterogeneity estimate +, +95% CI [0.00, 2.59].
+The third part, under the ‘Conditional estimates’ heading displays +the conditional effect size and heterogeneity estimates (i.e., estimates +assuming presence of the effect or heterogeneity) corresponding to the +one reported in Bartoš et al. (2023), OR = +4.24, 95% CI [0.78, 17.61], and a heterogeneity estimate +, +95% CI [0.10, 3.23].
+We can also visualize the posterior distributions of the effect size
+and heterogeneity parameters using the plot()
function.
+Here, we set the conditional = TRUE
argument to display the
+conditional effect size estimate and prior = TRUE
to
+include the prior distribution in the plot.
+plot(fit, parameter = "mu", prior = TRUE, conditional = TRUE)
Additional visualizations and summaries are demonstrated in the Reproducing BMA and Informed Bayesian Model-Averaged Meta-Analysis +in Medicine vignettes.
+vignettes/MetaRegression.Rmd
+ MetaRegression.Rmd
Robust Bayesian model-averaged meta-regression (RoBMA-reg) extends
+the robust Bayesian model-averaged meta-analysis (RoBMA) by including
+covariates in the meta-analytic model. RoBMA-reg allows for estimating
+and testing the moderating effects of study-level covariates on the
+meta-analytic effect in a unified framework (e.g., accounting for
+uncertainty in the presence vs. absence of the effect, heterogeneity,
+and publication bias). This vignette illustrates how to fit a robust
+Bayesian model-averaged meta-regression using the RoBMA
R
+package. We reproduce the example from Bartoš,
+Maier, Stanley, et al. (2023), who re-analyzed a meta-analysis of
+the effect of household chaos on child executive functions with the mean
+age and assessment type covariates based on Andrews et al. (2021)’s meta-analysis.
First, we fit a frequentist meta-regression using the
+metafor
R package. Second, we explain the Bayesian
+meta-regression model specification, the default prior distributions for
+continuous and categorical moderators, and standardized effect sizes
+input specification. Third, we estimate Bayesian model-averaged
+meta-regression (without publication bias adjustment). Finally, we
+estimate the complete robust Bayesian model-averaged
+meta-regression.
We start by loading the Andrews2021
dataset included in
+the RoBMA
R package, which contains 36 estimates of the
+effect of household chaos on child executive functions with the mean age
+and assessment type covariates. The dataset includes correlation
+coefficients (r
), standard errors of the correlation
+coefficients (se
), the type of executive function
+assessment (measure
), and the mean age of the children
+(age
) in each study.
+library(RoBMA)
+data("Andrews2021", package = "RoBMA")
+head(Andrews2021)
+#> r se measure age
+#> 1 0.070 0.04743416 direct 4.606660
+#> 2 0.033 0.04371499 direct 2.480833
+#> 3 0.170 0.10583005 direct 7.750000
+#> 4 0.208 0.08661986 direct 4.000000
+#> 5 0.270 0.02641969 direct 4.000000
+#> 6 0.170 0.05147815 direct 4.487500
We start by fitting a frequentist meta-regression using the
+metafor
R package (Wolfgang,
+2010). While Andrews et al. (2021)
+estimated univariate meta-regressions for each moderator, we directly
+proceed by analyzing both moderators simultaneously. For consistency
+with original reporting, we estimate the meta-regression using the
+correlation coefficients and the standard errors provided by (Andrews et al., 2021); however, note that
+Fisher’s z transformation is recommended for estimating meta-analytic
+models (e.g., Stanley et al. (2024)).
+fit_rma <- metafor::rma(yi = r, sei = se, mods = ~ measure + age, data = Andrews2021)
+fit_rma
+#>
+#> Mixed-Effects Model (k = 36; tau^2 estimator: REML)
+#>
+#> tau^2 (estimated amount of residual heterogeneity): 0.0150 (SE = 0.0045)
+#> tau (square root of estimated tau^2 value): 0.1226
+#> I^2 (residual heterogeneity / unaccounted variability): 91.28%
+#> H^2 (unaccounted variability / sampling variability): 11.47
+#> R^2 (amount of heterogeneity accounted for): 15.24%
+#>
+#> Test for Residual Heterogeneity:
+#> QE(df = 33) = 340.7613, p-val < .0001
+#>
+#> Test of Moderators (coefficients 2:3):
+#> QM(df = 2) = 7.5445, p-val = 0.0230
+#>
+#> Model Results:
+#>
+#> estimate se zval pval ci.lb ci.ub
+#> intrcpt 0.0898 0.0467 1.9232 0.0545 -0.0017 0.1813 .
+#> measureinformant 0.1202 0.0466 2.5806 0.0099 0.0289 0.2115 **
+#> age 0.0030 0.0062 0.4867 0.6265 -0.0091 0.0151
+#>
+#> ---
+#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The results reveal a statistically significant moderation effect of
+the executive function assessment type on the effect of household chaos
+on child executive functions
+().
+To explore the moderation effect further, we estimate the estimated
+marginal means for the executive function assessment type using the
+emmeans
R package (Lenth et al.,
+2017).
+emmeans::emmeans(metafor::emmprep(fit_rma), specs = "measure")
+#> measure emmean SE df asymp.LCL asymp.UCL
+#> direct 0.109 0.0305 Inf 0.0492 0.169
+#> informant 0.229 0.0347 Inf 0.1612 0.297
+#>
+#> Confidence level used: 0.95
Studies using the informant-completed questionnaires show a stronger +effect of household chaos on child executive functions, r = 0.229, 95% +CI [0.161, 0.297], than the direct assessment, r = 0.109, 95% CI [0.049, +0.169]; both types of studies show statistically significant +effects.
+The mean age of the children does not significantly moderate the +effect +() +with the estimated regression coefficient of b = 0.003, 95% CI [-0.009, +0.015]. As usual, frequentist inference limits us to failing to reject +the null hypothesis. Here, we try to overcome this limitation with +Bayesian model-averaged meta-regression.
+Before we proceed with the Bayesian model-averaged meta-regression, +we provide a quick overview of the regression model specification. In +contrast to frequentist meta-regression, we need to specify prior +distributions on the regression coefficients, which encode the tested +hypotheses about the presence vs. absence of the moderation (specifying +different prior distributions corresponds to different hypotheses and +results in different conclusions). Importantly, the treatment of +continuous and categorical covariates differs in the Bayesian +model-averaged meta-regression.
+The default prior distribution for continuous moderators is a normal
+prior distribution with mean of 0 and a standard deviation of 1/4. In
+other words, the default prior distribution assumes that the effect of
+the moderator is small and smaller moderation effects are more likely
+than larger effects. The default choice for continuous moderators can be
+overridden by the prior_covariates
argument (for all
+continuous covariates) or by the priors
argument (for
+specific covariates, see ?RoBMA.reg
for more information).
+The package automatically standardizes the continuous moderators. This
+achieves scale-invariance of the specified prior distributions and
+ensures that the prior distribution for the intercept correspond to the
+grand mean effect. This setting can be overridden by specifying the
+standardize_predictors = FALSE
argument.
The default prior distribution for the categorical moderators is a
+normal distribution with a mean of 0 and a standard deviation of 1/4,
+representing the deviation of each level from the grand mean effect. The
+package uses standardized orthonormal contrasts
+(contrast = "meandif"
) to model deviations of each category
+from the grand mean effect. The default choice for categorical
+moderators can be overridden by the prior_factors
argument
+(for all categorical covariates) or by the priors
argument
+(for specific covariates, see ?RoBMA.reg
for more
+information). The "meandif"
contrasts achieve label
+invariance (i.e., the coding of the categorical covariates does not
+affect the results) and the prior distribution for the intercept
+corresponds to the grand mean effect. Alternatively, the package also
+allows specifying "treatment"
contrasts, which result in a
+prior distribution on the difference between the default level and the
+remaining levels of the categorical covariate (with the intercept
+corresponding to the effect in the default factor level).
Prior distributions for Bayesian meta-analyses are calibrated for the
+standardized effect size measures. As such, the fitting function needs
+to know what kind of effect size was supplied as the input. In
+RoBMA()
function, this is achieved by the d
,
+r
, logOR
, OR
, z
,
+se
, v
, n
, lCI
, and
+uCI
arguments. The input is passed to the
+combine_data()
function in the background that combines the
+effect sizes and merges them into a single data.frame. The
+RoBMA.reg()
(and NoBMA.reg()
) function
+requires the dataset to be passed as a data.frame (without missing
+values) with column names identifying the - moderators passed using the
+formula interface (i.e., ~ measure + age
in our example) -
+and the effect sizes and standard errors (i.e., r
and
+se
in our example).
As such, it is crucial for the column names to correctly identify the +standardized effect sizes, standard errors, sample sizes, and +moderators.
+We fit the Bayesian model-averaged meta-regression using the
+NoBMA.reg()
function (the NoBMA.reg()
function
+is a wrapper around the RoBMA.reg()
function that
+automatically removes models adjusting for publication bias). We specify
+the model formula with the ~
operator similarly to the
+rma()
function and pass the dataset as a data.frame with
+named columns as outlined in the section above (the names need to
+identify the moderators and effect size measures). We also set the
+parallel = TRUE
argument to speed up the computation by
+running the chains in parallel and seed = 1
argument to
+ensure reproducibility.
+fit_BMA <- NoBMA.reg(~ measure + age, data = Andrews2021, parallel = TRUE, seed = 1)
Note that the NoBMA.reg()
function specifies the
+combination of all models assuming presence vs. absence of the effect,
+heterogeneity, moderation by measure
, and moderation by
+age
, which corresponds to
+
+models. Including each additional moderator doubles the number of
+models, leading to an exponential increase in model count and
+significantly longer fitting times.
Once the ensemble is estimated, we can use the summary()
+functions with the output_scale = "r"
argument, which
+produces meta-analytic estimates that are transformed to the correlation
+scale.
+summary(fit_BMA, output_scale = "r")
+#> Call:
+#> RoBMA.reg(formula = formula, data = data, test_predictors = test_predictors,
+#> study_names = study_names, study_ids = study_ids, transformation = transformation,
+#> prior_scale = prior_scale, standardize_predictors = standardize_predictors,
+#> effect_direction = "positive", priors = priors, model_type = model_type,
+#> priors_effect = priors_effect, priors_heterogeneity = priors_heterogeneity,
+#> priors_bias = NULL, priors_effect_null = priors_effect_null,
+#> priors_heterogeneity_null = priors_heterogeneity_null, priors_bias_null = prior_none(),
+#> priors_hierarchical = priors_hierarchical, priors_hierarchical_null = priors_hierarchical_null,
+#> prior_covariates = prior_covariates, prior_covariates_null = prior_covariates_null,
+#> prior_factors = prior_factors, prior_factors_null = prior_factors_null,
+#> chains = chains, sample = sample, burnin = burnin, adapt = adapt,
+#> thin = thin, parallel = parallel, autofit = autofit, autofit_control = autofit_control,
+#> convergence_checks = convergence_checks, save = save, seed = seed,
+#> silent = silent)
+#>
+#> Bayesian model-averaged meta-regression (normal-normal model)
+#> Components summary:
+#> Models Prior prob. Post. prob. Inclusion BF
+#> Effect 8/16 0.500 1.000 6.637645e+05
+#> Heterogeneity 8/16 0.500 1.000 3.439130e+40
+#>
+#> Meta-regression components summary:
+#> Models Prior prob. Post. prob. Inclusion BF
+#> measure 8/16 0.500 0.826 4.739
+#> age 8/16 0.500 0.197 0.245
+#>
+#> Model-averaged estimates:
+#> Mean Median 0.025 0.975
+#> mu 0.163 0.163 0.118 0.208
+#> tau 0.121 0.120 0.086 0.167
+#> The effect size estimates are summarized on the correlation scale and heterogeneity is summarized on the Fisher's z scale (priors were specified on the Cohen's d scale).
+#>
+#> Model-averaged meta-regression estimates:
+#> Mean Median 0.025 0.975
+#> intercept 0.163 0.163 0.118 0.208
+#> measure [dif: direct] -0.047 -0.051 -0.099 0.000
+#> measure [dif: informant] 0.047 0.051 0.000 0.099
+#> age 0.003 0.000 -0.011 0.043
+#> The effect size estimates are summarized on the correlation scale and heterogeneity is summarized on the Fisher's z scale (priors were specified on the Cohen's d scale).
The summary function produces output with multiple sections The first
+section contains the Components summary
with the hypothesis
+test results for the overall effect size and heterogeneity. We find
+overwhelming evidence for both with inclusion Bayes factors
+(Inclusion BF
) above 10,000.
The second section contains the
+Meta-regression components summary
with the hypothesis test
+results for the moderators. We find moderate evidence for the moderation
+by the executive function assessment type,
+.
+Furthermore, we find moderate evidence for the null hypothesis of no
+moderation by mean age of the children,
+
+(i.e., BF for the null is
+).
+These findings extend the frequentist meta-regression by disentangling
+the absence of evidence from the evidence of absence.
The third section contains the Model-averaged estimates
+with the model-averaged estimates for mean effect
+,
+95% CI [0.12, 0.21] and between-study heterogeneity
+,
+95% CI [0.09, 0.17].
The fourth section contains the
+Model-averaged meta-regression estimates
with the
+model-averaged regression coefficient estimates. The main difference
+from the usual frequentist meta-regression output is that the
+categorical predictors are summarized as a difference from the grand
+mean for each factor level. Here, the intercept
regression
+coefficient estimate corresponds to the grand mean effect and the
+measure [dif: direct]
regression coefficient estimate of
+-0.047, 95% CI [-0.099, 0.000] corresponds to the difference between the
+direct assessment and the grand mean. As such, the results suggest that
+the effect size in studies using direct assessment is lower in
+comparison to the grand mean of the studies. The age
+regression coefficient estimate is standardized, therefore, the increase
+of 0.003, 95% CI [-0.011, 0.043] corresponds to the increase in the mean
+effect when increasing mean age of children by one standard
+deviation.
Similarly to the frequentist meta-regression, we can use the
+marginal_summary()
function to obtain the marginal
+estimates for each of the factor levels.
+marginal_summary(fit_BMA, output_scale = "r")
+#> Call:
+#> RoBMA.reg(formula = formula, data = data, test_predictors = test_predictors,
+#> study_names = study_names, study_ids = study_ids, transformation = transformation,
+#> prior_scale = prior_scale, standardize_predictors = standardize_predictors,
+#> effect_direction = "positive", priors = priors, model_type = model_type,
+#> priors_effect = priors_effect, priors_heterogeneity = priors_heterogeneity,
+#> priors_bias = NULL, priors_effect_null = priors_effect_null,
+#> priors_heterogeneity_null = priors_heterogeneity_null, priors_bias_null = prior_none(),
+#> priors_hierarchical = priors_hierarchical, priors_hierarchical_null = priors_hierarchical_null,
+#> prior_covariates = prior_covariates, prior_covariates_null = prior_covariates_null,
+#> prior_factors = prior_factors, prior_factors_null = prior_factors_null,
+#> chains = chains, sample = sample, burnin = burnin, adapt = adapt,
+#> thin = thin, parallel = parallel, autofit = autofit, autofit_control = autofit_control,
+#> convergence_checks = convergence_checks, save = save, seed = seed,
+#> silent = silent)
+#>
+#> Robust Bayesian meta-analysis
+#> Model-averaged marginal estimates:
+#> Mean Median 0.025 0.975 Inclusion BF
+#> intercept 0.163 0.163 0.118 0.208 Inf
+#> measure[direct] 0.117 0.116 0.052 0.185 50.151
+#> measure[informant] 0.208 0.210 0.130 0.280 Inf
+#> age[-1SD] 0.160 0.161 0.106 0.208 Inf
+#> age[0SD] 0.163 0.163 0.118 0.208 Inf
+#> age[1SD] 0.166 0.165 0.117 0.220 Inf
+#> The estimates are summarized on the correlation scale (priors were specified on the Cohen's d scale).
+#> mu_intercept: Posterior samples do not span both sides of the null hypothesis. The Savage-Dickey density ratio is likely to be overestimated.
+#> mu_measure[informant]: Posterior samples do not span both sides of the null hypothesis. The Savage-Dickey density ratio is likely to be overestimated.
+#> mu_age[-1SD]: Posterior samples do not span both sides of the null hypothesis. The Savage-Dickey density ratio is likely to be overestimated.
+#> mu_age[0SD]: There is a considerable cluster of prior samples at the exact null hypothesis values. The Savage-Dickey density ratio is likely to be invalid.
+#> mu_age[0SD]: Posterior samples do not span both sides of the null hypothesis. The Savage-Dickey density ratio is likely to be overestimated.
+#> mu_age[1SD]: Posterior samples do not span both sides of the null hypothesis. The Savage-Dickey density ratio is likely to be overestimated.
The estimated marginal means are similar to the frequentist results. +Studies using the informant-completed questionnaires again show a +stronger effect of household chaos on child executive functions, +, +95% CI [0.130, 0.280], than the direct assessment, +, +95% CI [0.052, 0.185].
+The last column summarizes results from a test against a null +hypothesis of marginal means equals 0. Here, we find very strong +evidence for the effect size of studies using the informant-completed +questionnaires differing from zero, + +and extreme evidence for the effect size of studies using the direct +assessment differing from zero, +. +The test is performed using the change from prior to posterior +distribution at 0 (i.e., the Savage-Dickey density ratio) assuming the +presence of the overall effect or the presence of difference according +to the tested factor. Because the tests use prior and posterior samples, +calculating the Bayes factor can be problematic when the posterior +distribution is far from the tested value. In such cases, warning +messages are printed and + +returned (like here)—while the actual Bayes factor is less than +infinity, it is still too large to be computed precisely given the +posterior samples.
+The full model-averaged posterior marginal means distribution can be
+visualized by the marginal_plot()
function.
+marginal_plot(fit_BMA, parameter = "measure", output_scale = "r", lwd = 2)
Finally, we adjust the Bayesian model-averaged meta-regression model
+by fitting the robust Bayesian model-averaged meta-regression. In
+contrast to the previous publication bias unadjusted model ensemble,
+RoBMA-reg extends the model ensemble by the publication bias component
+specified via 6 weight functions and PET-PEESE (Bartoš, Maier, Wagenmakers, et al., 2023). We
+use the RoBMA.reg()
function with the same arguments as in
+the previous section. The estimation time further increases as the
+ensemble now contains 144 models.
+fit_RoBMA <- RoBMA.reg(~ measure + age, data = Andrews2021, parallel = TRUE, seed = 1)
+summary(fit_RoBMA, output_scale = "r")
+#> Call:
+#> RoBMA.reg(formula = ~measure + age, data = Andrews2021, chains = 1,
+#> parallel = TRUE, seed = 1)
+#>
+#> Robust Bayesian meta-regression
+#> Components summary:
+#> Models Prior prob. Post. prob. Inclusion BF
+#> Effect 72/144 0.500 0.334 5.020000e-01
+#> Heterogeneity 72/144 0.500 1.000 1.043816e+23
+#> Bias 128/144 0.500 0.965 2.795800e+01
+#>
+#> Meta-regression components summary:
+#> Models Prior prob. Post. prob. Inclusion BF
+#> measure 72/144 0.500 0.950 19.086
+#> age 72/144 0.500 0.154 0.182
+#>
+#> Model-averaged estimates:
+#> Mean Median 0.025 0.975
+#> mu 0.031 0.000 0.000 0.164
+#> tau 0.106 0.104 0.074 0.147
+#> omega[0,0.025] 1.000 1.000 1.000 1.000
+#> omega[0.025,0.05] 0.999 1.000 1.000 1.000
+#> omega[0.05,0.5] 0.998 1.000 1.000 1.000
+#> omega[0.5,0.95] 0.997 1.000 1.000 1.000
+#> omega[0.95,0.975] 0.997 1.000 1.000 1.000
+#> omega[0.975,1] 0.997 1.000 1.000 1.000
+#> PET 2.056 2.494 0.000 3.293
+#> PEESE 1.916 0.000 0.000 19.068
+#> The effect size estimates are summarized on the correlation scale and heterogeneity is summarized on the Fisher's z scale (priors were specified on the Cohen's d scale).
+#> (Estimated publication weights omega correspond to one-sided p-values.)
+#>
+#> Model-averaged meta-regression estimates:
+#> Mean Median 0.025 0.975
+#> intercept 0.031 0.000 0.000 0.164
+#> measure [dif: direct] -0.063 -0.064 -0.106 0.000
+#> measure [dif: informant] 0.063 0.064 0.000 0.106
+#> age 0.000 0.000 -0.024 0.022
+#> The effect size estimates are summarized on the correlation scale and heterogeneity is summarized on the Fisher's z scale (priors were specified on the Cohen's d scale).
All previously described functions for manipulating the fitted model +work identically with the publication bias adjusted model. As such, we +just briefly mention the main differences found after adjusting for +publication bias.
+RoBMA-reg reveals strong evidence of publication bias +. +Furthermore, accounting for publication bias turns the previously found +evidence for the overall effect into a weak evidence against the effect + +and notably reduces the mean effect estimate +, +95% CI [0.000, 0.164].
+
+marginal_summary(fit_RoBMA, output_scale = "r")
+#> Call:
+#> RoBMA.reg(formula = ~measure + age, data = Andrews2021, chains = 1,
+#> parallel = TRUE, seed = 1)
+#>
+#> Robust Bayesian meta-analysis
+#> Model-averaged marginal estimates:
+#> Mean Median 0.025 0.975 Inclusion BF
+#> intercept 0.031 0.000 0.000 0.164 0.516
+#> measure[direct] -0.031 -0.056 -0.105 0.121 0.575
+#> measure[informant] 0.093 0.077 0.000 0.223 7.643
+#> age[-1SD] 0.031 0.000 -0.015 0.163 0.732
+#> age[0SD] 0.031 0.000 0.000 0.164 1.013
+#> age[1SD] 0.031 0.000 -0.024 0.168 0.743
+#> The estimates are summarized on the correlation scale (priors were specified on the Cohen's d scale).
+#> mu_age[0SD]: There is a considerable cluster of prior samples at the exact null hypothesis values. The Savage-Dickey density ratio is likely to be invalid.
The estimated marginal means now suggest that studies using the +informant-completed questionnaires show a much smaller effect of +household chaos on child executive functions, +, +95% CI [0.000, 0.223] with only moderate evidence against no effect, +, +while studies using direct assessment even provide weak evidence against +the effect of household chaos on child executive functions, +, +with most likely effect sizes around zero, +, +95% CI [-0.105, 0.121].
+A visual summary of the estimated marginal means highlights the +considerably wider model-averaged posterior distributions of the +marginal means—a consequence of accounting and adjusting for publication +bias.
+
+marginal_plot(fit_RoBMA, parameter = "measure", output_scale = "r", lwd = 2)
The Bayesian model-averaged meta-regression models are compatible
+with the remaining custom specification, visualization, and summary
+functions included in the RoBMA
R package, highlighted in
+other vignettes. E.g., custom model specification is demonstrated in the
+vignette Fitting Custom Meta-Analytic
+Ensembles and visualizations and summaries are demonstrated in the
+Reproducing BMA and Informed Bayesian Model-Averaged Meta-Analysis
+in Medicine vignettes.
vignettes/ReproducingBMA.Rmd
+ ReproducingBMA.Rmd
By default, the RoBMA package estimates an ensemble of 36 +meta-analytic models and provides functions for convenient manipulation +of the fitted object. However, the package has been designed so it can +be used as a framework for estimating any combination of meta-analytic +models (or a single model). Here, we illustrate how to build a custom +ensemble of meta-analytic models - specifically the same ensemble that +is used in ‘classical’ Bayesian Model-Averaged Meta-Analysis (Bartoš et al., 2021; Gronau et al., 2017, +2021). See this vignette if +you are interested in building more customized ensembles or Bartoš et al. (2022) for a tutorial on fitting +(custom) models in JASP.
+We illustrate how to fit a classical BMA (not adjusting for
+publication bias) using RoBMA
. For this purpose, we
+reproduce a meta-analysis of registered reports on Power posing by Gronau et al. (2017). We focus only on the
+analysis of all reported results using a Cauchy prior distribution with
+scale
+
+for the effect size estimation (half-Cauchy for testing) and
+inverse-gamma distribution with shape = 1 and scale = 0.15 for the
+heterogeneity parameter. You can find the figure from the original
+publication here
+and the paper’s supplementary materials at https://osf.io/fxg32/.
First, we load the power posing data provided within the metaBMA +package and reproduce the analysis performed by Gronau et al. (2017).
+
+data("power_pose", package = "metaBMA")
+power_pose[,c("study", "effectSize", "SE")]
+#> study effectSize SE
+#> 1 Bailey et al. 0.2507640 0.2071399
+#> 2 Ronay et al. 0.2275180 0.1931046
+#> 3 Klaschinski et al. 0.3186069 0.1423228
+#> 4 Bombari et al. 0.2832082 0.1421356
+#> 5 Latu et al. 0.1463949 0.1416107
+#> 6 Keller et al. 0.1509773 0.1221166
+fit_BMA_test <- metaBMA::meta_bma(y = power_pose$effectSize, SE = power_pose$SE,
+ d = metaBMA::prior(family = "halfcauchy", param = 1/sqrt(2)),
+ tau = metaBMA::prior(family = "invgamma", param = c(1, .15)))
+
+fit_BMA_est <- metaBMA::meta_bma(y = power_pose$effectSize, SE = power_pose$SE,
+ d = metaBMA::prior(family = "cauchy", param = c(0, 1/sqrt(2))),
+ tau = metaBMA::prior(family = "invgamma", param = c(1, .15)))
+fit_BMA_test$inclusion
+#> ### Inclusion Bayes factor ###
+#> Model Prior Posterior included
+#> 1 fixed_H0 0.25 0.00868
+#> 2 fixed_H1 0.25 0.77745 x
+#> 3 random_H0 0.25 0.02061
+#> 4 random_H1 0.25 0.19325 x
+#>
+#> Inclusion posterior probability: 0.971
+#> Inclusion Bayes factor: 33.136
+
+round(fit_BMA_est$estimates,2)
+#> mean sd 2.5% 50% 97.5% hpd95_lower hpd95_upper n_eff Rhat
+#> averaged 0.22 0.06 0.09 0.22 0.34 0.09 0.34 NA NA
+#> fixed 0.22 0.06 0.10 0.22 0.34 0.10 0.34 3026.5 1
+#> random 0.22 0.08 0.07 0.22 0.37 0.07 0.37 6600.4 1
From the output, we can see the inclusion Bayes factor for the effect
+size was
+
+and the effect size estimate 0.22, 95% HDI [0.09, 0.34], which matches
+the reported results. Please note that the metaBMA
package
+model-averages only across the
+
+models, whereas the RoBMA
package model-averages across all
+models (assuming the presence and absence of the effect).
Now we reproduce the analysis with RoBMA
. We set the
+corresponding prior distributions for effect sizes
+()
+and heterogeneity
+(),
+and remove the alternative prior distributions for the publication bias
+by setting priors_bias = NULL
. To specify the half-Cauchy
+prior distribution with the RoBMA::prior()
function we use
+a regular Cauchy distribution and truncate it at zero (note that both
+metaBMA
and RoBMA
export their own
+prior()
functions that will clash when loading both
+packages simultaneously). The inverse-gamma prior distribution for the
+heterogeneity parameter is the default option (we specify it for
+completeness). We omit the specifications for the null prior
+distributions for the effect size, heterogeneity (both of which are set
+to a spike at 0 by default), and publication bias (which is set to no
+publication bias by default). Note that starting from version 3.1, the
+package includes the NoBMA()
function, which allows users
+to skip publication bias adjustment directly.
Since metaBMA
model-averages the effect size estimates
+only across the models assuming presence of the effect, we remove the
+models assuming absence of the effect from the estimation ensemble with
+priors_effect_null = NULL
. Finally, we set
+transformation = "cohens_d"
to estimate the models on
+Cohen’s d scale. RoBMA uses Fisher’s z scale by
+default and transforms the estimated coefficients back to the scale that
+is used for specifying the prior distributions. We speed up the
+computation by setting parallel = TRUE
, and set a seed for
+reproducibility.
+library(RoBMA)
+
+fit_RoBMA_test <- RoBMA(d = power_pose$effectSize, se = power_pose$SE, study_names = power_pose$study,
+ priors_effect = prior(
+ distribution = "cauchy",
+ parameters = list(location = 0, scale = 1/sqrt(2)),
+ truncation = list(0, Inf)),
+ priors_heterogeneity = prior(
+ distribution = "invgamma",
+ parameters = list(shape = 1, scale = 0.15)),
+ priors_bias = NULL,
+ transformation = "cohens_d", seed = 1, parallel = TRUE)
+
+fit_RoBMA_est <- RoBMA(d = power_pose$effectSize, se = power_pose$SE, study_names = power_pose$study,
+ priors_effect = prior(
+ distribution = "cauchy",
+ parameters = list(location = 0, scale = 1/sqrt(2))),
+ priors_heterogeneity = prior(
+ distribution = "invgamma",
+ parameters = list(shape = 1, scale = 0.15)),
+ priors_bias = NULL,
+ priors_effect_null = NULL,
+ transformation = "cohens_d", seed = 2, parallel = TRUE)
+summary(fit_RoBMA_test)
+#> Call:
+#> RoBMA(d = power_pose$effectSize, se = power_pose$SE, study_names = power_pose$study,
+#> transformation = "cohens_d", priors_effect = prior(distribution = "cauchy",
+#> parameters = list(location = 0, scale = 1/sqrt(2)), truncation = list(0,
+#> Inf)), priors_heterogeneity = prior(distribution = "invgamma",
+#> parameters = list(shape = 1, scale = 0.15)), priors_bias = NULL,
+#> parallel = TRUE, seed = 1)
+#>
+#> Robust Bayesian meta-analysis
+#> Components summary:
+#> Models Prior prob. Post. prob. Inclusion BF
+#> Effect 2/4 0.500 0.971 33.112
+#> Heterogeneity 2/4 0.500 0.214 0.273
+#>
+#> Model-averaged estimates:
+#> Mean Median 0.025 0.975
+#> mu 0.213 0.217 0.000 0.348
+#> tau 0.022 0.000 0.000 0.178
+#> The estimates are summarized on the Cohen's d scale (priors were specified on the Cohen's d scale).
+
+summary(fit_RoBMA_est)
+#> Call:
+#> RoBMA(d = power_pose$effectSize, se = power_pose$SE, study_names = power_pose$study,
+#> transformation = "cohens_d", priors_effect = prior(distribution = "cauchy",
+#> parameters = list(location = 0, scale = 1/sqrt(2))),
+#> priors_heterogeneity = prior(distribution = "invgamma", parameters = list(shape = 1,
+#> scale = 0.15)), priors_bias = NULL, priors_effect_null = NULL,
+#> parallel = TRUE, seed = 2)
+#>
+#> Robust Bayesian meta-analysis
+#> Components summary:
+#> Models Prior prob. Post. prob. Inclusion BF
+#> Effect 2/2 1.000 1.000 Inf
+#> Heterogeneity 1/2 0.500 0.200 0.250
+#>
+#> Model-averaged estimates:
+#> Mean Median 0.025 0.975
+#> mu 0.220 0.220 0.096 0.346
+#> tau 0.019 0.000 0.000 0.152
+#> The estimates are summarized on the Cohen's d scale (priors were specified on the Cohen's d scale).
The output from the summary.RoBMA()
function has 2
+parts. The first one under the “Robust Bayesian Meta-Analysis” heading
+provides a basic summary of the fitted models by component types
+(presence of the Effect and Heterogeneity). The table summarizes the
+prior and posterior probabilities and the inclusion Bayes factors of the
+individual components. The results for the half-Cauchy model specified
+for testing show that the inclusion BF is nearly identical to the one
+computed by the metaBMA
package,
+.
The second part under the ‘Model-averaged estimates’ heading displays
+the parameter estimates. The results for the unrestricted Cauchy model
+specified for estimation show the effect size estimate
+,
+95% CI [0.10, 0.35] that also mirrors the one obtained from
+metaBMA
package.
RoBMA provides extensive options for visualizing the results. Here, +we visualize the prior (grey) and posterior (black) distribution for the +mean parameter.
+ + +If we visualize the effect size from the model specified for testing,
+we notice a few more things. The function plots the model-averaged
+estimates across all models by default, including models assuming the
+absence of the effect. The arrows represents the probability of a spike,
+here, at the value 0. The secondary y-axis (right) shows the probability
+of the value 0 decreased from 0.50, to 0.03 (also obtainable from the
+“Robust Bayesian Meta-Analysis” field in the
+summary.RoBMA()
function). Furthermore, the continuous
+prior distributions for the effect size under the alternative hypothesis
+are truncated to only positive values, reflecting the assumption that
+the effect size cannot be negative.
We can also visualize the estimates from the individual models used
+in the ensemble. We do that with the plot_models()
+function, which visualizes the effect size estimates and 95% CI of each
+of the specified models from the estimation ensemble (Model 1
+corresponds to the fixed effect model and Model 2 to the random effect
+model). The size of the square representing the mean estimate reflects
+the posterior model probability of the model, which is also displayed in
+the right-hand side panel. The bottom part of the figure shows the
+model-averaged estimate that is a combination of the individual model
+posterior distributions weighted by the posterior model
+probabilities.
+plot_models(fit_RoBMA_est)
The last type of visualization that we show here is the forest plot.
+It displays both the effect sizes from the original studies and the
+overall meta-analytic estimate in a single figure. It can be requested
+by using the forest()
function.
+forest(fit_RoBMA_est)
For more options provided by the plotting function, see its
+documentation using ?plot.RoBMA()
,
+?plot_models()
, and ?forest()
.
vignettes/Tutorial.Rmd
+ Tutorial.Rmd
This R markdown file accompanies the tutorial Adjusting for +publication bias in JASP and R: Selection models, PET-PEESE, and robust +Bayesian meta-analysis published in Advances in Methods and +Practices in Psychological Science (Bartoš +et al., 2022).
+The following R-markdown file illustrates how to:
+See the full paper for additional details regarding the data set, +methods, and interpretation.
+Before we start, we need to install JAGS
(which is
+needed for installation of the RoBMA
package) and the R
+packages that we use in the analysis. Specifically the
+RoBMA
, weightr
, and metafor
R
+packages.
JAGS can be downloaded from the JAGS website.
+Subsequently, we install the R packages with the
+install.packages()
function.
{r install.packages(c("RoBMA", "weightr", "metafor"))
If
+you happen to use the new M1 Mac machines with Apple silicon, see this
+blogpost outlining how to install JAGS on M1. In short, you will
+have to install Intel version of R (Intel/x86-64) from CRAN, not the Arm64
+(Apple silicon) version. Note that there might have been some changes in
+the installation process since the blogpost was written and there might
+be a JAGS version compatible with Apple silicon available now.
Once all of the packages are installed, we can load them into the
+workspace with the library()
function.
Lui (2015) studied how the +acculturation mismatch (AM) that is the result of the contrast between +the collectivist cultures of Asian and Latin immigrant groups and the +individualist culture in the United States correlates with +intergenerational cultural conflict (ICC). Lui +(2015) meta-analyzed 18 independent studies correlating AM with +ICC. A standard reanalysis indicates a significant effect of AM on +increased ICC, r = 0.250, p < .001.
+First, we load the Lui2015.csv file into R with the
+read.csv()
function and inspect the first six data entries
+with the head()
function (the data set is also included in
+the package and can be accessed via the
+data("Lui2015", package = "RoBMA")
call).
+df <- read.csv(file = "Lui2015.csv")
+head(df)
+#> r n study
+#> 1 0.21 115 Ahn, Kim, & Park (2008)
+#> 2 0.29 283 Basanez et al. (2013)
+#> 3 0.22 80 Bounkeua (2007)
+#> 4 0.26 109 Hajizadeh (2009)
+#> 5 0.23 61 Hamid (2007)
+#> 6 0.54 107 Hwang & Wood (2009a)
We see that the data set contains three columns. The first column
+called r
contains the effect sizes coded as correlation
+coefficients, the second column called n
contains the
+sample sizes, and the third column called study
contains
+names of the individual studies.
We can access the individual variables using the data set name and
+the dollar ($
) sign followed by the name of the column. For
+example, we can print all of the effect sizes with the df$r
+command.
+df$r
+#> [1] 0.21 0.29 0.22 0.26 0.23 0.54 0.56 0.29 0.26 0.02 -0.06 0.38
+#> [13] 0.25 0.08 0.17 0.33 0.36 0.13
The printed output shows that the data set contains mostly positive +effect sizes with the largest correlation coefficient r = 0.54.
+Before we start analyzing the data, we transform the effect sizes +from correlation coefficients + +to Fisher’s z. Correlation coefficients are not well suited for +meta-analysis because (1) they are bounded to a range (-1, 1) with +non-linear increases near the boundaries and (2) the standard error of +the correlation coefficients is related to the effect size. Fisher’s +z transformation mitigates both issues. It unwinds the (-1, 1) +range to +(, +), +makes the sampling distribution approximately normal, and breaks the +dependency between standard errors and effect sizes.
+To apply the transformation, we use the combine_data()
+function from the RoBMA
package. We pass the correlation
+coefficients into the r
argument, the sample sizes to the
+n
argument, and set the transformation
+argument to "fishers_z"
(the study_names
+argument is optional). The function combine_data()
then
+saves the transformed effect size estimates into a data frame called
+dfz
, where the y
column corresponds to
+Fisher’s z transformation of the correlation coefficient and
+se
column corresponds to the standard error of Fisher’s
+z.
+dfz <- combine_data(r = df$r, n = df$n, study_names = df$study, transformation = "fishers_z")
+head(dfz)
+#> y se study_names study_ids weight
+#> 1 0.2131713 0.09449112 Ahn, Kim, & Park (2008) NA NA
+#> 2 0.2985663 0.05976143 Basanez et al. (2013) NA NA
+#> 3 0.2236561 0.11396058 Bounkeua (2007) NA NA
+#> 4 0.2661084 0.09712859 Hajizadeh (2009) NA NA
+#> 5 0.2341895 0.13130643 Hamid (2007) NA NA
+#> 6 0.6041556 0.09805807 Hwang & Wood (2009a) NA NA
We can also transform the effect sizes according to Cohen’s +d transformation (which we utilize later to fit the selection +models).
+
+dfd <- combine_data(r = df$r, n = df$n, study_names = df$study, transformation = "cohens_d")
+head(dfd)
+#> y se study_names study_ids weight
+#> 1 0.4295790 0.1886397 Ahn, Kim, & Park (2008) NA NA
+#> 2 0.6060437 0.1215862 Basanez et al. (2013) NA NA
+#> 3 0.4510508 0.2264322 Bounkeua (2007) NA NA
+#> 4 0.5385205 0.1950065 Hajizadeh (2009) NA NA
+#> 5 0.4726720 0.2596249 Hamid (2007) NA NA
+#> 6 1.2831708 0.2123140 Hwang & Wood (2009a) NA NA
We now estimate a random effect meta-analysis with the
+rma()
function imported from the metafor
+package (Wolfgang, 2010) and verify that
+we arrive at the same results as reported in the Lui (2015) paper. The yi
argument
+is used to pass the column name containing effect sizes, the
+sei
argument is used to pass the column name containing
+standard errors, and the data
argument is used to pass the
+data frame containing both variables.
+fit_rma <- rma(yi = y, sei = se, data = dfz)
+fit_rma
+#>
+#> Random-Effects Model (k = 18; tau^2 estimator: REML)
+#>
+#> tau^2 (estimated amount of total heterogeneity): 0.0229 (SE = 0.0107)
+#> tau (square root of estimated tau^2 value): 0.1513
+#> I^2 (total heterogeneity / total variability): 77.79%
+#> H^2 (total variability / sampling variability): 4.50
+#>
+#> Test for Heterogeneity:
+#> Q(df = 17) = 73.5786, p-val < .0001
+#>
+#> Model Results:
+#>
+#> estimate se zval pval ci.lb ci.ub
+#> 0.2538 0.0419 6.0568 <.0001 0.1717 0.3359 ***
+#>
+#> ---
+#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Indeed, we find that the effect size estimate from the random effect
+meta-analysis corresponds to the one reported in Lui (2015). It is important to remember that we
+used Fisher’s z to estimate the models; therefore, the
+estimated results are on the Fisher’s z scale. To transform the
+effect size estimate to the correlation coefficients, we can use the
+z2r()
function from the RoBMA
package,
+z2r(fit_rma$b)
+#> [,1]
+#> intrcpt 0.2484877
Transforming the effect size estimate results in the correlation +coefficient + += 0.25.
+The first publication bias adjustment that we perform is PET-PEESE.
+PET-PEESE adjusts for the relationship between effect sizes and standard
+errors. To our knowledge, PET-PEESE is not currently implemented in any
+R-package. However, since PET and PEESE are weighted regressions of
+effect sizes on standard errors (PET) or standard errors squared
+(PEESE), we can estimate both PET and PEESE models with the
+lm()
function. Inside the lm()
function call,
+we specify that y
is the response variable (left hand side
+of the ~
sign) and se
is the predictor (the
+right-hand side). Furthermore, we specify the weights
+argument that allows us to weight the meta-regression by inverse
+variance and set the data = dfz
argument, which specifies
+that all of the variables come from the transformed, dfz
,
+data set.
+fit_PET <- lm(y ~ se, weights = 1/se^2, data = dfz)
+summary(fit_PET)
+#>
+#> Call:
+#> lm(formula = y ~ se, data = dfz, weights = 1/se^2)
+#>
+#> Weighted Residuals:
+#> Min 1Q Median 3Q Max
+#> -3.8132 -0.9112 -0.0139 0.5166 3.3151
+#>
+#> Coefficients:
+#> Estimate Std. Error t value Pr(>|t|)
+#> (Intercept) -0.0008722 0.1081247 -0.008 0.9937
+#> se 2.8549650 1.3593450 2.100 0.0519 .
+#> ---
+#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+#>
+#> Residual standard error: 1.899 on 16 degrees of freedom
+#> Multiple R-squared: 0.2161, Adjusted R-squared: 0.1671
+#> F-statistic: 4.411 on 1 and 16 DF, p-value: 0.05192
The summary()
function allows us to explore details of
+the fitted model. The (Intercept)
coefficient refers to the
+meta-analytic effect size (corrected for the correlation with standard
+errors). Again, it is important to keep in mind that the effect size
+estimate is on the Fisher’s z scale. We obtain the estimate on
+correlation scale with the z2r()
function (we pass the
+estimated effect size using the
+summary(fit_PET)$coefficients["(Intercept)", "Estimate"]
+command, which extracts the estimate from the fitted model, it is
+equivalent to simply pasting the value directly
+z2r(-0.0008722083)
).
Since the Fisher’s z transformation is almost linear around +zero, we obtain an almost identical estimate.
+More importantly, since the test for the effect size with PET was not
+significant at
+,
+we interpret the PET model. However, if the test for effect size were
+significant, we would fit and interpret the PEESE model. The PEESE model
+can be fitted in an analogous way, by replacing the predictor of
+standard errors with standard errors squared (we need to wrap the
+se^2
predictor in I()
that tells R to square
+the predictor prior to fitting the model).
+fit_PEESE <- lm(y ~ I(se^2), weights = 1/se^2, data = dfz)
+summary(fit_PEESE)
+#>
+#> Call:
+#> lm(formula = y ~ I(se^2), data = dfz, weights = 1/se^2)
+#>
+#> Weighted Residuals:
+#> Min 1Q Median 3Q Max
+#> -3.7961 -0.9581 -0.1156 0.6718 3.4608
+#>
+#> Coefficients:
+#> Estimate Std. Error t value Pr(>|t|)
+#> (Intercept) 0.11498 0.06201 1.854 0.0822 .
+#> I(se^2) 15.58064 7.96723 1.956 0.0682 .
+#> ---
+#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+#>
+#> Residual standard error: 1.927 on 16 degrees of freedom
+#> Multiple R-squared: 0.1929, Adjusted R-squared: 0.1425
+#> F-statistic: 3.824 on 1 and 16 DF, p-value: 0.06821
The second publication bias adjustment that we will perform is
+selection models. Selection models adjust for the different publication
+probabilities in different p-value intervals. Selection models
+are implemented in weightr
package
+(weightfunct()
function; Coburn et
+al. (2019)) and newly also in the metafor
package
+(selmodel()
function; Wolfgang
+(2010)). First, we use the weightr
implementation
+and fit the “4PSM” selection model that specifies three distinct
+p-value intervals: (1) covering the range of significant
+p-values for effect sizes in the expected direction
+(0.00-0.025), (2) covering the range of “marginally” significant
+p-values for effect sizes in the expected direction
+(0.025-0.05), and (3) covering the range of non-significant
+p-values (0.05-1). We use Cohen’s d transformation of
+the correlation coefficients since it is better at maintaining the
+distribution of test statistics. To fit the model, we need to pass the
+effect sizes (dfd$y
) into the effect
argument
+and variances (dfd$se^2
) into the v
argument
+(note that we need to pass the vector of values directly since the
+weightfunct()
function does not allow us to pass the data
+frame directly as did the previous functions). We further set
+steps = c(0.025, 0.05)
to specify the appropriate
+cut-points (note that the steps correspond to one-sided
+p-values), and we set table = TRUE
to obtain the
+frequency of p values in each of the specified intervals.
+fit_4PSM <- weightfunct(effect = dfd$y, v = dfd$se^2, steps = c(0.025, 0.05), table = TRUE)
+#> Warning in weightfunct(effect = dfd$y, v = dfd$se^2, steps = c(0.025, 0.05), :
+#> At least one of the p-value intervals contains three or fewer effect sizes,
+#> which may lead to estimation problems. Consider re-specifying the cutpoints.
+fit_4PSM
+#>
+#> Unadjusted Model (k = 18):
+#>
+#> tau^2 (estimated amount of total heterogeneity): 0.0920 (SE = 0.0423)
+#> tau (square root of estimated tau^2 value): 0.3034
+#>
+#> Test for Heterogeneity:
+#> Q(df = 17) = 75.4999, p-val = 5.188348e-09
+#>
+#> Model Results:
+#>
+#> estimate std.error z-stat p-val ci.lb ci.ub
+#> Intercept 0.516 0.08473 6.09 1.1283e-09 0.35 0.6821
+#>
+#> Adjusted Model (k = 18):
+#>
+#> tau^2 (estimated amount of total heterogeneity): 0.1289 (SE = 0.0682)
+#> tau (square root of estimated tau^2 value): 0.3590
+#>
+#> Test for Heterogeneity:
+#> Q(df = 17) = 75.4999, p-val = 5.188348e-09
+#>
+#> Model Results:
+#>
+#> estimate std.error z-stat p-val ci.lb ci.ub
+#> Intercept 0.2675 0.2009 1.3311 0.18316 -0.1264 0.6613
+#> 0.025 < p < 0.05 0.5008 0.5449 0.9191 0.35803 -0.5671 1.5688
+#> 0.05 < p < 1 0.1535 0.1570 0.9777 0.32821 -0.1542 0.4611
+#>
+#> Likelihood Ratio Test:
+#> X^2(df = 2) = 3.844252, p-val = 0.1463
+#>
+#> Number of Effect Sizes per Interval:
+#>
+#> Frequency
+#> p-values <0.025 14
+#> 0.025 < p-values < 0.05 1
+#> 0.05 < p-values < 1 3
Note the warning message informing us about the fact that our data do
+not contain a sufficient number of p-values in one of the
+p-value intervals. The model output obtained by printing the
+fitted model object fit_4PSM
shows that there is only one
+p-value in the (0.025, 0.05) interval. We can deal with this
+issue by joining the “marginally” significant and non-significant
+p-value interval, resulting in the “3PSM” model.
+fit_3PSM <- weightfunct(effect = dfd$y, v = dfd$se^2, steps = c(0.025), table = TRUE)
+fit_3PSM
+#>
+#> Unadjusted Model (k = 18):
+#>
+#> tau^2 (estimated amount of total heterogeneity): 0.0920 (SE = 0.0423)
+#> tau (square root of estimated tau^2 value): 0.3034
+#>
+#> Test for Heterogeneity:
+#> Q(df = 17) = 75.4999, p-val = 5.188348e-09
+#>
+#> Model Results:
+#>
+#> estimate std.error z-stat p-val ci.lb ci.ub
+#> Intercept 0.516 0.08473 6.09 1.1283e-09 0.35 0.6821
+#>
+#> Adjusted Model (k = 18):
+#>
+#> tau^2 (estimated amount of total heterogeneity): 0.1148 (SE = 0.0577)
+#> tau (square root of estimated tau^2 value): 0.3388
+#>
+#> Test for Heterogeneity:
+#> Q(df = 17) = 75.4999, p-val = 5.188348e-09
+#>
+#> Model Results:
+#>
+#> estimate std.error z-stat p-val ci.lb ci.ub
+#> Intercept 0.3220 0.1676 1.921 0.054698 -0.006484 0.6504
+#> 0.025 < p < 1 0.2275 0.2004 1.135 0.256293 -0.165324 0.6204
+#>
+#> Likelihood Ratio Test:
+#> X^2(df = 1) = 3.107176, p-val = 0.077948
+#>
+#> Number of Effect Sizes per Interval:
+#>
+#> Frequency
+#> p-values <0.025 14
+#> 0.025 < p-values < 1 4
The new model does not suffer from the estimation problem due to the
+limited number of p-values in the intervals, so we can now
+interpret the results with more confidence. First, we check the test for
+heterogeneity that clearly rejects the null hypothesis
+Q(df = 17) = 75.4999, $p$ = 5.188348e-09
(if we did not
+find evidence for heterogeneity, we could have proceeded by fitting the
+fixed effects version of the model by specifying the
+fe = TRUE
argument). We follow by checking the test for
+publication bias which is a likelihood ratio test comparing the
+unadjusted and adjusted estimate
+X^2(df = 1) = 3.107176, $p$ = 0.077948
. The result of the
+test is slightly ambiguous – we would reject the null hypothesis of no
+publication bias with
+
+but not with
+.
If we decide to interpret the estimated effect size, we have to again
+transform it back to the correlation scale. However, this time we need
+to use the d2r()
function since we supplied the effect
+sizes as Cohen’s d (note that the effect size estimate
+corresponds to the second value in the fit_3PSM$adj_est
+object for the random effect model, alternatively, we could simply use
+d2r(0.3219641)
).
+d2r(fit_3PSM$adj_est[2])
+#> [1] 0.1589358
Alternatively, we could have conducted the analysis analogously but
+with the metafor
package. First, we would fit a random
+effect meta-analysis with the Cohen’s d transformed effect
+sizes.
+fit_rma_d <- rma(yi = y, sei = se, data = dfd)
Subsequently, we would have used the selmodel
function,
+passing the estimated random effect meta-analysis object and specifying
+the type = "stepfun"
argument to obtain a step weight
+function and setting the appropriate steps with the
+steps = c(0.025)
argument.
+fit_sel_d <- selmodel(fit_rma_d, type = "stepfun", steps = c(0.025))
+fit_sel_d
+#>
+#> Random-Effects Model (k = 18; tau^2 estimator: ML)
+#>
+#> tau^2 (estimated amount of total heterogeneity): 0.1148 (SE = 0.0577)
+#> tau (square root of estimated tau^2 value): 0.3388
+#>
+#> Test for Heterogeneity:
+#> LRT(df = 1) = 32.7499, p-val < .0001
+#>
+#> Model Results:
+#>
+#> estimate se zval pval ci.lb ci.ub
+#> 0.3220 0.1676 1.9214 0.0547 -0.0065 0.6504 .
+#>
+#> Test for Selection Model Parameters:
+#> LRT(df = 1) = 3.1072, p-val = 0.0779
+#>
+#> Selection Model Results:
+#>
+#> k estimate se zval pval ci.lb ci.ub
+#> 0 < p <= 0.025 14 1.0000 --- --- --- --- ---
+#> 0.025 < p <= 1 4 0.2275 0.2004 -3.8537 0.0001 0.0000 0.6204 ***
+#>
+#> ---
+#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The output verifies the results obtained in the previous +analysis.
+The third and final publication bias adjustment that we will perform
+is robust Bayesian meta-analysis (RoBMA). RoBMA uses Bayesian
+model-averaging to combine inference from both PET-PEESE and selection
+models. We use the RoBMA
R package (and the
+RoBMA()
function; Bartoš & Maier
+(2020)) to fit the default 36 model ensemble (called RoBMA-PSMA)
+based on an orthogonal combination of models assuming the presence and
+absence of the effect size, heterogeneity, and publication bias. The
+models assuming the presence of publication bias are further split into
+six weight function models and models utilizing the PET and PEESE
+publication bias adjustment. To fit the model, we can directly pass the
+original correlation coefficients into the r
argument and
+sample sizes into the n
argument – the RoBMA()
+function will internally transform them to the Fisher’s z scale
+and, by default, return the estimates on a Cohen’s d scale
+which is used to specify the prior distributions (both of these settings
+can be changed with the prior_scale
and
+transformation
arguments, and the output can be
+conveniently transformed later). We further set the model
+argument to "PSMA"
to fit the 36 model ensemble and use the
+seed
argument to make the analysis reproducible (it uses
+MCMC sampling in contrast to the previous methods). We turn on parallel
+estimation by setting the parallel = TRUE
argument (the
+parallel processing might in some cases fail, try rerunning the model
+one more time or turning the parallel processing off in that case).
+fit_RoBMA <- RoBMA(r = df$r, n = df$n, seed = 1, model = "PSMA", parallel = TRUE)
This step can take some time depending on your CPU. For example, this +will take approximately 1 minute on a fast CPU (e.g., AMD Ryzen 3900x +12c/24t) and up to ten minutes or longer on slower CPUs (e.g., 2.7 GHz +Intel Core i5).
+We use the summary()
function to explore details of the
+fitted model.
+summary(fit_RoBMA)
+#> Call:
+#> RoBMA(r = df$r, n = df$n, model_type = "PSMA", parallel = TRUE,
+#> save = "min", seed = 1)
+#>
+#> Robust Bayesian meta-analysis
+#> Components summary:
+#> Models Prior prob. Post. prob. Inclusion BF
+#> Effect 18/36 0.500 0.552 1.232
+#> Heterogeneity 18/36 0.500 1.000 19168.311
+#> Bias 32/36 0.500 0.845 5.436
+#>
+#> Model-averaged estimates:
+#> Mean Median 0.025 0.975
+#> mu 0.195 0.087 -0.008 0.598
+#> tau 0.330 0.307 0.166 0.597
+#> omega[0,0.025] 1.000 1.000 1.000 1.000
+#> omega[0.025,0.05] 0.936 1.000 0.438 1.000
+#> omega[0.05,0.5] 0.740 1.000 0.065 1.000
+#> omega[0.5,0.95] 0.697 1.000 0.028 1.000
+#> omega[0.95,0.975] 0.704 1.000 0.028 1.000
+#> omega[0.975,1] 0.713 1.000 0.028 1.000
+#> PET 0.828 0.000 0.000 3.291
+#> PEESE 0.802 0.000 0.000 10.805
+#> The estimates are summarized on the Cohen's d scale (priors were specified on the Cohen's d scale).
+#> (Estimated publication weights omega correspond to one-sided p-values.)
The printed output consists of two parts. The first table called
+Components summary
contains information about the fitted
+models. It tells us that we estimated the ensemble with 18/36 models
+assuming the presence of an effect, 18/36 models assuming the presence
+of heterogeneity, and 32/36 models assuming the presence of the
+publication bias. The second column summarizes the prior model
+probabilities of models assuming either presence of the individual
+components – here, we see that the presence and absence of the
+components are balanced a priori. The third column contains information
+about the posterior probability of models assuming the presence of the
+components – we can observe that the posterior model probabilities of
+models assuming the presence of an effect slightly increased to 0.552.
+The last column contains information about the evidence in favor of the
+presence of any of those components. Evidence for the presence of an
+effect is undecided; the models assuming the presence of an effect are
+only 1.232 times more likely given the data than the models assuming the
+absence of an effect. However, we find overwhelming evidence in favor of
+heterogeneity, with the models assuming the presence of heterogeneity
+being 19,168 times more likely given the data than models assuming the
+absence of heterogeneity, and moderate evidence in favor of publication
+bias.
As the name indicates, the second table called
+Model-averaged estimates
contains information about the
+model-averaged estimates. The first row labeled mu
+corresponds to the model-averaged effect size estimate (on Cohen’s
+d scale) and the second row labeled tau
+corresponds to the model-averaged heterogeneity estimates. Below are the
+estimated model-averaged weights for the different p-value
+intervals and the PET and PEESE regression coefficients. We convert the
+estimates to the correlation coefficients by adding the
+output_scale = "r"
argument to the summary function.
+summary(fit_RoBMA, output_scale = "r")
+#> Call:
+#> RoBMA(r = df$r, n = df$n, model_type = "PSMA", parallel = TRUE,
+#> save = "min", seed = 1)
+#>
+#> Robust Bayesian meta-analysis
+#> Components summary:
+#> Models Prior prob. Post. prob. Inclusion BF
+#> Effect 18/36 0.500 0.552 1.232
+#> Heterogeneity 18/36 0.500 1.000 19168.311
+#> Bias 32/36 0.500 0.845 5.436
+#>
+#> Model-averaged estimates:
+#> Mean Median 0.025 0.975
+#> mu 0.095 0.043 -0.004 0.286
+#> tau 0.165 0.154 0.083 0.299
+#> omega[0,0.025] 1.000 1.000 1.000 1.000
+#> omega[0.025,0.05] 0.936 1.000 0.438 1.000
+#> omega[0.05,0.5] 0.740 1.000 0.065 1.000
+#> omega[0.5,0.95] 0.697 1.000 0.028 1.000
+#> omega[0.95,0.975] 0.704 1.000 0.028 1.000
+#> omega[0.975,1] 0.713 1.000 0.028 1.000
+#> PET 0.828 0.000 0.000 3.291
+#> PEESE 1.603 0.000 0.000 21.610
+#> The effect size estimates are summarized on the correlation scale and heterogeneity is summarized on the Fisher's z scale (priors were specified on the Cohen's d scale).
+#> (Estimated publication weights omega correspond to one-sided p-values.)
Now, we have obtained the model-averaged effect size estimate on the
+correlation scale. If we were interested in the estimates
+model-averaging only across the models assuming the presence of an
+effect (for the effect size estimate), heterogeneity (for the
+heterogeneity estimate), and publication bias (for the publication bias
+weights and PET and PEESE regression coefficients), we could have added
+the conditional = TRUE
argument to the summary function. A
+quick textual summary of the model can also be generated with the
+interpret()
function.
+interpret(fit_RoBMA, output_scale = "r")
+#> [1] "Robust Bayesian meta-analysis found weak evidence in favor of the effect, BF_10 = 1.23, with mean model-averaged estimate correlation = 0.095, 95% CI [-0.004, 0.286]. Robust Bayesian meta-analysis found strong evidence in favor of the heterogeneity, BF^rf = 19168.31, with mean model-averaged estimate tau = 0.165, 95% CI [0.083, 0.299]. Robust Bayesian meta-analysis found moderate evidence in favor of the publication bias, BF_pb = 5.44."
We can also obtain summary information about the individual models by
+specifying the type = "models"
option. The resulting table
+shows the prior and posterior model probabilities and inclusion Bayes
+factors for the individual models (we also set the
+short_name = TRUE
argument reducing the width of the output
+by abbreviating names of the prior distributions).
+summary(fit_RoBMA, type = "models", short_name = TRUE)
+#> Call:
+#> RoBMA(r = df$r, n = df$n, model_type = "PSMA", parallel = TRUE,
+#> save = "min", seed = 1)
+#>
+#> Robust Bayesian meta-analysis
+#> Models overview:
+#> Model Prior Effect Prior Heterogeneity
+#> 1 S(0) S(0)
+#> 2 S(0) S(0)
+#> 3 S(0) S(0)
+#> 4 S(0) S(0)
+#> 5 S(0) S(0)
+#> 6 S(0) S(0)
+#> 7 S(0) S(0)
+#> 8 S(0) S(0)
+#> 9 S(0) S(0)
+#> 10 S(0) Ig(1, 0.15)
+#> 11 S(0) Ig(1, 0.15)
+#> 12 S(0) Ig(1, 0.15)
+#> 13 S(0) Ig(1, 0.15)
+#> 14 S(0) Ig(1, 0.15)
+#> 15 S(0) Ig(1, 0.15)
+#> 16 S(0) Ig(1, 0.15)
+#> 17 S(0) Ig(1, 0.15)
+#> 18 S(0) Ig(1, 0.15)
+#> 19 N(0, 1) S(0)
+#> 20 N(0, 1) S(0)
+#> 21 N(0, 1) S(0)
+#> 22 N(0, 1) S(0)
+#> 23 N(0, 1) S(0)
+#> 24 N(0, 1) S(0)
+#> 25 N(0, 1) S(0)
+#> 26 N(0, 1) S(0)
+#> 27 N(0, 1) S(0)
+#> 28 N(0, 1) Ig(1, 0.15)
+#> 29 N(0, 1) Ig(1, 0.15)
+#> 30 N(0, 1) Ig(1, 0.15)
+#> 31 N(0, 1) Ig(1, 0.15)
+#> 32 N(0, 1) Ig(1, 0.15)
+#> 33 N(0, 1) Ig(1, 0.15)
+#> 34 N(0, 1) Ig(1, 0.15)
+#> 35 N(0, 1) Ig(1, 0.15)
+#> 36 N(0, 1) Ig(1, 0.15)
+#> Prior Bias Prior prob. log(marglik)
+#> 0.125 -74.67
+#> omega[2s: .05] ~ CumD(1, 1) 0.010 -49.60
+#> omega[2s: .1, .05] ~ CumD(1, 1, 1) 0.010 -47.53
+#> omega[1s: .05] ~ CumD(1, 1) 0.010 -41.70
+#> omega[1s: .05, .025] ~ CumD(1, 1, 1) 0.010 -38.03
+#> omega[1s: .5, .05] ~ CumD(1, 1, 1) 0.010 -44.41
+#> omega[1s: .5, .05, .025] ~ CumD(1, 1, 1, 1) 0.010 -40.79
+#> PET ~ C(0, 1)[0, Inf] 0.031 -5.01
+#> PEESE ~ C(0, 5)[0, Inf] 0.031 -12.17
+#> 0.125 -6.95
+#> omega[2s: .05] ~ CumD(1, 1) 0.010 -5.96
+#> omega[2s: .1, .05] ~ CumD(1, 1, 1) 0.010 -5.09
+#> omega[1s: .05] ~ CumD(1, 1) 0.010 2.72
+#> omega[1s: .05, .025] ~ CumD(1, 1, 1) 0.010 2.93
+#> omega[1s: .5, .05] ~ CumD(1, 1, 1) 0.010 2.91
+#> omega[1s: .5, .05, .025] ~ CumD(1, 1, 1, 1) 0.010 3.30
+#> PET ~ C(0, 1)[0, Inf] 0.031 3.62
+#> PEESE ~ C(0, 5)[0, Inf] 0.031 1.62
+#> 0.125 -13.17
+#> omega[2s: .05] ~ CumD(1, 1) 0.010 -13.10
+#> omega[2s: .1, .05] ~ CumD(1, 1, 1) 0.010 -12.87
+#> omega[1s: .05] ~ CumD(1, 1) 0.010 -12.75
+#> omega[1s: .05, .025] ~ CumD(1, 1, 1) 0.010 -12.86
+#> omega[1s: .5, .05] ~ CumD(1, 1, 1) 0.010 -13.29
+#> omega[1s: .5, .05, .025] ~ CumD(1, 1, 1, 1) 0.010 -13.25
+#> PET ~ C(0, 1)[0, Inf] 0.031 -7.07
+#> PEESE ~ C(0, 5)[0, Inf] 0.031 -7.58
+#> 0.125 1.79
+#> omega[2s: .05] ~ CumD(1, 1) 0.010 1.75
+#> omega[2s: .1, .05] ~ CumD(1, 1, 1) 0.010 2.16
+#> omega[1s: .05] ~ CumD(1, 1) 0.010 3.11
+#> omega[1s: .05, .025] ~ CumD(1, 1, 1) 0.010 3.01
+#> omega[1s: .5, .05] ~ CumD(1, 1, 1) 0.010 2.98
+#> omega[1s: .5, .05, .025] ~ CumD(1, 1, 1, 1) 0.010 3.06
+#> PET ~ C(0, 1)[0, Inf] 0.031 2.75
+#> PEESE ~ C(0, 5)[0, Inf] 0.031 2.55
+#> Post. prob. Inclusion BF
+#> 0.000 0.000
+#> 0.000 0.000
+#> 0.000 0.000
+#> 0.000 0.000
+#> 0.000 0.000
+#> 0.000 0.000
+#> 0.000 0.000
+#> 0.000 0.001
+#> 0.000 0.000
+#> 0.000 0.000
+#> 0.000 0.001
+#> 0.000 0.001
+#> 0.033 3.231
+#> 0.041 4.025
+#> 0.040 3.919
+#> 0.059 5.927
+#> 0.243 9.957
+#> 0.033 1.055
+#> 0.000 0.000
+#> 0.000 0.000
+#> 0.000 0.000
+#> 0.000 0.000
+#> 0.000 0.000
+#> 0.000 0.000
+#> 0.000 0.000
+#> 0.000 0.000
+#> 0.000 0.000
+#> 0.155 1.287
+#> 0.012 1.201
+#> 0.019 1.822
+#> 0.048 4.831
+#> 0.044 4.347
+#> 0.043 4.223
+#> 0.046 4.617
+#> 0.102 3.504
+#> 0.083 2.797
To obtain a summary of the individual model diagnostics, we set the
+type = "diagnostics"
argument. The resulting table provides
+information about the maximum MCMC error, relative MCMC error, minimum
+ESS, and maximum R-hat when aggregating over the parameters of each
+model. As we can see, we obtain acceptable ESS and R-hat diagnostic
+values.
+summary(fit_RoBMA, type = "diagnostics")
+#> Call:
+#> RoBMA(r = df$r, n = df$n, model_type = "PSMA", parallel = TRUE,
+#> save = "min", seed = 1)
+#>
+#> Robust Bayesian meta-analysis
+#> Diagnostics overview:
+#> Model Prior Effect Prior Heterogeneity
+#> 1 Spike(0) Spike(0)
+#> 2 Spike(0) Spike(0)
+#> 3 Spike(0) Spike(0)
+#> 4 Spike(0) Spike(0)
+#> 5 Spike(0) Spike(0)
+#> 6 Spike(0) Spike(0)
+#> 7 Spike(0) Spike(0)
+#> 8 Spike(0) Spike(0)
+#> 9 Spike(0) Spike(0)
+#> 10 Spike(0) InvGamma(1, 0.15)
+#> 11 Spike(0) InvGamma(1, 0.15)
+#> 12 Spike(0) InvGamma(1, 0.15)
+#> 13 Spike(0) InvGamma(1, 0.15)
+#> 14 Spike(0) InvGamma(1, 0.15)
+#> 15 Spike(0) InvGamma(1, 0.15)
+#> 16 Spike(0) InvGamma(1, 0.15)
+#> 17 Spike(0) InvGamma(1, 0.15)
+#> 18 Spike(0) InvGamma(1, 0.15)
+#> 19 Normal(0, 1) Spike(0)
+#> 20 Normal(0, 1) Spike(0)
+#> 21 Normal(0, 1) Spike(0)
+#> 22 Normal(0, 1) Spike(0)
+#> 23 Normal(0, 1) Spike(0)
+#> 24 Normal(0, 1) Spike(0)
+#> 25 Normal(0, 1) Spike(0)
+#> 26 Normal(0, 1) Spike(0)
+#> 27 Normal(0, 1) Spike(0)
+#> 28 Normal(0, 1) InvGamma(1, 0.15)
+#> 29 Normal(0, 1) InvGamma(1, 0.15)
+#> 30 Normal(0, 1) InvGamma(1, 0.15)
+#> 31 Normal(0, 1) InvGamma(1, 0.15)
+#> 32 Normal(0, 1) InvGamma(1, 0.15)
+#> 33 Normal(0, 1) InvGamma(1, 0.15)
+#> 34 Normal(0, 1) InvGamma(1, 0.15)
+#> 35 Normal(0, 1) InvGamma(1, 0.15)
+#> 36 Normal(0, 1) InvGamma(1, 0.15)
+#> Prior Bias max[error(MCMC)]
+#> NA
+#> omega[two-sided: .05] ~ CumDirichlet(1, 1) 0.00024
+#> omega[two-sided: .1, .05] ~ CumDirichlet(1, 1, 1) 0.00295
+#> omega[one-sided: .05] ~ CumDirichlet(1, 1) 0.00014
+#> omega[one-sided: .05, .025] ~ CumDirichlet(1, 1, 1) 0.00326
+#> omega[one-sided: .5, .05] ~ CumDirichlet(1, 1, 1) 0.00033
+#> omega[one-sided: .5, .05, .025] ~ CumDirichlet(1, 1, 1, 1) 0.00309
+#> PET ~ Cauchy(0, 1)[0, Inf] 0.00236
+#> PEESE ~ Cauchy(0, 5)[0, Inf] 0.01223
+#> 0.00118
+#> omega[two-sided: .05] ~ CumDirichlet(1, 1) 0.00296
+#> omega[two-sided: .1, .05] ~ CumDirichlet(1, 1, 1) 0.00295
+#> omega[one-sided: .05] ~ CumDirichlet(1, 1) 0.00110
+#> omega[one-sided: .05, .025] ~ CumDirichlet(1, 1, 1) 0.00331
+#> omega[one-sided: .5, .05] ~ CumDirichlet(1, 1, 1) 0.00357
+#> omega[one-sided: .5, .05, .025] ~ CumDirichlet(1, 1, 1, 1) 0.00307
+#> PET ~ Cauchy(0, 1)[0, Inf] 0.00454
+#> PEESE ~ Cauchy(0, 5)[0, Inf] 0.02470
+#> 0.00038
+#> omega[two-sided: .05] ~ CumDirichlet(1, 1) 0.00303
+#> omega[two-sided: .1, .05] ~ CumDirichlet(1, 1, 1) 0.00290
+#> omega[one-sided: .05] ~ CumDirichlet(1, 1) 0.00309
+#> omega[one-sided: .05, .025] ~ CumDirichlet(1, 1, 1) 0.00278
+#> omega[one-sided: .5, .05] ~ CumDirichlet(1, 1, 1) 0.00332
+#> omega[one-sided: .5, .05, .025] ~ CumDirichlet(1, 1, 1, 1) 0.00293
+#> PET ~ Cauchy(0, 1)[0, Inf] 0.03247
+#> PEESE ~ Cauchy(0, 5)[0, Inf] 0.05228
+#> 0.00090
+#> omega[two-sided: .05] ~ CumDirichlet(1, 1) 0.00308
+#> omega[two-sided: .1, .05] ~ CumDirichlet(1, 1, 1) 0.00293
+#> omega[one-sided: .05] ~ CumDirichlet(1, 1) 0.00477
+#> omega[one-sided: .05, .025] ~ CumDirichlet(1, 1, 1) 0.00340
+#> omega[one-sided: .5, .05] ~ CumDirichlet(1, 1, 1) 0.00543
+#> omega[one-sided: .5, .05, .025] ~ CumDirichlet(1, 1, 1, 1) 0.00499
+#> PET ~ Cauchy(0, 1)[0, Inf] 0.04070
+#> PEESE ~ Cauchy(0, 5)[0, Inf] 0.07238
+#> max[error(MCMC)/SD] min(ESS) max(R-hat)
+#> NA NA NA
+#> 0.016 4158 1.000
+#> 0.016 3793 1.000
+#> 0.015 4622 1.000
+#> 0.017 3357 1.000
+#> 0.017 3509 1.001
+#> 0.018 3064 1.001
+#> 0.010 9917 1.001
+#> 0.010 9589 1.000
+#> 0.010 9632 1.001
+#> 0.013 5518 1.002
+#> 0.015 4565 1.001
+#> 0.015 4395 1.001
+#> 0.015 4502 1.002
+#> 0.018 3206 1.001
+#> 0.017 3480 1.001
+#> 0.012 7342 1.001
+#> 0.012 7051 1.000
+#> 0.010 9712 1.001
+#> 0.013 5522 1.000
+#> 0.015 4382 1.001
+#> 0.013 5771 1.000
+#> 0.014 4859 1.001
+#> 0.015 4430 1.000
+#> 0.016 4135 1.001
+#> 0.042 565 1.005
+#> 0.024 1678 1.001
+#> 0.011 7736 1.000
+#> 0.014 5254 1.001
+#> 0.016 4103 1.001
+#> 0.021 2240 1.001
+#> 0.020 2527 1.001
+#> 0.026 1529 1.007
+#> 0.024 1756 1.000
+#> 0.038 692 1.001
+#> 0.024 1765 1.005
Finally, we can also plot the model-averaged posterior distribution
+with the plot()
function. We set the
+prior = TRUE
argument to include the prior distribution as
+a grey line (and arrow for the point density at zero) and
+output_scale = "r"
to transform the posterior distribution
+to the correlation scale (the default figure output would be on Cohen’s
+d scale). (The par(mar = c(4, 4, 1, 4))
call
+increases the left margin of the figure, so the secondary y-axis text is
+not cut off.)
The RoBMA
package allows us to fit ensembles of highly
+customized meta-analytic models. Here we reproduce the ensemble for the
+perinull directional hypothesis test from the Appendix (see the R
+package vignettes for more examples and details). Instead of using the
+fully pre-specified model with the model = "PSMA"
argument,
+we explicitly specify the prior distribution for models assuming
+presence of the effect with the
+priors_effect = prior("normal", parameters = list(mean = 0.60, sd = 0.20), truncation = list(0, Inf))
+argument, which assigns Normal(0.60, 0.20) distribution bounded to the
+positive numbers to the
+
+parameter (note that the prior distribution is specified on the Cohen’s
+d scale, corresponding to 95% prior probability mass contained
+approximately in the
+
+= (0.10, 0.45) interval). Similarly, we also replace the default prior
+distribution for the models assuming absence of the effect with a
+perinull hypothesis with the
+priors_effect_null = prior("normal", parameters = list(mean = 0, sd = 0.10))
+argument that sets 95% prior probability mass to values in the
+
+= (-0.10, 0.10) interval.
+fit_RoBMA2 <- RoBMA(r = df$r, n = df$n, seed = 2, parallel = TRUE,
+ priors_effect = prior("normal", parameters = list(mean = 0.60, sd = 0.20), truncation = list(0, Inf)),
+ priors_effect_null = prior("normal", parameters = list(mean = 0, sd = 0.10)))
As previously, we can use the summary()
function to
+inspect the model fit and verify that the specified models correspond to
+the settings.
+summary(fit_RoBMA2, type = "models")
+#> Call:
+#> RoBMA(r = df$r, n = df$n, priors_effect = prior("normal", parameters = list(mean = 0.6,
+#> sd = 0.2), truncation = list(0, Inf)), priors_effect_null = prior("normal",
+#> parameters = list(mean = 0, sd = 0.1)), parallel = TRUE,
+#> save = "min", seed = 2)
+#>
+#> Robust Bayesian meta-analysis
+#> Models overview:
+#> Model Prior Effect Prior Heterogeneity
+#> 1 Normal(0, 0.1) Spike(0)
+#> 2 Normal(0, 0.1) Spike(0)
+#> 3 Normal(0, 0.1) Spike(0)
+#> 4 Normal(0, 0.1) Spike(0)
+#> 5 Normal(0, 0.1) Spike(0)
+#> 6 Normal(0, 0.1) Spike(0)
+#> 7 Normal(0, 0.1) Spike(0)
+#> 8 Normal(0, 0.1) Spike(0)
+#> 9 Normal(0, 0.1) Spike(0)
+#> 10 Normal(0, 0.1) InvGamma(1, 0.15)
+#> 11 Normal(0, 0.1) InvGamma(1, 0.15)
+#> 12 Normal(0, 0.1) InvGamma(1, 0.15)
+#> 13 Normal(0, 0.1) InvGamma(1, 0.15)
+#> 14 Normal(0, 0.1) InvGamma(1, 0.15)
+#> 15 Normal(0, 0.1) InvGamma(1, 0.15)
+#> 16 Normal(0, 0.1) InvGamma(1, 0.15)
+#> 17 Normal(0, 0.1) InvGamma(1, 0.15)
+#> 18 Normal(0, 0.1) InvGamma(1, 0.15)
+#> 19 Normal(0.6, 0.2)[0, Inf] Spike(0)
+#> 20 Normal(0.6, 0.2)[0, Inf] Spike(0)
+#> 21 Normal(0.6, 0.2)[0, Inf] Spike(0)
+#> 22 Normal(0.6, 0.2)[0, Inf] Spike(0)
+#> 23 Normal(0.6, 0.2)[0, Inf] Spike(0)
+#> 24 Normal(0.6, 0.2)[0, Inf] Spike(0)
+#> 25 Normal(0.6, 0.2)[0, Inf] Spike(0)
+#> 26 Normal(0.6, 0.2)[0, Inf] Spike(0)
+#> 27 Normal(0.6, 0.2)[0, Inf] Spike(0)
+#> 28 Normal(0.6, 0.2)[0, Inf] InvGamma(1, 0.15)
+#> 29 Normal(0.6, 0.2)[0, Inf] InvGamma(1, 0.15)
+#> 30 Normal(0.6, 0.2)[0, Inf] InvGamma(1, 0.15)
+#> 31 Normal(0.6, 0.2)[0, Inf] InvGamma(1, 0.15)
+#> 32 Normal(0.6, 0.2)[0, Inf] InvGamma(1, 0.15)
+#> 33 Normal(0.6, 0.2)[0, Inf] InvGamma(1, 0.15)
+#> 34 Normal(0.6, 0.2)[0, Inf] InvGamma(1, 0.15)
+#> 35 Normal(0.6, 0.2)[0, Inf] InvGamma(1, 0.15)
+#> 36 Normal(0.6, 0.2)[0, Inf] InvGamma(1, 0.15)
+#> Prior Bias Prior prob.
+#> 0.125
+#> omega[two-sided: .05] ~ CumDirichlet(1, 1) 0.010
+#> omega[two-sided: .1, .05] ~ CumDirichlet(1, 1, 1) 0.010
+#> omega[one-sided: .05] ~ CumDirichlet(1, 1) 0.010
+#> omega[one-sided: .05, .025] ~ CumDirichlet(1, 1, 1) 0.010
+#> omega[one-sided: .5, .05] ~ CumDirichlet(1, 1, 1) 0.010
+#> omega[one-sided: .5, .05, .025] ~ CumDirichlet(1, 1, 1, 1) 0.010
+#> PET ~ Cauchy(0, 1)[0, Inf] 0.031
+#> PEESE ~ Cauchy(0, 5)[0, Inf] 0.031
+#> 0.125
+#> omega[two-sided: .05] ~ CumDirichlet(1, 1) 0.010
+#> omega[two-sided: .1, .05] ~ CumDirichlet(1, 1, 1) 0.010
+#> omega[one-sided: .05] ~ CumDirichlet(1, 1) 0.010
+#> omega[one-sided: .05, .025] ~ CumDirichlet(1, 1, 1) 0.010
+#> omega[one-sided: .5, .05] ~ CumDirichlet(1, 1, 1) 0.010
+#> omega[one-sided: .5, .05, .025] ~ CumDirichlet(1, 1, 1, 1) 0.010
+#> PET ~ Cauchy(0, 1)[0, Inf] 0.031
+#> PEESE ~ Cauchy(0, 5)[0, Inf] 0.031
+#> 0.125
+#> omega[two-sided: .05] ~ CumDirichlet(1, 1) 0.010
+#> omega[two-sided: .1, .05] ~ CumDirichlet(1, 1, 1) 0.010
+#> omega[one-sided: .05] ~ CumDirichlet(1, 1) 0.010
+#> omega[one-sided: .05, .025] ~ CumDirichlet(1, 1, 1) 0.010
+#> omega[one-sided: .5, .05] ~ CumDirichlet(1, 1, 1) 0.010
+#> omega[one-sided: .5, .05, .025] ~ CumDirichlet(1, 1, 1, 1) 0.010
+#> PET ~ Cauchy(0, 1)[0, Inf] 0.031
+#> PEESE ~ Cauchy(0, 5)[0, Inf] 0.031
+#> 0.125
+#> omega[two-sided: .05] ~ CumDirichlet(1, 1) 0.010
+#> omega[two-sided: .1, .05] ~ CumDirichlet(1, 1, 1) 0.010
+#> omega[one-sided: .05] ~ CumDirichlet(1, 1) 0.010
+#> omega[one-sided: .05, .025] ~ CumDirichlet(1, 1, 1) 0.010
+#> omega[one-sided: .5, .05] ~ CumDirichlet(1, 1, 1) 0.010
+#> omega[one-sided: .5, .05, .025] ~ CumDirichlet(1, 1, 1, 1) 0.010
+#> PET ~ Cauchy(0, 1)[0, Inf] 0.031
+#> PEESE ~ Cauchy(0, 5)[0, Inf] 0.031
+#> log(marglik) Post. prob. Inclusion BF
+#> -18.84 0.000 0.000
+#> -17.66 0.000 0.000
+#> -17.06 0.000 0.000
+#> -17.35 0.000 0.000
+#> -17.04 0.000 0.000
+#> -18.11 0.000 0.000
+#> -17.69 0.000 0.000
+#> -5.24 0.000 0.000
+#> -7.61 0.000 0.000
+#> -3.20 0.000 0.003
+#> -1.45 0.000 0.022
+#> -0.42 0.001 0.061
+#> 3.01 0.020 1.939
+#> 3.19 0.024 2.317
+#> 3.09 0.022 2.104
+#> 3.46 0.031 3.062
+#> 3.64 0.112 3.909
+#> 2.35 0.031 0.986
+#> -11.84 0.000 0.000
+#> -11.88 0.000 0.000
+#> -11.71 0.000 0.000
+#> -11.54 0.000 0.000
+#> -11.70 0.000 0.000
+#> -12.05 0.000 0.000
+#> -12.07 0.000 0.000
+#> -8.38 0.000 0.000
+#> -7.36 0.000 0.000
+#> 3.35 0.337 3.564
+#> 3.13 0.023 2.190
+#> 3.42 0.030 2.951
+#> 4.12 0.061 6.123
+#> 3.85 0.046 4.602
+#> 3.94 0.050 5.027
+#> 3.84 0.046 4.572
+#> 3.23 0.074 2.492
+#> 3.44 0.092 3.132