Multilevel (#17)

- multivariate (weighted) normal distributions (both R & C++) - removing macOS from actions (new C++ version does not compile the package)
FBartos · Mar 2, 2022 · abc4e74 · abc4e74
1 parent 89b2486
commit abc4e74
Show file tree

Hide file tree

Showing 98 changed files with 6,071 additions and 176 deletions.
diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml
@@ -25,10 +25,10 @@ jobs:
       matrix:
         config:
           - {os: windows-latest, r: 'release'}
-          - {os: macOS-latest,   r: 'release'}
           - {os: ubuntu-20.04,   r: 'release', rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest"}
           - {os: ubuntu-20.04,   r: 'devel',   rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest", http-user-agent: "R/4.1.0 (ubuntu-20.04) R (4.1.0 x86_64-pc-linux-gnu x86_64 linux-gnu) on GitHub Actions" }
-
+        #  - {os: macOS-11.6.2,   r: 'release'} # the release C++ version does not compile properly
+
     env:
       R_REMOTES_NO_ERRORS_FROM_WARNINGS: true
       RSPM: ${{ matrix.config.rspm }}
@@ -104,6 +104,7 @@ jobs:
         env:
           _R_CHECK_CRAN_INCOMING_REMOTE_: false
           JAGS_ROOT: "/c/progra~1/JAGS/JAGS-4.2.0"
+          JAGS_MAJOR_VERSION: 4
         run: |
           options(crayon.enabled = TRUE)
           rcmdcheck::rcmdcheck(

diff --git a/.github/workflows/R-CMD-tests.yaml b/.github/workflows/R-CMD-tests.yaml
@@ -25,8 +25,8 @@ jobs:
       matrix:
         config:
           - {os: windows-latest, r: 'release'}
-          - {os: macOS-latest,   r: 'release'}
           - {os: ubuntu-20.04,   r: 'release', rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest"}
+        #  - {os: macOS-11.6.2,   r: 'release'} # the release C++ version does not compile properly
 
     env:
       R_REMOTES_NO_ERRORS_FROM_WARNINGS: true
@@ -87,7 +87,7 @@ jobs:
           remotes::install_cran("rcmdcheck")
           remotes::install_version("BayesTools", "0.1.3")
           install.packages("devtools")
-          install.packages("testhat")
+          install.packages("testthat")
           install.packages("vdiffr")
           devtools::install()
           library(vdiffr)
@@ -100,12 +100,13 @@ jobs:
         env:
           _R_CHECK_CRAN_INCOMING_REMOTE_: false
           JAGS_ROOT: "/c/progra~1/JAGS/JAGS-4.2.0"
+          JAGS_MAJOR_VERSION: 4
         run: |
           remotes::install_deps(dependencies = TRUE)
           remotes::install_cran("rcmdcheck")
           remotes::install_version("BayesTools", "0.1.3")
           install.packages("devtools")
-          install.packages("testhat")
+          install.packages("testthat")
           install.packages("vdiffr")
           devtools::install()
           library(vdiffr)

diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml
@@ -10,9 +10,19 @@ name: pkgdown
 
 jobs:
   pkgdown:
-    runs-on: macOS-latest
+    runs-on: ${{ matrix.config.os }}
+
+    name: ${{ matrix.config.os }} (${{ matrix.config.r }})
+
+    strategy:
+      fail-fast: false
+      matrix:
+        config:
+          - {os: ubuntu-20.04,   r: 'release', rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest"}
+
     env:
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+
     steps:
       - uses: actions/checkout@v2
 

diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml
@@ -14,9 +14,19 @@ name: test-coverage
 
 jobs:
   test-coverage:
-    runs-on: macOS-latest
+    runs-on: ${{ matrix.config.os }}
+
+    name: ${{ matrix.config.os }} (${{ matrix.config.r }})
+
+    strategy:
+      fail-fast: false
+      matrix:
+        config:
+          - {os: ubuntu-20.04,   r: 'release', rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest"}
+
     env:
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+
     steps:
       - uses: actions/checkout@v2
 

diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: RoBMA
 Title: Robust Bayesian Meta-Analyses
-Version: 2.1.2
+Version: 2.2.0
 Maintainer: František Bartoš <[email protected]>
 Authors@R: c( 
     person("František", "Bartoš",     role = c("aut", "cre"),
@@ -49,6 +49,7 @@ Imports:
     stats,
     graphics,
     extraDistr,
+    mvtnorm,
     scales,
     callr,
     Rdpack,
@@ -62,5 +63,7 @@ Suggests:
     knitr,
     rmarkdown,
     covr
+LinkingTo:
+    mvtnorm
 RdMacros: Rdpack
 VignetteBuilder: knitr
diff --git a/NEWS.md b/NEWS.md
@@ -1,3 +1,7 @@
+## version 2.2
+### Features
+- three-level meta-analysis by specifying `study_ids` argument in `RoBMA`. However, note that this is (1) an experimental feature and (2) the computational expense of fitting selection models with clustering is extreme. As of now, it is almost impossible to have more than 2-3 estimates clustered within a single study).
+
 ## version 2.1.2
 ### Fixes
 - adding Windows ucrt patch (thanks to Tomas Kalibera)

diff --git a/R/check-input-and-settings.R b/R/check-input-and-settings.R
@@ -29,11 +29,13 @@ check_setup <- function(
   priors_effect_null         = prior(distribution = "point", parameters = list(location = 0)),
   priors_heterogeneity_null  = prior(distribution = "point", parameters = list(location = 0)),
   priors_bias_null           = prior_none(),
+  priors_rho                 = prior("beta", parameters = list(alpha = 1, beta = 1)),
+  priors_rho_null            = NULL,
   models = FALSE, silent = FALSE){
 
   object <- list()
-  object$priors   <- .check_and_list_priors(tolower(model_type), priors_effect_null, priors_effect, priors_heterogeneity_null, priors_heterogeneity, priors_bias_null, priors_bias, object$add_info[["prior_scale"]])
-  object$models   <- .make_models(object[["priors"]])
+  object$priors   <- .check_and_list_priors(tolower(model_type), priors_effect_null, priors_effect, priors_heterogeneity_null, priors_heterogeneity, priors_bias_null, priors_bias, priors_rho, priors_rho_null, object$add_info[["prior_scale"]])
+  object$models   <- .make_models(object[["priors"]], multivariate = FALSE)
 
 
   ### model types overview

diff --git a/R/check-priors-and-models.R b/R/check-priors-and-models.R
@@ -1,5 +1,5 @@
 ### functions for creating model objects
-.check_and_list_priors            <- function(model_type, priors_effect_null, priors_effect, priors_heterogeneity_null, priors_heterogeneity, priors_bias_null, priors_bias, prior_scale){
+.check_and_list_priors            <- function(model_type, priors_effect_null, priors_effect, priors_heterogeneity_null, priors_heterogeneity, priors_bias_null, priors_bias, priors_rho_null, priors_rho, prior_scale){
 
   if(!is.null(model_type) & length(model_type == 1)){
     # precanned models
@@ -13,32 +13,38 @@
         prior_weightfunction(distribution = "one.sided", parameters = list(alpha = c(1, 1, 1),    steps = c(0.025, 0.05)),      prior_weights = 1/12),
         prior_weightfunction(distribution = "one.sided", parameters = list(alpha = c(1, 1, 1),    steps = c(0.05, 0.5)),        prior_weights = 1/12),
         prior_weightfunction(distribution = "one.sided", parameters = list(alpha = c(1, 1, 1, 1), steps = c(0.025, 0.05, 0.5)), prior_weights = 1/12),
-        prior_PET(distribution = "Cauchy", parameters = list(0,1), truncation = list(0, Inf),             prior_weights = 1/4),
-        prior_PEESE(distribution = "Cauchy", parameters = list(0,5), truncation = list(0, Inf),             prior_weights = 1/4)
+        prior_PET(distribution = "Cauchy", parameters = list(0,1), truncation = list(0, Inf),   prior_weights = 1/4),
+        prior_PEESE(distribution = "Cauchy", parameters = list(0,5), truncation = list(0, Inf), prior_weights = 1/4)
       )
+      priors_rho                <- NULL
       priors_effect_null        <- prior(distribution = "point", parameters = list(location = 0))
       priors_heterogeneity_null <- prior(distribution = "point", parameters = list(location = 0))
       priors_bias_null          <- prior_none()
+      priors_rho_null           <- NULL
     }else if(model_type == "pp"){
       priors_effect          <- prior(distribution = "normal",    parameters = list(mean = 0, sd = 1))
       priors_heterogeneity   <- prior(distribution = "invgamma",  parameters = list(shape = 1, scale = .15))
       priors_bias            <- list(
-        prior_PET(distribution = "Cauchy", parameters = list(0,1), truncation = list(0, Inf),  prior_weights = 1/2),
+        prior_PET(distribution = "Cauchy", parameters = list(0,1), truncation = list(0, Inf),    prior_weights = 1/2),
         prior_PEESE(distribution = "Cauchy", parameters = list(0,5), truncation = list(0, Inf),  prior_weights = 1/2)
       )
+      priors_rho                 <- NULL
       priors_effect_null         <- prior(distribution = "point", parameters = list(location = 0))
       priors_heterogeneity_null  <- prior(distribution = "point", parameters = list(location = 0))
       priors_bias_null           <- prior_none()
+      priors_rho_null            <- NULL
     }else if(model_type == "2w"){
       priors_effect              <- prior(distribution = "normal",    parameters = list(mean = 0, sd = 1))
       priors_heterogeneity       <- prior(distribution = "invgamma",  parameters = list(shape = 1, scale = .15))
       priors_bias                <- list(
         prior_weightfunction(distribution = "two.sided", parameters = list(alpha = c(1, 1),       steps = c(0.05)),        prior_weights = 1/2),
         prior_weightfunction(distribution = "two.sided", parameters = list(alpha = c(1, 1, 1),    steps = c(0.05, 0.10)),  prior_weights = 1/2)
       )
+      priors_rho                 <- NULL
       priors_effect_null         <- prior(distribution = "point", parameters = list(location = 0))
       priors_heterogeneity_null  <- prior(distribution = "point", parameters = list(location = 0))
       priors_bias_null           <- prior_none()
+      priors_rho_null            <- NULL
     }else{
       stop("Unknown 'model_type'.")
     }
@@ -48,13 +54,14 @@
   priors$effect         <- .check_and_list_component_priors(priors_effect_null,          priors_effect,         "effect")
   priors$heterogeneity  <- .check_and_list_component_priors(priors_heterogeneity_null,   priors_heterogeneity,  "heterogeneity")
   priors$bias           <- .check_and_list_component_priors(priors_bias_null,            priors_bias,           "bias")
+  priors$rho            <- .check_and_list_component_priors(priors_rho_null,             priors_rho,            "rho")
 
   return(priors)
 }
 .check_and_list_component_priors  <- function(priors_null, priors_alt, component){
 
   # check that at least one prior is specified (either null or alternative)
-  if(is.null(priors_null) & is.null(priors_alt))
+  if(component != "rho" && (is.null(priors_null) & is.null(priors_alt)))
     stop(paste0("At least one prior needs to be specified for the ", component," parameter (either null or alternative)."))
 
   # create an empty list if user didn't specified priors
@@ -132,28 +139,45 @@
       if(!(is.prior.PET(priors[[p]]) | is.prior.PEESE(priors[[p]]) | is.prior.weightfunction(priors[[p]]) | is.prior.none(priors[[p]])))
         stop(paste0("'", print(priors[[p]], silent = TRUE),"' prior distribution is not supported for the bias component."))
     }
+  }else if(component == "rho"){
+
+    for(p in seq_along(priors)){
+
+      # check for allowed priors
+      if(!(priors[[p]][["distribution"]] == "beta"))
+        stop(paste0("'", print(priors[[p]], silent = TRUE),"' prior distribution is not supported for the rho component."))
+    }
   }
 
   return(priors)
 }
-.make_models  <- function(priors){
+.make_models  <- function(priors, multivariate){
 
   # create models according to the set priors
   models <- NULL
   for(effect in priors[["effect"]]){
     for(heterogeneity in priors[["heterogeneity"]]){
       for(bias in priors[["bias"]]){
-        models <- c(
-          models,
-          list(.make_model(effect, heterogeneity, bias, effect[["prior_weights"]] * heterogeneity[["prior_weights"]] * bias[["prior_weights"]]))
-        )
+        if(!is.null(priors[["rho"]]) && multivariate){
+          for(rho in priors[["rho"]]){
+            models <- c(
+              models,
+              list(.make_model(effect, heterogeneity, bias, rho, effect[["prior_weights"]] * heterogeneity[["prior_weights"]] * bias[["prior_weights"]] * rho[["prior_weights"]]))
+            )
+          }
+        }else{
+          models <- c(
+            models,
+            list(.make_model(effect, heterogeneity, bias, NULL, effect[["prior_weights"]] * heterogeneity[["prior_weights"]] * bias[["prior_weights"]]))
+          )
+        }
       }
     }
   }
 
   return(models)
 }
-.make_model   <- function(prior_effect, prior_heterogeneity, prior_bias, prior_weights){
+.make_model   <- function(prior_effect, prior_heterogeneity, prior_bias, prior_rho, prior_weights){
 
   priors <- list()
 
@@ -166,13 +190,18 @@
   }else if(is.prior.weightfunction(prior_bias)){
     priors$omega  <- prior_bias
   }
+  # add 3 level structure only if there is heterogeneity
+  if(!(prior_heterogeneity[["distribution"]] == "point" && prior_heterogeneity$parameters[["location"]] == 0) && !is.null(prior_rho)){
+    priors$rho    <- prior_rho
+  }
 
   model <- list(
     priors         = priors,
     prior_weights     = prior_weights,
     prior_weights_set = prior_weights
   )
   class(model) <- "RoBMA.model"
+  attr(model, "multivariate") <- !is.null(priors$rho)
 
   return(model)
 }