Merge branch 'master' into re-predictors

CRAN-SUBMISSION

@@ -1,3 +1,3 @@
-Version: 2.21.0
-Date: 2024-03-19 21:32:37 UTC
-SHA: 05cfbb6450960ac042be60a23293a646cf83a94d
+Version: 2.22.0
+Date: 2024-09-23 05:31:28 UTC
+SHA: 6c551b09bf3144ab049124f1147a0fb3c5605dd0

DESCRIPTION

@@ -2,8 +2,8 @@ Package: brms
 Encoding: UTF-8
 Type: Package
 Title: Bayesian Regression Models using 'Stan'
-Version: 2.21.11
-Date: 2024-09-19
+Version: 2.22.2
+Date: 2024-10-02
 Authors@R:
     c(person("Paul-Christian", "Bürkner", email = "paul.buerkner@gmail.com",
              role = c("aut", "cre")),
@@ -94,7 +94,7 @@ Description: Fit Bayesian generalized (non-)linear multivariate multilevel model
 LazyData: true
 NeedsCompilation: no
 License: GPL-2
-URL: https://github.com/paul-buerkner/brms, https://discourse.mc-stan.org/, https://paul-buerkner.github.io/brms/
+URL: https://github.com/paul-buerkner/brms, https://discourse.mc-stan.org/, https://paulbuerkner.com/brms/
 BugReports: https://github.com/paul-buerkner/brms/issues
 Additional_repositories:
     https://stan-dev.r-universe.dev/

NAMESPACE

@@ -285,6 +285,7 @@ S3method(stancode,default)
 S3method(standata,brmsfit)
 S3method(standata,default)
 S3method(stanplot,brmsfit)
+S3method(subset,psis)
 S3method(summarise_families,brmsformula)
 S3method(summarise_families,mvbrmsformula)
 S3method(summarise_links,brmsformula)

NEWS.md

@@ -1,4 +1,4 @@
-# brms 2.21.0++
+# brms 2.22.0
 
 ### New Features
 
@@ -9,6 +9,9 @@
 * Add priorsense support via `create_priorsense_data.brmsfit` 
 thanks to Noa Kallioinen. (#1354)
 * Vectorize censored log likelihoods in the Stan code when possible. (#1657)
+* Force Stan to activate threading without altering the Stan code
+via argument `force` of function `threading`. (#1549)
+* Support moment matching `loo` prediction methods. (#1674)
 
 ### Bug Fixes
 
@@ -17,11 +20,14 @@ thanks to Henrik Singmann. (#1651)
 * Fix problems with parallel executions of post-processing functions 
 sometimes leaving unused R instances behind. Thanks to Andrew Johnson,
 Aki Vehtari, and Noa Kallioinen. (#1658)
+* Fix several minor bugs. (#1648, #1644, #1672, #1642, #1634, #1666, #1664)
 
 ### Other Changes
 
 * Refactor some of the internal code base to avoid evaluating
 many data-dependent quantities several times. (#1653)
+* Smartly access internal functions when evaluating non-linear formulas. (#1635)
+* Improve the documentation in several places.
 * Make argument `loo` optional in `loo_moment_match`.
 * Change the output format of `loo_predict` and `loo_linpred` to be
 more consistent with other post-processing functions.

R/backends.R

@@ -70,7 +70,7 @@ compile_model <- function(model, backend, ...) {
   if (silent < 2) {
     message("Compiling Stan program...")
   }
-  if (use_threading(threads)) {
+  if (use_threading(threads, force = TRUE)) {
     if (utils::packageVersion("rstan") >= "2.26") {
       threads_per_chain_def <- rstan::rstan_options("threads_per_chain")
       on.exit(rstan::rstan_options(threads_per_chain = threads_per_chain_def))
@@ -100,7 +100,7 @@ compile_model <- function(model, backend, ...) {
   # if (cmdstanr::cmdstan_version() >= "2.29.0") {
   #   .canonicalize_stan_model(args$stan_file, overwrite_file = TRUE)
   # }
-  if (use_threading(threads)) {
+  if (use_threading(threads, force = TRUE)) {
     args$cpp_options$stan_threads <- TRUE
   }
   if (use_opencl(opencl)) {
@@ -147,7 +147,7 @@ fit_model <- function(model, backend, ...) {
                              seed, control, silent, future, ...) {
 
   # some input checks and housekeeping
-  if (use_threading(threads)) {
+  if (use_threading(threads, force = TRUE)) {
     if (utils::packageVersion("rstan") >= "2.26") {
       threads_per_chain_def <- rstan::rstan_options("threads_per_chain")
       on.exit(rstan::rstan_options(threads_per_chain = threads_per_chain_def))
@@ -265,6 +265,7 @@ fit_model <- function(model, backend, ...) {
   if (silent < 2) {
     message("Start sampling")
   }
+  use_threading <- use_threading(threads, force = TRUE)
   if (algorithm %in% c("sampling", "fixed_param")) {
     c(args) <- nlist(
       iter_sampling = iter - warmup,
@@ -275,7 +276,7 @@ fit_model <- function(model, backend, ...) {
       show_exceptions = silent == 0,
       fixed_param = algorithm == "fixed_param"
     )
-    if (use_threading(threads)) {
+    if (use_threading) {
       args$threads_per_chain <- threads$threads
     }
     if (future) {
@@ -304,17 +305,17 @@ fit_model <- function(model, backend, ...) {
     }
   } else if (algorithm %in% c("fullrank", "meanfield")) {
     c(args) <- nlist(iter, algorithm)
-    if (use_threading(threads)) {
+    if (use_threading) {
       args$threads <- threads$threads
     }
     out <- do_call(model$variational, args)
   } else if (algorithm %in% c("pathfinder")) {
-    if (use_threading(threads)) {
+    if (use_threading) {
       args$num_threads <- threads$threads
     }
     out <- do_call(model$pathfinder, args)
   } else if (algorithm %in% c("laplace")) {
-    if (use_threading(threads)) {
+    if (use_threading) {
       args$threads <- threads$threads
     }
     out <- do_call(model$laplace, args)
@@ -487,6 +488,10 @@ require_backend <- function(backend, x) {
 #'   \code{reduce_sum}? Defaults to \code{FALSE}. Setting it to \code{TRUE}
 #'   is required to achieve exact reproducibility of the model results
 #'   (if the random seed is set as well).
+#' @param force Logical. Defaults to \code{FALSE}. If \code{TRUE}, this will
+#'   force the Stan model to compile with threading enabled without altering the
+#'   Stan code generated by brms. This can be useful if your own custom Stan
+#'   functions use threading internally.
 #'
 #' @return A \code{brmsthreads} object which can be passed to the
 #'   \code{threads} argument of \code{brm} and related functions.
@@ -515,7 +520,8 @@ require_backend <- function(backend, x) {
 #' }
 #'
 #' @export
-threading <- function(threads = NULL, grainsize = NULL, static = FALSE) {
+threading <- function(threads = NULL, grainsize = NULL, static = FALSE,
+                      force = FALSE) {
   out <- list(threads = NULL, grainsize = NULL)
   class(out) <- "brmsthreads"
   if (!is.null(threads)) {
@@ -533,6 +539,7 @@ threading <- function(threads = NULL, grainsize = NULL, static = FALSE) {
     out$grainsize <- grainsize
   }
   out$static <- as_one_logical(static)
+  out$force <- as_one_logical(force)
   out
 }
 
@@ -555,8 +562,14 @@ validate_threads <- function(threads) {
 }
 
 # is threading activated?
-use_threading <- function(threads) {
-  isTRUE(validate_threads(threads)$threads > 0)
+use_threading <- function(threads, force = FALSE) {
+  threads <- validate_threads(threads)
+  out <- isTRUE(threads$threads > 0)
+  if (!force) {
+    # Stan code will only be altered in non-forced mode
+    out <- out && !isTRUE(threads$force)
+  }
+  out
 }
 
 #' GPU support in Stan via OpenCL

R/brm_multiple.R

@@ -158,6 +158,8 @@ brm_multiple <- function(formula, data, family = gaussian(), prior = NULL,
   if (combine) {
     fits <- combine_models(mlist = fits, check_data = FALSE)
     attr(fits$data, "data_name") <- data_name
+    # attribute to remember how many imputed datasets where used
+    attr(fits, "nimp") <- length(data)
     class(fits) <- c("brmsfit_multiple", class(fits))
   }
   if (!is.null(file)) {

R/brmsterms.R

@@ -131,7 +131,7 @@ brmsterms.brmsformula <- function(formula, check_response = TRUE,
       y$dpars[[dp]]$respform <- y$respform
       y$dpars[[dp]]$adforms <- y$adforms
     }
-    y$dpars[[dp]]$transform <- stan_eta_transform(y$dpars[[dp]]$family, y)
+    y$dpars[[dp]]$transform <- stan_eta_transform(y, y$dpars[[dp]]$family)
     check_cs(y$dpars[[dp]])
   }
 

R/conditional_effects.R

@@ -97,16 +97,18 @@
 #'   If \code{NULL} (default), \code{ncol} is computed internally based
 #'   on the number of rows of \code{conditions}.
 #' @param points Logical. Indicates if the original data points should be added
-#'   via \code{\link{geom_jitter}}. Default is \code{FALSE}. Can be controlled
-#'   globally via the \code{brms.plot_points} option. Note that only those data
-#'   points will be added that match the specified conditions defined in
-#'   \code{conditions}. For categorical predictors, the conditions have to match
-#'   exactly. For numeric predictors, argument \code{select_points} is used to
-#'   determine, which points do match a condition.
+#'   via \code{\link[ggplot2:geom_jitter]{geom_jitter}}. Default is
+#'   \code{FALSE}. Can be controlled globally via the \code{brms.plot_points}
+#'   option. Note that only those data points will be added that match the
+#'   specified conditions defined in \code{conditions}. For categorical
+#'   predictors, the conditions have to match exactly. For numeric predictors,
+#'   argument \code{select_points} is used to determine, which points do match a
+#'   condition.
 #' @param rug Logical. Indicates if a rug representation of predictor values
-#'   should be added via \code{\link{geom_rug}}. Default is \code{FALSE}.
-#'   Depends on \code{select_points} in the same way as \code{points} does. Can
-#'   be controlled globally via the \code{brms.plot_rug} option.
+#'   should be added via \code{\link[ggplot2:geom_rug]{geom_rug}}. Default is
+#'   \code{FALSE}. Depends on \code{select_points} in the same way as
+#'   \code{points} does. Can be controlled globally via the \code{brms.plot_rug}
+#'   option.
 #' @param mean Logical. Only relevant for spaghetti plots.
 #'   If \code{TRUE} (the default), display the mean regression
 #'   line on top of the regression lines for each sample.
@@ -118,30 +120,30 @@
 #'   Either \code{"contour"} or \code{"raster"}.
 #' @param line_args Only used in plots of continuous predictors:
 #'   A named list of arguments passed to
-#'   \code{\link{geom_smooth}}.
+#'   \code{\link[ggplot2:geom_smooth]{geom_smooth}}.
 #' @param cat_args Only used in plots of categorical predictors:
 #'   A named list of arguments passed to
-#'   \code{\link{geom_point}}.
+#'   \code{\link[ggplot2:geom_point]{geom_point}}.
 #' @param errorbar_args Only used in plots of categorical predictors:
 #'   A named list of arguments passed to
-#'   \code{\link{geom_errorbar}}.
+#'   \code{\link[ggplot2:geom_errorbar]{geom_errorbar}}.
 #' @param surface_args Only used in surface plots:
 #'   A named list of arguments passed to
-#'   \code{\link{geom_contour}} or
-#'   \code{\link{geom_raster}}
+#'   \code{\link[ggplot2:geom_contour]{geom_contour}} or
+#'   \code{\link[ggplot2:geom_raster]{geom_raster}}
 #'   (depending on argument \code{stype}).
 #' @param spaghetti_args Only used in spaghetti plots:
 #'   A named list of arguments passed to
-#'   \code{\link{geom_smooth}}.
+#'   \code{\link[ggplot2:geom_smooth]{geom_smooth}}.
 #' @param point_args Only used if \code{points = TRUE}:
 #'   A named list of arguments passed to
-#'   \code{\link{geom_jitter}}.
+#'   \code{\link[ggplot2:geom_jitter]{geom_jitter}}.
 #' @param rug_args Only used if \code{rug = TRUE}:
 #'   A named list of arguments passed to
-#'   \code{\link{geom_rug}}.
-#' @param facet_args Only used if if multiple condtions are provided:
+#'   \code{\link[ggplot2:geom_rug]{geom_rug}}.
+#' @param facet_args Only used if if multiple conditions are provided:
 #'   A named list of arguments passed to
-#'   \code{\link{facet_wrap}}.
+#'   \code{\link[ggplot2:facet_wrap]{facet_wrap}}.
 #'
 #' @return An object of class \code{'brms_conditional_effects'} which is a
 #'   named list with one data.frame per effect containing all information
@@ -154,7 +156,7 @@
 #'   rows).
 #'
 #'   The corresponding \code{plot} method returns a named
-#'   list of \code{\link{ggplot}} objects, which can be further
+#'   list of \code{\link[ggplot2:ggplot]{ggplot}} objects, which can be further
 #'   customized using the \pkg{ggplot2} package.
 #'
 #' @details When creating \code{conditional_effects} for a particular predictor
@@ -170,7 +172,8 @@
 #'
 #'   To fully change colors of the created plots, one has to amend both
 #'   \code{scale_colour} and \code{scale_fill}. See
-#'   \code{\link{scale_colour_grey}} or \code{\link{scale_colour_gradient}} for
+#'   \code{\link[ggplot2:scale_color_grey]{scale_colour_grey}} or
+#'   \code{\link[ggplot2:scale_color_gradient]{scale_colour_gradient}} for
 #'   more details.
 #'
 #' @examples

R/conditional_smooths.R

@@ -9,6 +9,10 @@
 #' @param smooths Optional character vector of smooth terms
 #'   to display. If \code{NULL} (the default) all smooth terms
 #'   are shown.
+#' @param surface Logical. Indicates if interactions or
+#'   two-dimensional smooths should be visualized as a surface.
+#'   Defaults to \code{TRUE}. The surface type can be controlled
+#'   via argument \code{stype} of the related plotting method.
 #' @param ndraws Positive integer indicating how many
 #'   posterior draws should be used.
 #'   If \code{NULL} (the default) all draws are used.
@@ -49,12 +53,14 @@
 conditional_smooths.brmsfit <- function(x, smooths = NULL,
                                         int_conditions = NULL,
                                         prob = 0.95, spaghetti = FALSE,
+                                        surface = TRUE,
                                         resolution = 100, too_far = 0,
                                         ndraws = NULL, draw_ids = NULL,
                                         nsamples = NULL, subset = NULL,
                                         probs = NULL, ...) {
   probs <- validate_ci_bounds(prob, probs = probs)
   spaghetti <- as_one_logical(spaghetti)
+  surface <- as_one_logical(surface)
   draw_ids <- use_alias(draw_ids, subset)
   ndraws <- use_alias(ndraws, nsamples)
   contains_draws(x)
@@ -68,7 +74,7 @@ conditional_smooths.brmsfit <- function(x, smooths = NULL,
     bterms, fit = x, smooths = smooths,
     conditions = conditions, int_conditions = int_conditions,
     too_far = too_far, resolution = resolution, probs = probs,
-    spaghetti = spaghetti, draw_ids = draw_ids
+    spaghetti = spaghetti, surface = surface, draw_ids = draw_ids
   )
   if (!length(out)) {
     stop2("No valid smooth terms found in the model.")
@@ -118,7 +124,7 @@ conditional_smooths.brmsterms <- function(x, ...) {
 #' @export
 conditional_smooths.btl <- function(x, fit, smooths, conditions, int_conditions,
                                     probs, resolution, too_far, spaghetti,
-                                    ...) {
+                                    surface, ...) {
   stopifnot(is.brmsfit(fit))
   out <- list()
   mf <- model.frame(fit)
@@ -149,17 +155,24 @@ conditional_smooths.btl <- function(x, fit, smooths, conditions, int_conditions,
     is_numeric <- setNames(rep(FALSE, ncovars), covars)
     for (cv in covars) {
       is_numeric[cv] <- is.numeric(mf[[cv]])
+      is_second_covar <- isTRUE(cv == covars[2])
       if (cv %in% names(int_conditions)) {
         int_cond <- int_conditions[[cv]]
         if (is.function(int_cond)) {
           int_cond <- int_cond(mf[[cv]])
         }
         values[[cv]] <- int_cond
       } else if (is_numeric[cv]) {
-        values[[cv]] <- seq(
-          min(mf[[cv]]), max(mf[[cv]]),
-          length.out = resolution
-        )
+        if (!surface && is_second_covar) {
+          mean2 <- mean(mf[[cv]], na.rm = TRUE)
+          sd2 <- sd(mf[[cv]], na.rm = TRUE)
+          values[[cv]] <- (-1:1) * sd2 + mean2
+        } else {
+          values[[cv]] <- seq(
+            min(mf[[cv]]), max(mf[[cv]]),
+            length.out = resolution
+          )
+        }
       } else {
         values[[cv]] <- levels(factor(mf[[cv]]))
       }
@@ -180,7 +193,8 @@ conditional_smooths.btl <- function(x, fit, smooths, conditions, int_conditions,
       }
     }
     newdata <- expand.grid(values)
-    if (ncovars == 2L && too_far > 0) {
+    need_surface <- surface && ncovars == 2L && all(is_numeric)
+    if (need_surface && too_far > 0) {
       # exclude prediction grid points too far from data
       ex_too_far <- mgcv::exclude.too.far(
         g1 = newdata[[covars[1]]],
@@ -196,6 +210,18 @@ conditional_smooths.btl <- function(x, fit, smooths, conditions, int_conditions,
     eta <- posterior_smooths(x, fit, smooth, newdata, ...)
     effects <- na.omit(sub_smframe$covars[[1]][1:2])
     cond_data <- add_effects__(newdata[, vars, drop = FALSE], effects)
+    second_numeric <- isTRUE(is_numeric[2])
+    if (second_numeric && !surface) {
+      # only convert 'effect2__' to factor so that the original
+      # second effect variable remains unchanged in the data
+      mde2 <- round(cond_data[[effects[2]]], 2)
+      levels2 <- sort(unique(mde2), TRUE)
+      cond_data$effect2__ <- factor(mde2, levels = levels2)
+      labels2 <- names(int_conditions[[effects[2]]])
+      if (length(labels2) == length(levels2)) {
+        levels(cond_data$effect2__) <- labels2
+      }
+    }
     if (length(byvars)) {
       # byvars will be plotted as facets
       cond_data$cond__ <- rows2labels(cond_data[, byvars, drop = FALSE])
@@ -218,7 +244,7 @@ conditional_smooths.btl <- function(x, fit, smooths, conditions, int_conditions,
     points <- add_effects__(points, covars)
     attr(eta, "response") <- response
     attr(eta, "effects") <- effects
-    attr(eta, "surface") <- all(is_numeric) && ncovars == 2L
+    attr(eta, "surface") <- need_surface
     attr(eta, "spaghetti") <- spa_data
     attr(eta, "points") <- points
     out[[response]] <- eta