JASP updates (#5)

* `check_setup()` function

* `diagnostics()` function

* `diagnostics()` function (II)

* `rho2logsdr` transformations

* forward BayesTools transformation changes

* update figures

* small fixes

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: RoBTT
 Title: Robust Bayesian T-Test
-Version: 1.0.3
+Version: 1.1.0
 Maintainer: František Bartoš <f.bartos96@gmail.com>
 Authors@R: c( 
     person("František", "Bartoš",     role = c("aut", "cre"),
@@ -26,7 +26,7 @@ Encoding: UTF-8
 LazyData: true
 SystemRequirements: GNU make
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.1
+RoxygenNote: 7.2.3
 NeedsCompilation: yes
 ByteCompile: true
 LinkingTo: 
@@ -43,9 +43,10 @@ Imports:
     rstan(>= 2.21.2),
     rstantools(>= 1.5.0),
     RcppParallel (>= 5.0.1),
-    BayesTools (>= 0.2.12),
+    BayesTools (>= 0.2.14),
     bridgesampling,
     methods,
+    ggplot2,
     Rdpack
 Suggests:
     parallel,

NAMESPACE

@@ -8,9 +8,15 @@ export(RoBTT)
 export(RoBTT.get_option)
 export(RoBTT.options)
 export(check_RoBTT)
+export(check_setup)
+export(diagnostics)
+export(diagnostics_autocorrelation)
+export(diagnostics_density)
+export(diagnostics_trace)
 export(interpret)
 export(is.RoBTT)
 export(prior)
+export(rho2logsdr)
 export(set_control)
 export(set_convergence_checks)
 import(Rcpp)

NEWS.md

@@ -1,3 +1,9 @@
+## version 1.1.0
+### Features
+- `check_setup()` function
+- `diagnostics()` function
+- `rho2logsdr` transformations
+
 ## version 1.0.3
 Compatibility update for the rstan 2.31.
 (correct CRAN commit)

R/RoBTT-package.R

@@ -27,3 +27,6 @@
 ##' @import rstantools
 ##' @useDynLib RoBTT,.registration = TRUE 
 "_PACKAGE"
+
+# hack for removing note from using quote in '.plot.RoBTT_par_names'
+utils::globalVariables(c("delta", "mu", "nu", "rho", "sigma"))

R/check-input-and-settings.R

@@ -1,3 +1,123 @@
+#' @title Prints summary of \code{"RoBTT"} ensemble implied by the specified priors
+#'
+#' @description \code{check_setup} prints summary of \code{"RoBTT"} ensemble
+#' implied by the specified prior distributions. It is useful for checking
+#' the ensemble configuration prior to fitting all of the models.
+#'
+#' @inheritParams RoBTT
+#' @param models should the models' details be printed.
+#' @param silent do not print the results.
+#'
+#' @return \code{check_setup} invisibly returns list of summary tables.
+#'
+#' @seealso [RoBTT()]
+#' @export
+check_setup <- function(
+    prior_delta  = prior(distribution = "cauchy",  parameters = list(location = 0, scale = sqrt(2)/2)),
+    prior_rho    = prior(distribution = "beta",    parameters = list(alpha = 1, beta = 1)),
+    prior_nu     = prior(distribution = "exp",     parameters = list(rate = 1)),
+
+    prior_delta_null  = prior(distribution = "spike",  parameters = list(location = 0)),
+    prior_rho_null    = prior(distribution = "spike",  parameters = list(location = 0.5)),
+    prior_nu_null     = NULL,
+
+    likelihood = c("normal", if(!is.null(prior_nu)) "t"),
+
+    models = FALSE, silent = FALSE){
+
+  object <- list()
+  object$priors      <- .set_priors(prior_delta, prior_rho, prior_nu, prior_delta_null, prior_rho_null, prior_nu_null)
+  object$models      <- .get_models(object$priors, likelihood)
+
+  ### model types overview
+  effect        <- sapply(object$models, function(m)!.is_parameter_null(m$priors, "delta"))
+  heterogeneity <- sapply(object$models, function(m)!.is_parameter_null(m$priors, "rho"))
+  outliers      <- sapply(object$models, function(m)!.is_parameter_null(m$priors, "nu"))
+
+  # number of model types
+  n_models    <- c(
+    effect        = sum(effect),
+    heterogeneity = sum(heterogeneity),
+    outliers      = sum(outliers)
+  )
+
+  # extract model weights
+  prior_weights   <- sapply(object$models, function(m) m$prior_weights)
+  # standardize model weights
+  prior_weights   <- prior_weights / sum(prior_weights)
+  # conditional model weights
+  models_prior <- c(
+    effect         = sum(prior_weights[effect]),
+    heterogeneity  = sum(prior_weights[heterogeneity]),
+    outliers       = sum(prior_weights[outliers])
+  )
+
+  # create overview table
+  components <- data.frame(
+    "models"     = n_models,
+    "prior_prob" = models_prior
+  )
+  rownames(components) <- c("Effect", "Heterogeneity", "Outliers")
+
+  class(components)             <- c("BayesTools_table", "BayesTools_ensemble_summary", class(components))
+  attr(components, "type")      <- c("n_models", "prior_prob")
+  attr(components, "rownames")  <- TRUE
+  attr(components, "n_models")  <- length(object$models)
+  attr(components, "title")     <- "Components summary:"
+  attr(components, "footnotes") <- NULL
+  attr(components, "warnings")  <- NULL
+
+  object$components <- components
+
+  ### model details
+  if(models){
+    priors_effect        <- sapply(1:length(object$models), function(i) print(object$models[[i]]$priors$delta, silent = TRUE))
+    priors_heterogeneity <- sapply(1:length(object$models), function(i) print(object$models[[i]]$priors$rho,   silent = TRUE))
+    priors_outliers      <- sapply(1:length(object$models), function(i) {
+      if(is.null(object$models[[i]]$priors$nu)){
+        return("")
+      }else{
+        print(object$models[[i]]$priors$nu,    silent = TRUE)
+      }
+    })
+
+    prior_weights    <- sapply(1:length(object$models), function(i)object$models[[i]]$prior_weights)
+    prior_prob       <- prior_weights / sum(prior_weights)
+    model_likelihood <- sapply(object[["models"]], function(m) m[["likelihood"]])
+
+    summary <- cbind.data.frame(
+      "Model"         = 1:length(object$models),
+      "Distribution"  = model_likelihood,
+      "delta"         = priors_effect,
+      "rho"           = priors_heterogeneity,
+      "nu"            = priors_outliers,
+      "prior_prob"    = prior_prob
+    )
+    class(summary)             <- c("BayesTools_table", "BayesTools_ensemble_inference", class(summary))
+    attr(summary, "type")      <- c("integer", "string",  rep("prior", 3), "prior_prob")
+    attr(summary, "rownames")  <- FALSE
+    attr(summary, "title")     <- "Models overview:"
+    attr(summary, "footnotes") <- NULL
+    attr(summary, "warnings")  <- NULL
+
+    object$summary <- summary
+  }
+
+
+  if(!silent){
+    cat("Robust Bayesian t-test (set-up)\n")
+    print(components, quote = FALSE, right = TRUE)
+
+    if(models){
+      cat("\n")
+      print(summary, quote = FALSE, right = TRUE)
+    }
+  }
+
+  return(invisible(object))
+}
+
+
 #' @title Convergence checks of the fitting process
 #'
 #' @description Set values for the convergence checks of the fitting process.

R/diagnostics.R

@@ -0,0 +1,180 @@
+#' @title Checks a fitted RoBTT object
+#'
+#' @description \code{diagnostics} creates visual
+#' checks of individual models convergence. Numerical
+#' overview of individual models can be obtained by
+#' \code{summary(object, type = "models", diagnostics = TRUE)},
+#' or even more detailed information by
+#' \code{summary(object, type = "individual")}.
+#'
+#' @param fit a fitted RoBTT object
+#' @param parameter a parameter to be plotted. Either
+#' \code{"delta"}, \code{"rho"}, \code{"nu"}, \code{"mu"},
+#' or \code{"sigma"}.
+#' @param type type of MCMC diagnostic to be plotted.
+#' Options are \code{"chains"} for the chains' trace plots,
+#' \code{"autocorrelation"} for autocorrelation of the
+#' chains, and \code{"densities"} for the overlaying
+#' densities of the individual chains. Can be abbreviated to
+#' first letters.
+#' @param show_models MCMC diagnostics of which models should be
+#' plotted. Defaults to \code{NULL} which plots MCMC diagnostics
+#' for a specified parameter for every model that is part of the
+#' ensemble.
+#' @param title whether the model number should be displayed in title.
+#' Defaults to \code{TRUE} when more than one model is selected.
+#' @param lags number of lags to be shown for
+#' \code{type = "autocorrelation"}. Defaults to \code{30}.
+#' @param ... additional arguments to be passed to
+#' \link[graphics]{par} if \code{plot_type = "base"}.
+#'
+#' @details The visualization functions are based on
+#' \link[rstan]{stan_plot} function and its color schemes.
+#'
+#' @examples \dontrun{
+#' # using the example data from Darwin
+#' data("fertilization", package = "RoBTT")
+#' fit <- RoBTT(
+#'   x1       = fertilization$Self,
+#'   x2       = fertilization$Crossed,
+#'   prior_delta = prior("cauchy", list(0, 1/sqrt(2))),
+#'   prior_rho   = prior("beta",   list(3, 3)),
+#'   likelihood  = "normal",
+#'   seed        = 1, 
+#'   chains      = 1,
+#'   warmup      = 1000,
+#'   iter        = 2000,
+#'   control     = set_control(adapt_delta = 0.95)
+#' )
+#'
+#' ### ggplot2 version of all of the plots can be obtained by adding 'model_type = "ggplot"
+#' # diagnostics function allows to visualize diagnostics of a fitted RoBTT object, for example,
+#' # the trace plot for the mean parameter in each model model
+#' diagnostics(fit, parameter = "delta", type = "chain")
+#'
+#' # in order to show the trace plot only for the 11th model, add show_models parameter
+#' diagnostics(fit, parameter = "delta", type = "chain", show_models = 11)
+#'
+#' # furthermore, the autocorrelations
+#' diagnostics(fit, parameter = "delta", type = "autocorrelation")
+#'
+#' # and overlying densities for each plot can also be visualize
+#' diagnostics(fit, parameter = "delta", type = "densities")
+#' }
+#'
+#'
+#' @return \code{diagnostics} returns either \code{NULL} if \code{plot_type = "base"}
+#' or an object/list of objects (depending on the number of parameters to be plotted)
+#' of class 'ggplot2' if \code{plot_type = "ggplot2"}.
+#'
+#' @seealso [RoBTT()], [summary.RoBTT()]
+#' @name diagnostics
+#' @aliases diagnostics_autocorrelation diagnostics_trace diagnostics_density
+#' @export diagnostics
+#' @export diagnostics_density
+#' @export diagnostics_autocorrelation
+#' @export diagnostics_trace
+
+#' @rdname diagnostics
+diagnostics <- function(fit, parameter, type, show_models = NULL,
+                        lags = 30, title = is.null(show_models) | length(show_models) > 1, ...){
+
+  # check settings
+  if(!is.RoBTT(fit))
+    stop("Diagnostics are available only for RoBTT models.")
+  BayesTools::check_char(parameter, "parameter")
+  BayesTools::check_char(type, "type")
+
+  # deal with bad type names
+  if(substr(type, 1, 1) == "c"){
+    type <- "trace"
+  }else if(substr(type, 1, 1) == "t"){
+    type <- "trace" # for trace
+  }else if(substr(type, 1, 1) == "d"){
+    type <- "density"
+  }else if(substr(type, 1, 1) == "a"){
+    type <- "autocorrelation"
+  }else{
+    stop("Unsupported diagnostics type. See '?diagnostics' for more details.")
+  }
+
+  if(parameter %in% c("delta", "rho", "nu", "mu")){
+    parameter         <- parameter
+    parameter_samples <- parameter
+  }else if(parameter == "sigma"){
+    parameter         <- "pooled_sigma"
+    parameter_samples <- "pooled_sigma"
+
+  }else{
+    stop(paste0("The passed parameter does not correspond to any of main model parameter ('delta', 'rho', 'nu', 'mu', 'sigma'). See '?plot.RoBTT' for more details."))
+  }
+
+  # do the plotting
+  models_ind <- 1:length(fit$models)
+  if(!is.null(show_models)){
+    models_ind <- models_ind[show_models]
+  }
+
+  plots <- list()
+
+  for(i in models_ind){
+
+    prior_names <- names(attr(fit$models[[i]][["fit"]], "prior_list"))
+    prior_names <- prior_names[!sapply(attr(fit$models[[i]][["fit"]], "prior_list"), BayesTools::is.prior.point)]
+    model_parameters <- c(prior_names, "mu", "pooled_sigma")
+
+    if(!parameter_samples %in% model_parameters){
+
+      plots[[i]] <- NULL
+
+    }else if(inherits(fit$models[[i]][["fit"]], "null_model")){
+
+      plots[[i]] <- NULL
+
+    }else{
+
+      if(type == "density"){
+        plots[[i]] <- rstan::stan_dens(fit$models[[i]]$fit, pars = parameter) + 
+          ggplot2::ylab("Density") + 
+          ggplot2::xlab(.plot.RoBTT_par_names(parameter))
+      }else if(type == "trace"){
+        plots[[i]] <- rstan::traceplot(fit$models[[i]]$fit, pars = parameter) + 
+          ggplot2::ylab(.plot.RoBTT_par_names(parameter)) + 
+          ggplot2::xlab("Iteration")
+      }else if(type == "autocorrelation"){
+        plots[[i]] <- rstan::stan_ac(fit$models[[i]]$fit, pars = parameter, lags = lags, separate_chains = TRUE) + 
+          ggplot2::ylab(substitute("Autocorrelation"~(x), list(x = .plot.RoBTT_par_names(parameter, quote = TRUE)))) +
+          ggplot2::xlab("Lag")
+      }
+    }
+  }
+
+
+  # return the plots
+  if(length(plots) == 0){
+    plots <- NULL
+  }else if(length(models_ind) == 1){
+    plots <- plots[[models_ind]]
+  }
+  return(plots)
+}
+
+
+#' @rdname diagnostics
+diagnostics_autocorrelation <- function(fit, parameter = NULL, show_models = NULL,
+                                        lags = 30, title = is.null(show_models) | length(show_models) > 1, ...){
+  diagnostics(fit = fit, parameter = parameter, type = "autocorrelation", show_models = show_models, lags = lags, title = title, ...)
+}
+
+#' @rdname diagnostics
+diagnostics_trace           <- function(fit, parameter = NULL, show_models = NULL,
+                                        title = is.null(show_models) | length(show_models) > 1, ...){
+  diagnostics(fit = fit, parameter = parameter, type = "trace", show_models = show_models, title = title, ...)
+}
+
+#' @rdname diagnostics
+diagnostics_density         <- function(fit, parameter = NULL, show_models = NULL,
+                                        title = is.null(show_models) | length(show_models) > 1, ...){
+  diagnostics(fit = fit, parameter = parameter, type = "density", show_models = show_models, title = title, ...)
+}
+