Add Bayesian meta-analysis model

[LukasWallrich]

Can take inspiration - but not code directly from https://github.com/neuromadlab/taVNSHRVmeta

[LukasWallrich]

RoBMA estimate - https://cran.r-project.org/web/packages/RoBMA/readme/README.html#:~:text=This%20package%20estimates%20an%20ensemble,model%20averaging%20to%20combine%20them.

[LukasWallrich]

Thanks @LukasRoeseler for providing some initial code to integrate this. However, at present, the calculation is too slow to include it in the model ensemble ... so parking this for now. It might be substantially faster with weakly informative rather than uniform priors ... but I am too unfamiliar with bayesian analysis to propose that.

# Bayesian Meta-Analyses --------------------------------------------------

## Installing packages
# install.packages("bayesmeta")
# install.packages("psymetadata")
# install.packages("RoBMA")
# devtools::install_github("fbartos/RoBMA")

library(psymetadata)
library(RoBMA)
library(bayesmeta)

## NOTES
# JAGS is required: https://sourceforge.net/projects/mcmc-jags/
# both analyses currently work for simple structures, only (1 effect per study)

ds <- psymetadata::aksayli2019
ds <- ds[1:50, ]

# Fitting the standard RoBMA model ----------------------------------------

t0 <- Sys.time()
fit <- RoBMA::RoBMA(d = ds$yi
                    , se = sqrt(ds$vi)
                    , study_names = ds$test
                    , seed = 1
                    , effect_direction = "positive" # this is the default direction and would need to be adapted similarly to the effect size sign p-curve
                    , silent = FALSE) # Effect size is Hedges's g here
t1 <- Sys.time()
t1-t0 # I stopped this after it had fited 14 models, which took 1.12 hours on my machine



# Fitting standard bayesian model -----------------------------------------
t0 <- Sys.time()
bm <- bayesmeta::bayesmeta(y = ds$yi
                     , sigma = sqrt(ds$vi)
                     , label = ds$test
                     , tau.prior = "uniform", mu.prior = c("mean" = NA, "sd" = NA) # options to modify priors (more available)
                     )
t1 <- Sys.time()
t1-t0 # this took >11 minutes (cancelled for whole dataset with k = 637) and 1.23 for only 50 effects

# Values for plotting
# model, es, lcl, ucl, k
bm
forestplot::forestplot(bm)
print(bm)

## Code to facilitate integration in the App
estimates_explo_agg <- data.frame("Model" = as.character(), "g" = as.numeric(), "LCL" = as.numeric(), "UCL" = as.numeric(), "k" = as.numeric()
                                  , stringsAsFactors = FALSE)

estimates_explo_agg[nrow(estimates_explo_agg)+1,] <- c("Bayesian Random-Effects"
                                                       , bm$summary["mean", "mu"]
                                                       , bm$summary["95% lower", "mu"]
                                                       , bm$summary["95% upper", "mu"]
                                                       , bm$k)