diff --git a/.gitignore b/.gitignore index 9b2a8e65..b07cabf0 100644 --- a/.gitignore +++ b/.gitignore @@ -15,3 +15,5 @@ Meta RoBMA.Rcheck .Rhistory .Rprofile +/doc/ +/Meta/ diff --git a/DESCRIPTION b/DESCRIPTION index bd2ee829..5fee6d7b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: RoBMA Title: Robust Bayesian Meta-Analyses -Version: 2.1.0 +Version: 2.1.1 Maintainer: František Bartoš Authors@R: c( person("František", "Bartoš", role = c("aut", "cre"), @@ -40,7 +40,7 @@ NeedsCompilation: yes Depends: R (>= 4.0.0) Imports: - BayesTools (>= 0.1.2), + BayesTools (>= 0.1.3), runjags, bridgesampling, rjags, diff --git a/NAMESPACE b/NAMESPACE index 3acc7aae..058aa834 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -35,6 +35,7 @@ export(plot_models) export(prior) export(prior_PEESE) export(prior_PET) +export(prior_informed) export(prior_none) export(prior_weightfunction) export(pwnorm) diff --git a/NEWS.md b/NEWS.md index e8bd7b5b..bbe03a40 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,13 @@ +## version 2.1.1 +### Fixes +- incorrectly formatted citations in vignettes and capitalization + +### Features +- adding `informed_prior()` function (from the BayesTools package) that allows specification of various informed prior distributions from the field of medicine and psychology +- adding a vignette reproducing the example of dentine sensitivity with the informed Bayesian model-averaged meta-analysis from Bartoš et al., 2021 ([open-access](https://onlinelibrary.wiley.com/doi/10.1002/sim.9170)), +- further reductions of fitted object size when setting `save = "min"` + + ## version 2.1 ### Fixes - more informative error message when the JAGS module fails to load diff --git a/R/RoBMA-package.R b/R/RoBMA-package.R index a6a72cd2..f75a8f08 100644 --- a/R/RoBMA-package.R +++ b/R/RoBMA-package.R @@ -17,8 +17,9 @@ ##' methodology. ##' ##' More details regarding customization of the model ensembles are provided in the -##' \href{../doc/ReproducingBMA.html}{\bold{Reproducing BMA}} and -##' \href{../doc/CustomEnsembles.html}{\bold{Fitting custom meta-analytic ensembles}} +##' \href{../doc/ReproducingBMA.html}{\bold{Reproducing BMA}}, +##' \href{../doc/MedicineBMA.html}{\bold{BMA in Medicine}}, and +##' \href{../doc/CustomEnsembles.html}{\bold{Fitting Custom Meta-Analytic Ensembles}} ##' vignettes. Please, use the "Issues" section in the GitHub repository to ask any ##' further questions. ##' diff --git a/R/data.R b/R/data.R index 054f66ae..e2567332 100644 --- a/R/data.R +++ b/R/data.R @@ -1,4 +1,4 @@ -#' @title 27 experimental studies from from +#' @title 27 experimental studies from #' \insertCite{anderson2010violent;textual}{RoBMA} that meet the best practice criteria #' #' @description The data set contains correlation coefficients, sample @@ -16,7 +16,7 @@ "Anderson2010" -#' @title 9 experimental studies from from +#' @title 9 experimental studies from #' \insertCite{bem2011feeling;textual}{RoBMA} as described in #' \insertCite{bem2011must;textual}{RoBMA} #' @@ -32,3 +32,22 @@ #' @references #' \insertAllCited{} "Bem2011" + + +#' @title 5 studies with a tactile outcome assessment from +#' \insertCite{poulsen2006potassium;textual}{RoBMA} of the effect of potassium-containing toothpaste +#' on dentine hypersensitivity +#' +#' @description The data set contains Cohen's d effect sizes, standard errors, +#' and labels for 5 studies assessing the tactile outcome from a meta-analysis of +#' the effect of potassium-containing toothpaste on dentine hypersensitivity +#' \insertCite{poulsen2006potassium}{RoBMA} which was used as an example in +#' \insertCite{bartos2021bayesian;textual}{RoBMA}. +#' +#' @format A data.frame with 3 columns and 5 observations. +#' +#' @return a data.frame. +#' +#' @references +#' \insertAllCited{} +"Poulsen2006" diff --git a/R/main.R b/R/main.R index e8b1d663..90656238 100644 --- a/R/main.R +++ b/R/main.R @@ -297,6 +297,7 @@ RoBMA <- function( ### remove model posteriors if asked to if(save == "min"){ object <- .remove_model_posteriors(object) + object <- .remove_model_margliks(object) } @@ -499,6 +500,7 @@ update.RoBMA <- function(object, refit_failed = TRUE, ### remove model posteriors if asked to if(save == "min"){ object <- .remove_model_posteriors(object) + object <- .remove_model_margliks(object) } return(object) diff --git a/R/priors.R b/R/priors.R index c87a0c85..452f2aa9 100644 --- a/R/priors.R +++ b/R/priors.R @@ -18,9 +18,17 @@ prior_PET <- BayesTools::prior_PET #' @export prior_PEESE <- BayesTools::prior_PEESE - #' @name prior_weightfunction #' @inherit BayesTools::prior_weightfunction #' @export prior_weightfunction <- BayesTools::prior_weightfunction +#' @name prior_informed +#' @inherit BayesTools::prior_informed +#' @details Further details can be found in \insertCite{erp2017estimates;textual}{RoBMA}, +#' \insertCite{gronau2017bayesian;textual}{RoBMA}, and +#' \insertCite{bartos2021bayesian;textual}{RoBMA}. +#' @references +#' \insertAllCited{} +#' @export +prior_informed <- BayesTools::prior_informed diff --git a/R/tools.R b/R/tools.R index 2d680a67..e663859f 100644 --- a/R/tools.R +++ b/R/tools.R @@ -17,7 +17,7 @@ check_RoBMA <- function(fit){ .is_model_constant <- function(priors){ # checks whether there is at least one non-nill prior - return(all(sapply(priors, function(prior)is.prior.point(prior) | is.prior.none(prior))) & is.null(priors[["omega"]])) + return(all(sapply(priors, function(prior) is.prior.point(prior) | is.prior.none(prior))) & is.null(priors[["omega"]])) } .remove_model_posteriors <- function(object){ for(i in seq_along(object[["models"]])){ @@ -27,6 +27,17 @@ check_RoBMA <- function(fit){ } return(object) } +.remove_model_margliks <- function(object){ + for(i in seq_along(object[["models"]])){ + if(inherits(object$models[[i]][["marglik"]], "bridge")){ + object$models[[i]]$marglik[["q11"]] <- NULL + object$models[[i]]$marglik[["q12"]] <- NULL + object$models[[i]]$marglik[["q21"]] <- NULL + object$models[[i]]$marglik[["q22"]] <- NULL + } + } + return(object) +} .print_errors_and_warnings <- function(object, max_print = 5){ errors_and_warnings <- .collect_errors_and_warnings(object, max_print = max_print) diff --git a/data/Poulsen2006.RData b/data/Poulsen2006.RData new file mode 100644 index 00000000..3dd008ae Binary files /dev/null and b/data/Poulsen2006.RData differ diff --git a/inst/REFERENCES.bib b/inst/REFERENCES.bib index cc37915b..b5505961 100644 --- a/inst/REFERENCES.bib +++ b/inst/REFERENCES.bib @@ -6,7 +6,7 @@ @article{anderson2010violent journal = {Psychological Bulletin}, year = {2010}, pages = {151}, - url = {https://doi.org/10.1037/a0018251} + doi = {10.1037/a0018251} } @article{gronau2017bayesian, @@ -18,15 +18,19 @@ @article{gronau2017bayesian pages = {123--138}, year = {2017}, publisher = {Taylor \& Francis}, - url = {https://doi.org/10.1080/23743603.2017.1326760} + doi = {10.1080/23743603.2017.1326760} } @article{gronau2020primer, title = {A primer on {B}ayesian model-averaged meta-analysis}, - author = {Gronau, Quentin Frederik and Heck, Daniel W and Berkhout, Sophie W and Haaf, Julia M and Wagenmakers, Eric-Jan}, + author = {Gronau, Quentin F and Heck, Daniel W and Berkhout, Sophie W and Haaf, Julia M and Wagenmakers, Eric-Jan}, journal = {Advances in Methods and Practices in Psychological Science}, - year = {in press}, - url = {https://doi.org/10.31234/osf.io/97qup} + volume = {4}, + number = {3}, + pages = {25152459211031256}, + year = {2021}, + publisher = {SAGE Publications Sage CA: Los Angeles, CA}, + doi = {10.1177/25152459211031256} } @article{erp2017estimates, @@ -37,7 +41,7 @@ @article{erp2017estimates number = {1}, year = {2017}, publisher = {Ubiquity Press}, - url = {https://doi.org/10.5334/jopd.33} + doi = {10.5334/jopd.33} } @article{maier2020robust, @@ -45,7 +49,7 @@ @article{maier2020robust author = {Maier, Maximilian and Barto{\v{s}}, Franti{\v{s}}ek and Wagenmakers, Eric-Jan}, year = {in press}, journal = {Psychological Methods}, - url = {https://doi.org/10.31234/osf.io/u4cns} + doi = {10.31234/osf.io/u4cns} } @@ -55,16 +59,16 @@ @unpublished{bartos2021no year = {2021}, publisher = {PsyArXiv}, note = {preprint at \url{https://doi.org/10.31234/osf.io/kvsp7}}, - url = {https://doi.org/10.31234/osf.io/kvsp7} + doi = {10.31234/osf.io/kvsp7} } @unpublished{bartos2020adjusting, - title = {Adjusting for Publication Bias in {JASP} — Selection Models and Robust {B}ayesian Meta-Analysis}, - author = {Barto{\v{s}}, Franti{\v{s}}ek and Maier, Maximilian and Wagenmakers, Eric-Jan}, - year = {2020}, + title = {Adjusting for publication bias in {JASP} & {R} -- selection models, {PET-PEESE}, and robust {B}ayesian meta-analysis}, + author = {Barto{\v{s}}, Franti{\v{s}}ek and Maier, Maximilian and Quintana, Daniel S and Wagenmakers, Eric-Jan}, + year = {2021}, journal = {PsyArXiv}, note = {preprint at \url{https://doi.org/10.31234/osf.io/75bqn}}, - url = {https://doi.org/10.31234/osf.io/75bqn} + doi = {10.31234/osf.io/75bqn} } @misc{jasp14, @@ -83,7 +87,7 @@ @article{bem2011feeling pages = {407}, year = {2011}, publisher = {American Psychological Association}, - url = {https://doi.org/10.1037/a0021524} + doi = {10.1037/a0021524} } @article{bem2011must, @@ -95,7 +99,7 @@ @article{bem2011must number = {4}, pages = {716}, year = {2011}, - url = {https://doi.org/10.1037/a0024777} + doi = {10.1037/a0024777} } @article{wagenmakers2011why, @@ -106,7 +110,7 @@ @article{wagenmakers2011why volume = {100}, number = {3}, pages = {426--432}, - url = {https://doi.org/10.1037/a0022790} + doi = {10.1037/a0022790} } @article{rouder2011bayes, @@ -118,7 +122,7 @@ @article{rouder2011bayes pages = {682--689}, year = {2011}, publisher = {Springer}, - url = {https://doi.org/10.3758/s13423-011-0088-7} + doi = {10.3758/s13423-011-0088-7} } @misc{schimmack2018my, @@ -136,3 +140,56 @@ @book{borenstein2011introduction year = {2011}, publisher = {John Wiley \& Sons} } + +@article{bartos2021bayesian, + title = {Bayesian model-averaged meta-analysis in medicine}, + author = {Barto{\v{s}}, Franti{\v{s}}ek and Gronau, Quentin F and Timmers, Bram and Otte, Willem M. and Ly, Alexander and Wagenmakers, Eric-Jan}, + journal = {Statistics in Medicine}, + doi = {10.1002/sim.9170}, + year = {2021} +} + +@article{poulsen2006potassium, + title = {Potassium containing toothpastes for dentine hypersensitivity}, + author = {Poulsen, Sven and Errboe, Marie and Mevil, Yamila Lescay and Glenny, Anne-Marie}, + journal = {Cochrane Database of Systematic Reviews}, + number = {3}, + year = {2006}, + doi = {10.1002/14651858.cd001476.pub2}, + publisher = {John Wiley \& Sons, Ltd} +} + +@article{stanley2014meta, + title = {Meta-regression approximations to reduce publication selection bias}, + author = {Stanley, Tom D and Doucouliagos, Hristos}, + journal = {Research Synthesis Methods}, + volume = {5}, + number = {1}, + pages = {60--78}, + year = {2014}, + publisher = {Wiley Online Library}, + doi = {10.1002/jrsm.1095} +} + +@article{stanley2017limitations, + title = {Limitations of {PET-PEESE} and other meta-analysis methods}, + author = {Stanley, Tom D}, + journal = {Social Psychological and Personality Science}, + volume = {8}, + number = {5}, + pages = {581--591}, + year = {2017}, + publisher = {Sage Publications Sage CA: Los Angeles, CA}, + doi = {10.1177/1948550617693062} +} + +@article{vevea1995general, + title = {A general linear model for estimating effect size in the presence of publication bias}, + author = {Vevea, Jack L and Hedges, Larry V}, + year = 1995, + journal = {Psychometrika}, + volume = 60, + number = 3, + pages = {419--435}, + doi = {10.1007/BF02294384} +} diff --git a/man/Anderson2010.Rd b/man/Anderson2010.Rd index d7e38c30..ceb95d26 100644 --- a/man/Anderson2010.Rd +++ b/man/Anderson2010.Rd @@ -3,7 +3,7 @@ \docType{data} \name{Anderson2010} \alias{Anderson2010} -\title{27 experimental studies from from +\title{27 experimental studies from \insertCite{anderson2010violent;textual}{RoBMA} that meet the best practice criteria} \format{ A data.frame with 3 columns and 23 observations. diff --git a/man/Bem2011.Rd b/man/Bem2011.Rd index 2b942c65..9361627a 100644 --- a/man/Bem2011.Rd +++ b/man/Bem2011.Rd @@ -3,7 +3,7 @@ \docType{data} \name{Bem2011} \alias{Bem2011} -\title{9 experimental studies from from +\title{9 experimental studies from \insertCite{bem2011feeling;textual}{RoBMA} as described in \insertCite{bem2011must;textual}{RoBMA}} \format{ diff --git a/man/Poulsen2006.Rd b/man/Poulsen2006.Rd new file mode 100644 index 00000000..46ca5aab --- /dev/null +++ b/man/Poulsen2006.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data.R +\docType{data} +\name{Poulsen2006} +\alias{Poulsen2006} +\title{5 studies with a tactile outcome assessment from +\insertCite{poulsen2006potassium;textual}{RoBMA} of the effect of potassium-containing toothpaste +on dentine hypersensitivity} +\format{ +A data.frame with 3 columns and 5 observations. +} +\usage{ +Poulsen2006 +} +\value{ +a data.frame. +} +\description{ +The data set contains Cohen's d effect sizes, standard errors, +and labels for 5 studies assessing the tactile outcome from a meta-analysis of +the effect of potassium-containing toothpaste on dentine hypersensitivity +\insertCite{poulsen2006potassium}{RoBMA} which was used as an example in +\insertCite{bartos2021bayesian;textual}{RoBMA}. +} +\references{ +\insertAllCited{} +} +\keyword{datasets} diff --git a/man/RoBMA-package.Rd b/man/RoBMA-package.Rd index d31e18ea..ed70c62e 100644 --- a/man/RoBMA-package.Rd +++ b/man/RoBMA-package.Rd @@ -19,8 +19,9 @@ inclusion Bayes factors. methodology. More details regarding customization of the model ensembles are provided in the -\href{../doc/ReproducingBMA.html}{\bold{Reproducing BMA}} and -\href{../doc/CustomEnsembles.html}{\bold{Fitting custom meta-analytic ensembles}} +\href{../doc/ReproducingBMA.html}{\bold{Reproducing BMA}}, +\href{../doc/MedicineBMA.html}{\bold{BMA in Medicine}}, and +\href{../doc/CustomEnsembles.html}{\bold{Fitting Custom Meta-Analytic Ensembles}} vignettes. Please, use the "Issues" section in the GitHub repository to ask any further questions. } diff --git a/man/prior_informed.Rd b/man/prior_informed.Rd new file mode 100644 index 00000000..70d832e4 --- /dev/null +++ b/man/prior_informed.Rd @@ -0,0 +1,69 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/priors.R +\name{prior_informed} +\alias{prior_informed} +\title{Creates an informed prior distribution based on research} +\usage{ +prior_informed(name, parameter = NULL, type = "smd") +} +\arguments{ +\item{name}{name of the prior distribution. There are many options based on prior psychological +or medical research. +For psychology, the possible options are +\describe{ +\item{\code{"van Erp"}}{for an informed prior distribution for the heterogeneity parameter tau +of meta-analytic effect size estimates based on standardized mean differences +(van Erp et al. 2017),} +\item{\code{"Oosterwijk"}}{for an informed prior distribution for the effect sizes expected in +social psychology based on prior elicitation with dr. Oosterwijk +(Gronau et al. 2017).} +} +For medicine, the possible options are based on BartoÅ¡ et al. (2021) +who developed empirical prior distributions for the effect size and heterogeneity parameters of the +continuous standardized outcomes based on the Cochrane database of systematic reviews. +Use \code{"Cochrane"} for a prior distribution based on the whole database or call +\code{print(prior_informed_medicine_names)} to inspect the names of +all 46 subfields and set the appropriate \code{parameter} and \code{type}.} + +\item{parameter}{parameter name describing what prior distribution is supposed to be produced in cases +where the \code{name} corresponds to multiple prior distributions. Relevant only for the empirical medical +prior distributions.} + +\item{type}{prior type describing what prior distribution is supposed to be produced in cases +where the \code{name} and \code{parameter} correspond to multiple prior distributions. Relevant only for +the empirical medical prior distributions.} +} +\value{ +\code{prior_informed} returns an object of class 'prior'. +} +\description{ +\code{prior_informed} creates an informed prior distribution based on past +research. The prior can be visualized by the \code{plot} function. +} +\details{ +Further details can be found in \insertCite{erp2017estimates;textual}{RoBMA}, +\insertCite{gronau2017bayesian;textual}{RoBMA}, and +\insertCite{bartos2021bayesian;textual}{RoBMA}. +} +\examples{ +# prior distribution representing expected effect sizes in social psychology +# based on prior elicitation with dr. Oosterwijk +p1 <- prior_informed("Oosterwijk") + +# the prior distribution can be visualized using the plot function +# (see ?plot.prior for all options) +plot(p1) + + +# empirical prior distribution for the standardized mean differences from the oral health +# medical subfield based on meta-analytic effect size estimates from the +# Cochrane database of systematic reviews +p2 <- prior_informed("Oral Health", parameter = "effect", type = "smd") +print(p2) +} +\references{ +\insertAllCited{} +} +\seealso{ +\code{\link[BayesTools:prior]{prior()}}, \link[BayesTools]{prior_informed_medicine_names} +} diff --git a/models/MedicineBMA/fit_BMA.RDS b/models/MedicineBMA/fit_BMA.RDS new file mode 100644 index 00000000..e1834c72 Binary files /dev/null and b/models/MedicineBMA/fit_BMA.RDS differ diff --git a/models/MedicineBMA/fit_BMAb.RDS b/models/MedicineBMA/fit_BMAb.RDS new file mode 100644 index 00000000..13051c00 Binary files /dev/null and b/models/MedicineBMA/fit_BMAb.RDS differ diff --git a/models/MedicineBMA/fit_RoBMA.RDS b/models/MedicineBMA/fit_RoBMA.RDS new file mode 100644 index 00000000..fb74b5c1 Binary files /dev/null and b/models/MedicineBMA/fit_RoBMA.RDS differ diff --git a/vignettes/CustomEnsembles.Rmd b/vignettes/CustomEnsembles.Rmd index b650f2dc..bc0be460 100644 --- a/vignettes/CustomEnsembles.Rmd +++ b/vignettes/CustomEnsembles.Rmd @@ -1,5 +1,5 @@ --- -title: "Fitting custom meta-analytic ensembles" +title: "Fitting Custom Meta-Analytic Ensembles" author: "František Bartoš" date: "`r Sys.Date()`" output: @@ -8,7 +8,7 @@ output: bibliography: ../inst/REFERENCES.bib csl: ../inst/apa.csl vignette: > - %\VignetteIndexEntry{Fitting custom meta-analytic ensembles} + %\VignetteIndexEntry{Fitting Custom Meta-Analytic Ensembles} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown_notangle} @@ -29,11 +29,11 @@ if(.Platform$OS.type == "windows"){ } ``` -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 [@bartos2021no]$^1$. In [another vignette](ReproducingBMA.html), we illustrated that RoBMA can be also utilized for reproducing Bayesian Model-Averaged Meta-Analysis (BMA) [@gronau2017bayesian]. 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 @bartos2020adjusting for a tutorial paper on customizing the model ensemble with JASP). +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 [@bartos2021no]$^1$. In [another vignette](ReproducingBMA.html), we illustrated that RoBMA can be also utilized for reproducing Bayesian Model-Averaged Meta-Analysis (BMA) [@gronau2017bayesian; @gronau2020primer; @bartos2021bayesian]. 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 @bartos2020adjusting 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 completely discussion with the relevant sources in the Example section of @bartos2021no. -### The dataset +### The Dataset To illustrate the custom model building procedure, we use data from the infamous @bem2011feeling "Feeling the future" pre-cognition study. We use coding of the results as summarized by Bem in one of his later replies [@bem2011must]. @@ -44,7 +44,7 @@ data("Bem2011", package = "RoBMA") Bem2011 ``` -### The custom ensemble +### The Custom Ensemble We consider the following scenarios as plausible explanations for the data, and decide to include only those models into the meta-analytic ensemble: @@ -154,7 +154,7 @@ fit <- readRDS(file = "../models/CustomEnsembles/Bem_update2.RDS") summary(fit, type = "models") ``` -### Using the fitted ensemble +### Using the Fitted Ensemble 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 pre-cognition effect, $\text{BF}_{10} = 0.584$ -> $\text{BF}_{01} = 1.71$, moderate evidence for the absence of heterogeneity, $\text{BF}_{\text{rf}} = 0.132$ -> $\text{BF}_{\text{fr}} = 7.58$, and moderate evidence for the presence of the publication bias, $\text{BF}_{\text{pb}} = 3.21$. diff --git a/vignettes/MedicineBMA.Rmd b/vignettes/MedicineBMA.Rmd new file mode 100644 index 00000000..578b4cfe --- /dev/null +++ b/vignettes/MedicineBMA.Rmd @@ -0,0 +1,167 @@ +--- +title: "Informed Bayesian Model-Averaged Meta-Analysis in Medicine" +author: "František Bartoš" +date: "`r Sys.Date()`" +output: + rmarkdown::html_vignette: + self_contained: yes +bibliography: ../inst/REFERENCES.bib +csl: ../inst/apa.csl +vignette: > + %\VignetteIndexEntry{Informed Bayesian Model-Averaged Meta-Analysis in Medicine} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} + %\VignetteEngine{knitr::rmarkdown_notangle} +--- + +```{r setup, include = FALSE} +is_check <- ("CheckExEnv" %in% search()) || + any(c("_R_CHECK_TIMINGS_", "_R_CHECK_LICENSE_") %in% names(Sys.getenv())) || + !file.exists("../models/ReproducingBMA/BMA_PowerPoseTest.RDS") +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + eval = !is_check, + dev = "png") +if(.Platform$OS.type == "windows"){ + knitr::opts_chunk$set(dev.args = list(type = "cairo")) +} +``` +```{r include = FALSE} +library(RoBMA) +# we pre-load the RoBMA models, the fitting time is around 2-5 minutes +fit_BMA <- readRDS(file = "../models/MedicineBMA/fit_BMA.RDS") +fit_BMAb <- readRDS(file = "../models/MedicineBMA/fit_BMAb.RDS") +fit_RoBMA <- readRDS(file = "../models/MedicineBMA/fit_RoBMA.RDS") +``` + + +Bayesian model-averaged meta-analysis allows researchers to seamlessly incorporate available prior information into the analysis [@gronau2017bayesian; @gronau2020primer; @bartos2021bayesian]. This vignette illustrates how to do this with an example from @bartos2021bayesian, 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 [@bartos2021no; @maier2020robust]. + + +### Reproducing Informed Bayesian Model-Averaged Meta-Analysis (BMA) + +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 @bartos2021bayesian, who reanalyzed five studies with a tactile outcome assessment that were subjected to a meta-analysis by @poulsen2006potassium. + +We load the dentine hypersensitivity data included in the package. + +```{r} +library(RoBMA) + +data("Poulsen2006", package = "RoBMA") +Poulsen2006 +``` + +To reproduce the analysis from the example, we need to set informed empirical prior distributions for the effect sizes ($\mu$) and heterogeneity ($\tau$) parameters that @bartos2021bayesian obtained from the Cochrane database of systematic reviews. We can either set them manually, +``` r +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 $\delta \sim T(0,0.51,5)$ and $\tau \sim InvGamma(1.79,0.28)$ informed prior distributions for the "oral health" subfield and removing the publication bias adjustment models by setting `priors_bias = NULL`$^1$. + +Alternatively, we can utilize the `informed_prior` function that prepares informed prior distributions for the individual medical subfields automatically. +``` r +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 (only standardized mean differences `smd` are available at the time of writing the vignette; 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$^2$. + +```{r} +summary(fit_BMA, conditional = TRUE) +``` + +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/Heterogeneity/Publication bias). We can see that there are no models correcting for publication bias (we disabled them by setting `priors_bias = NULL`). Furthermore, 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 @bartos2021bayesian package, $\text{BF}_{10} = 218.53$ and $\text{BF}_{\text{rf}} = 3.52$ (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 @bartos2021bayesian, $\delta = 1.082$, 95% CI [0.686,1.412], and a heterogeneity estimate (not reported previously). + + +### Visualizating the Results + +The `RoBMA` package provides extensive options for visualizing the results. Here, we visualize the prior (grey) and posterior (black) distribution for the mean parameter. + +```{r fig_mu_BMA, dpi = 300, fig.width = 6, fig.height = 4, out.width = "75%", fig.align = "center"} +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 .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, + +```{r fig_mu_BMA_cond, dpi = 300, fig.width = 6, fig.height = 4, out.width = "75%", fig.align = "center"} +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, $H_0^{\text{f}}$, Model 2 corresponds to the random effect model assuming the absence of the effect, $H_0^{\text{r}}$, Model 3 corresponds to the fixed effect model assuming the presence of the effect, $H_3^{\text{f}}$, and Model 4 corresponds to the random effect model assuming the presence of the effect, $H_1^{\text{r}}$,). 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. + +```{r fig_models, dpi = 300, fig.height = 4.5, fig.width = 7, out.width = '75%', fig.align = "center"} +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, $P(H_1^{\text{r}}) = 0.77$) + +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. + +```{r fig_forest, dpi = 300, fig.height = 4.5, fig.width = 7, out.width = '75%', fig.align = "center"} +forest(fit_BMA, conditional = TRUE) +``` + +For more options provided by the plotting function, see its documentation using `?plot.RoBMA()`, `?plot_models()`, and `?forest()`. + + +### Adjusting for Publication Bias with Robust Bayesian Meta-Analysis + +Finally, we illustrate how to adjust our informed BMA for publication bias with robust Bayesian meta-analysis [@bartos2021no, @maier2020robust]. 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 $p$-values [@vevea1995general] or by specifying a publication bias adjustment method correcting for the relationship between effect sizes and standard errors -- PET-PEESE [@stanley2017limitations; @stanley2014meta]. See @bartos2020adjusting 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 $z$). + +``` r +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) +``` +```{r} +summary(fit_BMAb, conditional = TRUE) +``` + +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 $H_1^f$. The increase in the posterior model probability of $H_1^f$ is, in our case, the fact that it predicted the data slightly better after removing the correlation between effect sizes and standard errors (a result of using Cohen's $d$ for model fitting). 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`, which allows us to obtain the default 36 models ensemble called RoBMA-PSMA [@bartos2021no]. + +``` r +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) +``` +```{r} +summary(fit_RoBMA, conditional = TRUE) +``` + +We notice that 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 $\text{BF}_{10} = 6.02$ vs. $\text{BF}_{10} = 347.93$, 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 $\text{BF}_{pb} = 2.30$, which can be interpreted as weak evidence in favor of the 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 $\delta = 0.838$, 95% CI [-0.035, 1.297] with a much wider confidence interval, as visualized in the prior and posterior conditional effect size estimate plot. + +```{r fig_mu_RoBMA_cond, dpi = 300, fig.width = 6, fig.height = 4, out.width = "75%", fig.align = "center"} +plot(fit_RoBMA, parameter = "mu", prior = TRUE, conditional = TRUE) +``` + + +### Footnotes + +$^1$ 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 $z$ for the fitting purposes and the `seed = 1` and `parallel = TRUE` options grant us exact reproducibility of the results and parallelization of the fitting process. + +$^2$ 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`. + +### References diff --git a/vignettes/ReproducingBMA.Rmd b/vignettes/ReproducingBMA.Rmd index 61b63e5e..d04ac87b 100644 --- a/vignettes/ReproducingBMA.Rmd +++ b/vignettes/ReproducingBMA.Rmd @@ -1,5 +1,5 @@ --- -title: "Reproducing BMA" +title: "Reproducing Bayesian Model-Averaged Meta-Analysis" author: "František Bartoš" date: "`r Sys.Date()`" output: @@ -8,7 +8,7 @@ output: bibliography: ../inst/REFERENCES.bib csl: ../inst/apa.csl vignette: > - %\VignetteIndexEntry{Reproducing BMA} + %\VignetteIndexEntry{Reproducing Bayesian Model-Averaged Meta-Analysis} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown_notangle} @@ -36,7 +36,7 @@ fit_RoBMA_test <- readRDS(file = "../models/ReproducingBMA/PowerPoseTest.RDS") fit_RoBMA_est <- readRDS(file = "../models/ReproducingBMA/PowerPoseEst.RDS") ``` -By default, the package estimates an ensemble of 36 meta-analytic models and provides functions for convenient manipulation with the fitted object. However, it has been built in a way that 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 [@gronau2017bayesian, @gronau2020primer]. See [this vignette](CustomEnsembles.html) if you are interested in building more customized ensembles or @bartos2020adjusting for a tutorial on fitting (custom) models in JASP. +By default, the package estimates an ensemble of 36 meta-analytic models and provides functions for convenient manipulation with the fitted object. However, it has been built in a way that 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 [@gronau2017bayesian; @gronau2020primer; @bartos2021bayesian]. See [this vignette](CustomEnsembles.html) if you are interested in building more customized ensembles or @bartos2020adjusting for a tutorial on fitting (custom) models in JASP. ### Reproducing Bayesian Model-Averaged Meta-Analysis (BMA) @@ -109,7 +109,7 @@ The output from the `summary.RoBMA()` function has 2 parts. The first one under 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 $\mu = 0.22$, 95% CI [0.10, 0.35] that also mirrors the one obtained from `metaBMA` package. -### Visualizating the results +### Visualizating the Results RoBMA provides extensive options for visualizing the results. Here, we visualize the prior (grey) and posterior (black) distribution for the mean parameter.