BayesTools 0.2.14 update (#3)

FBartos · May 30, 2023 · da46631 · da46631
1 parent 4c77b00
commit da46631
Show file tree

Hide file tree

Showing 29 changed files with 277 additions and 162 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: RoBSA
 Type: Package
 Title: Robust Bayesian Survival Analysis
-Version: 1.0.1
+Version: 1.0.2
 Maintainer: František Bartoš <f.bartos96@gmail.com>
 Authors@R: c( 
     person("František", "Bartoš", role = c("aut", "cre"),
@@ -19,7 +19,7 @@ Description: A framework for estimating ensembles of parametric survival models
     a model ensemble, weights the posterior parameter distributions based on 
     posterior model probabilities and uses Bayes factors to test for the 
     presence or absence of the individual predictors or preference for a 
-    parametric family (Bartoš, Aust & Haaf, 2021, <doi:10.48550/arXiv.2112.08311>).
+    parametric family (Bartoš, Aust & Haaf, 2022, <doi:10.1186/s12874-022-01676-9>).
     The user can define a wide range of informative priors for all parameters 
     of interest. The package provides convenient functions for summary, visualizations, 
     fit diagnostics, and prior distribution calibration.
@@ -34,7 +34,7 @@ SystemRequirements: JAGS >= 4.3.1 (https://mcmc-jags.sourceforge.io/)
 Depends:
     R (>= 4.0.0)
 Imports:
-    BayesTools (>= 0.2.10),
+    BayesTools (>= 0.2.14),
     survival,
     rjags,
     runjags,

diff --git a/NAMESPACE b/NAMESPACE
@@ -13,6 +13,7 @@ export(calibrate_meta_analytic)
 export(calibrate_quartiles)
 export(check_RoBSA)
 export(check_setup)
+export(contr.meandif)
 export(contr.orthonormal)
 export(diagnostics)
 export(diagnostics_autocorrelation)

diff --git a/NEWS.md b/NEWS.md
@@ -1,3 +1,10 @@
+## version 1.0.2
+### features
+- adding meandif contrasts to the `prior_factor()` function 
+
+### changes
+- update to BayesTools 0.2.14
+
 ## version 1.0.1
 ### changes
 - update to R 4.2 and JAGS 4.3.1

diff --git a/R/check-input-and-settings.R b/R/check-input-and-settings.R
@@ -481,7 +481,7 @@ get_default_prior_beta_alt  <- function(){
 }
 #' @rdname default_prior
 get_default_prior_factor_null  <- function(){
-  prior("spike", parameters = list("location" = 0))
+  prior_factor("spike", parameters = list(location = 0), contrast = "treatment")
 }
 #' @rdname default_prior
 get_default_prior_factor_alt   <- function(){
@@ -508,6 +508,10 @@ get_default_prior_aux       <- function(){
   )
 }
 
+.fix_default_prior_factor_null <- function(old_prior, contrast){
+  prior_factor(distribution = old_prior[["distribution"]], parameters = old_prior[["parameters"]], truncation = old_prior[["truncation"]], prior_weights = old_prior[["prior_weights"]], contrast = contrast)
+}
+
 # some functions for the JASP implementation
 .RoBSA_collect_dots      <- function(...){
 

diff --git a/R/check-priors-and-models.R b/R/check-priors-and-models.R
@@ -124,13 +124,25 @@
           null  = if(predictors_type[to_test[i]] == "factor") default_prior_factor_null else default_prior_beta_null,
           alt   = priors[[to_test[i]]]
         )
+        # fix possible missalignment of default and custom prior distribution
+        if(BayesTools::is.prior.treatment(priors[[to_test[i]]][["alt"]]) && !BayesTools::is.prior.treatment(priors[[to_test[i]]][["null"]])){
+          priors[[to_test[i]]][["null"]] <- .fix_default_prior_factor_null(priors[[to_test[i]]][["null"]], contrast = "treatment")
+        }else if(BayesTools::is.prior.orthonormal(priors[[to_test[i]]][["alt"]]) && !BayesTools::is.prior.orthonormal(priors[[to_test[i]]][["null"]])){
+          priors[[to_test[i]]][["null"]] <- .fix_default_prior_factor_null(priors[[to_test[i]]][["null"]], contrast = "orthonormal")
+        }else if(BayesTools::is.prior.meandif(priors[[to_test[i]]][["alt"]]) && !BayesTools::is.prior.meandif(priors[[to_test[i]]][["null"]])){
+          priors[[to_test[i]]][["null"]] <- .fix_default_prior_factor_null(priors[[to_test[i]]][["null"]], contrast = "meandif")
+        }
       }else if(length(priors[[to_test[i]]]) == 2 && all(names(priors[[to_test[i]]]) %in% c("null", "alt"))){
         priors[[to_test[i]]] <- list(
           null  = priors[[to_test[i]]][["null"]],
           alt   = priors[[to_test[i]]][["alt"]]
         )
+        if((BayesTools::is.prior.treatment(priors[[to_test[i]]][["alt"]])   && !BayesTools::is.prior.treatment(priors[[to_test[i]]][["null"]]))   |
+           (BayesTools::is.prior.orthonormal(priors[[to_test[i]]][["alt"]]) && !BayesTools::is.prior.orthonormal(priors[[to_test[i]]][["null"]])) |
+           (BayesTools::is.prior.meandif(priors[[to_test[i]]][["alt"]])     && !BayesTools::is.prior.meandif(priors[[to_test[i]]][["null"]])))
+          stop(paste0("The null and alternative prior distribution for '", to_test[i], "' factor must have matching contrasts."))
       }else{
-        stop(paste0("The predictor '", to_test[i], "' is supposed to be used for testing and the prior distributions are not specified properly"))
+        stop(paste0("The predictor '", to_test[i], "' is supposed to be used for testing and the prior distributions are not specified properly."))
       }
     }
 

diff --git a/R/diagnostics.R b/R/diagnostics.R
@@ -93,7 +93,7 @@ diagnostics <- function(fit, parameter = NULL, type, plot_type = "base", show_mo
       args$plot_type                <- plot_type
       args$lags                     <- lags
       args$transformations          <- NULL
-      args$transform_orthonormal    <- TRUE
+      args$transform_factors    <- TRUE
       args$short_name               <- FALSE
       args$parameter_names          <- FALSE
       args$formula_prefix           <- FALSE

diff --git a/R/fit-and-marglik.R b/R/fit-and-marglik.R
@@ -137,7 +137,7 @@
 
   # add model summaries
   if(has_posterior){
-    fit_summary <- BayesTools::runjags_estimates_table(fit = fit, warnings = warnings, transform_orthonormal = TRUE, formula_prefix = FALSE)
+    fit_summary <- BayesTools::runjags_estimates_table(fit = fit, warnings = warnings, transform_factors = TRUE, formula_prefix = FALSE)
   }else{
     fit_summary <- BayesTools::runjags_estimates_empty_table()
   }

diff --git a/R/priors.R b/R/priors.R
@@ -11,7 +11,9 @@ prior_none <- BayesTools::prior_none
 #' @name prior_factor
 #' @inherit BayesTools::prior_factor
 #' @export
-prior_factor <- BayesTools::prior_factor
+prior_factor <- function(distribution, parameters, truncation = list(lower = -Inf, upper = Inf), prior_weights = 1, contrast = "orthonormal"){
+  BayesTools::prior_factor(distribution = distribution, parameters = parameters, truncation = truncation, prior_weights = prior_weights, contrast = contrast)
+}
 
 #' @name prior_informed
 #' @inherit BayesTools::prior_informed
@@ -53,3 +55,30 @@ prior_informed_medicine_names <- BayesTools::prior_informed_medicine_names
 #'
 #' @export
 contr.orthonormal <- BayesTools::contr.orthonormal
+
+
+#' @title Mean difference contrast matrix
+#'
+#' @description Return a matrix of mean difference contrasts.
+#' This is an adjustment to the \code{contr.orthonormal} that ascertains that the prior
+#' distributions on difference between the gran mean and factor level are identical independent
+#' of the number of factor levels (which does not hold for the orthonormal contrast). Furthermore,
+#' the contrast is re-scaled so the specified prior distribution exactly corresponds to the prior
+#' distribution on difference between each factor level and the grand mean -- this is approximately
+#' twice the scale of \code{contr.orthonormal}.
+#'
+#' @param n a vector of levels for a factor, or the number of levels
+#' @param contrasts logical indicating whether contrasts should be computed
+#'
+#' @examples
+#' contr.meandif(c(1, 2))
+#' contr.meandif(c(1, 2, 3))
+#'
+#' @references
+#' \insertAllCited{}
+#'
+#' @return A matrix with n rows and k columns, with k = n - 1 if \code{contrasts = TRUE} and k = n
+#' if \code{contrasts = FALSE}.
+#'
+#' @export
+contr.meandif <- BayesTools::contr.meandif
diff --git a/R/summary.R b/R/summary.R
@@ -34,7 +34,7 @@ print.RoBSA <- function(x, ...){
 #' @param logBF show log of the BFs. Defaults to \code{FALSE}.
 #' @param BF01 show BF in support of the null hypotheses. Defaults to
 #' \code{FALSE}.
-#' @param transform_orthonormal Whether factors with orthonormal prior
+#' @param transform_factors Whether factors with orthonormal prior
 #' distributions should be transformed to differences from the grand mean. Defaults
 #' to \code{TRUE}.
 #' @param ... additional arguments
@@ -98,7 +98,7 @@ print.RoBSA <- function(x, ...){
 #' @export
 summary.RoBSA       <- function(object, type = "ensemble", conditional = FALSE,
                                 exp = FALSE, parameters = FALSE, probs = c(.025, .975), logBF = FALSE, BF01 = FALSE,
-                                transform_orthonormal = TRUE, short_name = FALSE, remove_spike_0 = FALSE, ...){
+                                transform_factors = TRUE, short_name = FALSE, remove_spike_0 = FALSE, ...){
 
   BayesTools::check_bool(conditional, "conditional")
   BayesTools::check_char(type, "type")
@@ -107,7 +107,7 @@ summary.RoBSA       <- function(object, type = "ensemble", conditional = FALSE,
   BayesTools::check_real(probs, "probs", allow_NULL = TRUE, check_length = 0)
   BayesTools::check_bool(BF01,  "BF01")
   BayesTools::check_bool(logBF, "logBF")
-  BayesTools::check_bool(transform_orthonormal, "transform_orthonormal")
+  BayesTools::check_bool(transform_factors, "transform_factors")
   BayesTools::check_bool(short_name, "short_name")
   BayesTools::check_bool(remove_spike_0, "remove_spike_0")
 
@@ -152,7 +152,7 @@ summary.RoBSA       <- function(object, type = "ensemble", conditional = FALSE,
         probs                 = probs,
         title                 = "Model-averaged estimates:",
         warnings              = .collect_errors_and_warnings(object),
-        transform_orthonormal = transform_orthonormal,
+        transform_factors = transform_factors,
         formula_prefix        = FALSE
       )
       if(exp){
@@ -259,16 +259,29 @@ summary.RoBSA       <- function(object, type = "ensemble", conditional = FALSE,
       components[[object$add_info[["predictors"]][i]]] <- .BayesTools_parameter_name(object$add_info[["predictors"]][i])
     }
 
-    summary <- BayesTools::ensemble_summary_table(
-      models         = object[["models"]],
-      parameters     = components,
-      title          = "Models overview:",
-      footnotes      = NULL,
-      warnings       = .collect_errors_and_warnings(object),
-      short_name     = short_name,
-      remove_spike_0 = remove_spike_0
-    )
-
+    # trick for dealing with empty models
+    if(length(components) == 0){
+      summary <- BayesTools::ensemble_summary_table(
+        models         = object[["models"]],
+        parameters     = "mu_intercept",
+        title          = "Models overview:",
+        footnotes      = NULL,
+        warnings       = .collect_errors_and_warnings(object),
+        short_name     = short_name,
+        remove_spike_0 = remove_spike_0
+      )
+      summary <- BayesTools::remove_column(summary, 2)
+    }else{
+      summary <- BayesTools::ensemble_summary_table(
+        models         = object[["models"]],
+        parameters     = components,
+        title          = "Models overview:",
+        footnotes      = NULL,
+        warnings       = .collect_errors_and_warnings(object),
+        short_name     = short_name,
+        remove_spike_0 = remove_spike_0
+      )
+    }
 
     # add distribution column to the summary
     summary <- BayesTools::add_column(
@@ -303,15 +316,30 @@ summary.RoBSA       <- function(object, type = "ensemble", conditional = FALSE,
       components[[object$add_info[["predictors"]][i]]] <- .BayesTools_parameter_name(object$add_info[["predictors"]][i])
     }
 
-    diagnostics <- BayesTools::ensemble_diagnostics_table(
-      models         = object[["models"]],
-      parameters     = components,
-      title          = "Diagnostics overview:",
-      footnotes      = NULL,
-      warnings       = .collect_errors_and_warnings(object),
-      short_name     = short_name,
-      remove_spike_0 = remove_spike_0
-    )
+
+    # trick for dealing with empty models
+    if(length(components) == 0){
+      diagnostics <- BayesTools::ensemble_diagnostics_table(
+        models         = object[["models"]],
+        parameters     = "mu_intercept",
+        title          = "Diagnostics overview:",
+        footnotes      = NULL,
+        warnings       = .collect_errors_and_warnings(object),
+        short_name     = short_name,
+        remove_spike_0 = remove_spike_0
+      )
+      diagnostics <- BayesTools::remove_column(diagnostics, 2)
+    }else{
+      diagnostics <- BayesTools::ensemble_diagnostics_table(
+        models         = object[["models"]],
+        parameters     = components,
+        title          = "Diagnostics overview:",
+        footnotes      = NULL,
+        warnings       = .collect_errors_and_warnings(object),
+        short_name     = short_name,
+        remove_spike_0 = remove_spike_0
+      )
+    }
 
 
     # add distribution column to the summary

diff --git a/R/utilities.R b/R/utilities.R
@@ -80,8 +80,9 @@ assign("distributions",   .distributions,                             envir = Ro
 
   BayesTools_required <- switch(
     paste0(RoBSA.version, collapse = "."),
-    "1.0.0" = c("0.2.10", "6.6.6"),
-    "1.0.1" = c("0.2.10", "6.6.6"),
+    "1.0.0" = c("0.2.10", "0.2.13"),
+    "1.0.1" = c("0.2.10", "0.2.13"),
+    "1.0.2" = c("0.2.14", "999.999.999"),
     stop("New RoBSA version needs to be defined in '.check_BayesTools' function!")
   )
 

diff --git a/README.Rmd b/README.Rmd
@@ -74,7 +74,7 @@ predicted_experimental <- readRDS(file = "models/README/predicted_experimental.t
 
 This package estimates an ensemble of parametric survival models with different parametric families and uses Bayesian model averaging to combine them. The RoBSA ensemble uses Bayes factors to test for the presence or absence of effects of the individual predictors and evaluates the support for each parametric family. The resulting model-averaged parameter estimates are based on posterior model probabilities. The user can define a wide range of prior distributions for the effect size, intercepts, and auxiliary parameters. The package provides convenient functions for summary, visualizations, and fit diagnostics.
 
-See our pre-print @bartos2021informed (https://arxiv.org/abs/2112.08311) introducing the methodology.
+See @bartos2021informed (https://doi.org/10.1186/s12874-022-01676-9) introducing the methodology.
 
 ## Installation
 

diff --git a/inst/REFERENCES.bib b/inst/REFERENCES.bib
@@ -1,8 +1,7 @@
-@unpublished{bartos2021informed,
+@article{bartos2021informed,
   title     = {Informed {B}ayesian survival analysis},
   author    = {Barto{\v{s}}, Franti{\v{s}}ek and Aust, Frederik and Haaf, Julia M.},
-  year      = {2021},
-  publisher = {arXiv},
-  note      = {preprint at \url{https://arxiv.org/abs/2112.08311}},
-  url       = {https://arxiv.org/abs/2112.08311}
+  year      = {2022},
+  journal   = {BMC Medical Research Methodology},
+  doi       = {10.1186/s12874-022-01676-9}
 }
diff --git a/man/contr.meandif.Rd b/man/contr.meandif.Rd
diff --git a/man/prior_factor.Rd b/man/prior_factor.Rd
diff --git a/man/summary.RoBSA.Rd b/man/summary.RoBSA.Rd
diff --git a/tests/results/fits/fit_1.RDS b/tests/results/fits/fit_1.RDS
diff --git a/tests/results/fits/fit_2.RDS b/tests/results/fits/fit_2.RDS
diff --git a/tests/results/fits/fit_3.RDS b/tests/results/fits/fit_3.RDS
diff --git a/tests/results/fits/fit_4.RDS b/tests/results/fits/fit_4.RDS
diff --git a/tests/results/fits/fit_5.RDS b/tests/results/fits/fit_5.RDS
diff --git a/tests/results/fits/fit_6.RDS b/tests/results/fits/fit_6.RDS
diff --git a/tests/results/fits/fit_7.RDS b/tests/results/fits/fit_7.RDS
diff --git a/tests/results/summary.diagnostics/4.txt b/tests/results/summary.diagnostics/4.txt
@@ -5,11 +5,11 @@
 "Robust Bayesian survival analysis"
 "Diagnostics overview:"
 " Model Distribution            Prior x_bin           max[error(MCMC)] max[error(MCMC)/SD] min(ESS) max(R-hat)"
-"     1      exp-aft                         Spike(0)          0.00076               0.011     8550      1.000"
-"     2  weibull-aft                         Spike(0)          0.00076               0.011     8993      1.000"
-"     3    lnorm-aft                         Spike(0)          0.00103               0.011     7831      1.000"
-"     4   llogis-aft                         Spike(0)          0.00089               0.011     8656      1.001"
-"     5    gamma-aft                         Spike(0)          0.00302               0.024     1725      1.001"
+"     1      exp-aft     treatment contrast: Spike(0)          0.00076               0.011     8550      1.000"
+"     2  weibull-aft     treatment contrast: Spike(0)          0.00076               0.011     8993      1.000"
+"     3    lnorm-aft     treatment contrast: Spike(0)          0.00103               0.011     7831      1.000"
+"     4   llogis-aft     treatment contrast: Spike(0)          0.00089               0.011     8656      1.001"
+"     5    gamma-aft     treatment contrast: Spike(0)          0.00302               0.024     1725      1.001"
 "     6      exp-aft treatment contrast: Normal(0, 1)          0.00268               0.019     2684      1.000"
 "     7  weibull-aft treatment contrast: Normal(0, 1)          0.00268               0.019     2804      1.001"
 "     8    lnorm-aft treatment contrast: Normal(0, 1)          0.00316               0.019     2903      1.000"