Merge branch 'master' into ridgeline-size

Author: jgabry

.github/FUNDING.yml

@@ -0,0 +1,2 @@
+github: stan-dev
+custom: https://mc-stan.org/support/

.github/workflows/R-CMD-check.yaml

@@ -19,11 +19,10 @@ jobs:
       fail-fast: false
       matrix:
         config:
-          - {os: macOS-latest,   r: 'devel'}
-          - {os: macOS-latest,   r: 'release'}
           - {os: windows-latest, r: 'release'}
-          - {os: ubuntu-16.04,   r: 'release', rspm: "https://packagemanager.rstudio.com/cran/__linux__/xenial/latest"}
-          - {os: ubuntu-16.04,   r: 'oldrel',  rspm: "https://packagemanager.rstudio.com/cran/__linux__/xenial/latest"}
+          - {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"}
     env:
       R_REMOTES_NO_ERRORS_FROM_WARNINGS: true
       RSPM: ${{ matrix.config.rspm }}

DESCRIPTION

@@ -49,7 +49,7 @@ Suggests:
     shinystan (>= 2.3.0),
     testthat (>= 2.0.0),
     vdiffr
-RoxygenNote: 7.1.0.9000
+RoxygenNote: 7.1.1
 VignetteBuilder: knitr
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)

NAMESPACE

@@ -2,10 +2,13 @@
 
 S3method("[",neff_ratio)
 S3method("[",rhat)
+S3method(log_posterior,CmdStanMCMC)
 S3method(log_posterior,stanfit)
 S3method(log_posterior,stanreg)
+S3method(neff_ratio,CmdStanMCMC)
 S3method(neff_ratio,stanfit)
 S3method(neff_ratio,stanreg)
+S3method(nuts_params,CmdStanMCMC)
 S3method(nuts_params,list)
 S3method(nuts_params,stanfit)
 S3method(nuts_params,stanreg)
@@ -15,6 +18,7 @@ S3method(pp_check,default)
 S3method(print,bayesplot_function_list)
 S3method(print,bayesplot_grid)
 S3method(print,bayesplot_scheme)
+S3method(rhat,CmdStanMCMC)
 S3method(rhat,stanfit)
 S3method(rhat,stanreg)
 export(abline_01)

NEWS.md

@@ -8,7 +8,14 @@
 * Items for next release go here
 -->
 
+* CmdStanMCMC objects (from CmdStanR) can now be used with extractor 
+  functions `nuts_params()`, `log_posterior()`, `rhat()`, and 
+  `neff_ratio()`. (#227)
+
 * Added missing `facet_args` argument to `mcmc_rank_overlay()`. (#221, @hhau)
+* Size of points and interval lines can set in 
+  `mcmc_intervals(..., outer_size, inner_size, point_size)`. (#215, #228, #229) 
+* `mcmc_areas()` tries to use less blank vertical blank space. (#218, #230) 
 
 
 # bayesplot 1.7.2

R/bayesplot-extractors.R

@@ -1,9 +1,9 @@
 #' Extract quantities needed for plotting from model objects
 #'
 #' Generics and methods for extracting quantities needed for plotting from
-#' various types of model objects. Currently methods are only provided for
-#' stanfit (**rstan**) and stanreg (**rstanarm**) objects, but adding new
-#' methods should be relatively straightforward.
+#' various types of model objects. Currently methods are provided for stanfit
+#' (**rstan**), CmdStanMCMC (**cmdstanr**), and stanreg (**rstanarm**) objects,
+#' but adding new methods should be relatively straightforward.
 #'
 #' @name bayesplot-extractors
 #' @param object The object to use.
@@ -87,7 +87,8 @@ log_posterior.stanfit <- function(object, inc_warmup = FALSE, ...) {
                                 ...)
   lp <- lapply(lp, as.array)
   lp <- set_names(reshape2::melt(lp), c("Iteration", "Value", "Chain"))
-  validate_df_classes(lp, c("integer", "numeric", "integer"))
+  validate_df_classes(lp[, c("Chain", "Iteration", "Value")],
+                      c("integer", "integer", "numeric"))
 }
 
 #' @rdname bayesplot-extractors
@@ -98,11 +99,22 @@ log_posterior.stanreg <- function(object, inc_warmup = FALSE, ...) {
   log_posterior.stanfit(object$stanfit, inc_warmup = inc_warmup, ...)
 }
 
+#' @rdname bayesplot-extractors
+#' @export
+#' @method log_posterior CmdStanMCMC
+log_posterior.CmdStanMCMC <- function(object, inc_warmup = FALSE, ...) {
+  lp <- object$draws("lp__", inc_warmup = inc_warmup)
+  lp <- reshape2::melt(lp)
+  lp$variable <- NULL
+  lp <- dplyr::rename_with(lp, capitalize_first)
+  validate_df_classes(lp[, c("Chain", "Iteration", "Value")],
+                      c("integer", "integer", "numeric"))
+}
+
 
 #' @rdname bayesplot-extractors
 #' @export
 #' @method nuts_params stanfit
-#'
 nuts_params.stanfit <-
   function(object,
            pars = NULL,
@@ -153,7 +165,23 @@ nuts_params.list <- function(object, pars = NULL, ...) {
 
   out <- reshape2::melt(object)
   out <- set_names(out, c("Iteration", "Parameter", "Value", "Chain"))
-  validate_df_classes(out, c("integer", "factor", "numeric", "integer"))
+  validate_df_classes(out[, c("Chain", "Iteration", "Parameter", "Value")],
+                      c("integer", "integer", "factor", "numeric"))
+}
+
+#' @rdname bayesplot-extractors
+#' @export
+#' @method nuts_params CmdStanMCMC
+nuts_params.CmdStanMCMC <- function(object, pars = NULL, ...) {
+  arr <- object$sampler_diagnostics()
+  if (!is.null(pars)) {
+    arr <- arr[,, pars]
+  }
+  out <- reshape2::melt(arr)
+  colnames(out)[colnames(out) == "variable"] <- "parameter"
+  out <- dplyr::rename_with(out, capitalize_first)
+  validate_df_classes(out[, c("Chain", "Iteration", "Parameter", "Value")],
+                      c("integer", "integer", "factor", "numeric"))
 }
 
 
@@ -188,6 +216,17 @@ rhat.stanreg <- function(object, pars = NULL, regex_pars = NULL, ...) {
   r[!names(r) %in% c("mean_PPD", "log-posterior")]
 }
 
+#' @rdname bayesplot-extractors
+#' @export
+#' @method rhat CmdStanMCMC
+rhat.CmdStanMCMC <- function(object, pars = NULL, ...) {
+  .rhat <- utils::getFromNamespace("rhat", "posterior")
+  s <- object$summary(pars, rhat = .rhat)[, c("variable", "rhat")]
+  r <- setNames(s$rhat, s$variable)
+  r <- validate_rhat(r)
+  r[!names(r) %in% "lp__"]
+}
+
 
 #' @rdname bayesplot-extractors
 #' @export
@@ -223,6 +262,18 @@ neff_ratio.stanreg <- function(object, pars = NULL, regex_pars = NULL, ...) {
   ratio[!names(ratio) %in% c("mean_PPD", "log-posterior")]
 }
 
+#' @rdname bayesplot-extractors
+#' @export
+#' @method neff_ratio CmdStanMCMC
+neff_ratio.CmdStanMCMC <- function(object, pars = NULL, ...) {
+  s <- object$summary(pars, "n_eff" = "ess_basic")[, c("variable", "n_eff")]
+  ess <- setNames(s$n_eff, s$variable)
+  tss <- prod(dim(object$draws())[1:2])
+  ratio <- ess / tss
+  ratio <- validate_neff_ratio(ratio)
+  ratio[!names(ratio) %in% "lp__"]
+}
+
 
 # internals ---------------------------------------------------------------
 
@@ -245,3 +296,10 @@ validate_df_classes <- function(x, classes = character()) {
   }
   x
 }
+
+# capitalize first letter in a string only
+capitalize_first <- function(name) {
+  name <- tolower(name) # in case whole string is capitalized
+  substr(name, 1, 1) <- toupper(substr(name, 1, 1))
+  name
+}

R/mcmc-diagnostics-nuts.R

@@ -513,8 +513,8 @@ validate_nuts_data_frame <- function(x, lp) {
     abort("NUTS parameters should be in a data frame.")
   }
 
-  valid_cols <- c("Iteration", "Parameter", "Value", "Chain")
-  if (!identical(colnames(x), valid_cols)) {
+  valid_cols <- sort(c("Iteration", "Parameter", "Value", "Chain"))
+  if (!identical(sort(colnames(x)), valid_cols)) {
     abort(paste(
       "NUTS parameter data frame must have columns:",
       paste(valid_cols, collapse = ", ")
@@ -529,8 +529,8 @@ validate_nuts_data_frame <- function(x, lp) {
       abort("lp should be in a data frame.")
     }
 
-    valid_lp_cols <- c("Iteration", "Value", "Chain")
-    if (!identical(colnames(lp), valid_lp_cols)) {
+    valid_lp_cols <- sort(c("Iteration", "Value", "Chain"))
+    if (!identical(sort(colnames(lp)), valid_lp_cols)) {
       abort(paste(
         "lp data frame must have columns:",
         paste(valid_lp_cols, collapse = ", ")

R/mcmc-intervals.R

@@ -25,6 +25,9 @@
 #'   `"equal height"` are scaled using `height*sqrt(height)`
 #' @param point_est The point estimate to show. Either `"median"` (the
 #'   default), `"mean"`, or `"none"`.
+#' @param inner_size,outer_size For `mcmc_intervals()`, the size of
+#'   the inner and interval segments, respectively.
+#' @param point_size For `mcmc_intervals()`, the size of point estimate.
 #' @param rhat An optional numeric vector of R-hat estimates, with one element
 #'   per parameter included in `x`. If `rhat` is provided, the intervals/areas
 #'   and point estimates in the resulting plot are colored based on R-hat value.
@@ -52,6 +55,10 @@
 #'
 #' @examples
 #' set.seed(9262017)
+#'
+#' # load ggplot2 to use its functions to modify our plots
+#' library(ggplot2)
+#'
 #' # some parameter draws to use for demonstration
 #' x <- example_mcmc_draws(params = 6)
 #' dim(x)
@@ -61,19 +68,31 @@
 #' mcmc_intervals(x)
 #' mcmc_intervals(x, pars = c("beta[1]", "beta[2]"))
 #' mcmc_areas(x, regex_pars = "beta\\[[1-3]\\]",  prob = 0.8) +
-#'  ggplot2::labs(
+#'  labs(
 #'    title = "Posterior distributions",
 #'    subtitle = "with medians and 80% intervals"
 #'  )
 #'
 #' color_scheme_set("red")
-#' mcmc_areas(
+#' p <- mcmc_areas(
 #'    x,
 #'    pars = c("alpha", "beta[4]"),
 #'    prob = 2/3,
 #'    prob_outer = 0.9,
 #'    point_est = "mean"
 #' )
+#' plot(p)
+#'
+#' # control spacing at top and bottom of plot
+#' # see ?ggplot2::expansion
+#' p + scale_y_discrete(
+#'   limits = c("beta[4]", "alpha"),
+#'   expand = expansion(add = c(1, 2))
+#' )
+#' p + scale_y_discrete(
+#'   limits = c("beta[4]", "alpha"),
+#'   expand = expansion(add = c(.1, .3))
+#' )
 #'
 #' # color by rhat value
 #' color_scheme_set("blue")
@@ -94,22 +113,22 @@
 #' b3 <- c("beta[1]", "beta[2]", "beta[3]")
 #'
 #' mcmc_areas(x, pars = b3, area_method = "equal area") +
-#'   ggplot2::labs(
+#'   labs(
 #'     title = "Curves have same area",
-#'     subtitle =
-#'       "A wide, uncertain interval is spread thin when areas are equal")
+#'     subtitle = "A wide, uncertain interval is spread thin when areas are equal"
+#'    )
 #'
 #' mcmc_areas(x, pars = b3, area_method = "equal height") +
-#'   ggplot2::labs(
+#'   labs(
 #'     title = "Curves have same maximum height",
-#'     subtitle =
-#'       "Local curvature is clearer but more uncertain curves use more area")
+#'     subtitle = "Local curvature is clearer but more uncertain curves use more area"
+#'   )
 #'
 #' mcmc_areas(x, pars = b3, area_method = "scaled height") +
-#'   ggplot2::labs(
+#'   labs(
 #'     title = "Same maximum heights but heights scaled by square-root",
-#'     subtitle =
-#'       "Compromise: Local curvature is accentuated and less area is used")
+#'     subtitle = "Compromise: Local curvature is accentuated and less area is used"
+#'    )
 #'
 #' \donttest{
 #' # apply transformations
@@ -145,7 +164,7 @@
 #' # plotted with ridgelines
 #' m <- shinystan::eight_schools@posterior_sample
 #' mcmc_areas_ridges(m, pars = "mu", regex_pars = "theta") +
-#'  ggplot2::ggtitle("Treatment effect on eight schools (Rubin, 1981)")
+#'   ggtitle("Treatment effect on eight schools (Rubin, 1981)")
 #' }
 #'
 NULL
@@ -160,6 +179,9 @@ mcmc_intervals <- function(x,
                            prob = 0.5,
                            prob_outer = 0.9,
                            point_est = c("median", "mean", "none"),
+                           outer_size = 0.5,
+                           inner_size = 2,
+                           point_size = 4,
                            rhat = numeric()) {
   check_ignored_arguments(...)
 
@@ -184,17 +206,18 @@ mcmc_intervals <- function(x,
 
   args_outer <- list(
     mapping = aes_(x = ~ ll, xend = ~ hh, y = ~ parameter, yend = ~ parameter),
-    color = get_color("mid")
+    color = get_color("mid"),
+    size = outer_size
   )
   args_inner <- list(
     mapping = aes_(x = ~ l, xend = ~ h, y = ~ parameter, yend = ~ parameter),
-    size = 2,
+    size = inner_size,
     show.legend = FALSE
   )
   args_point <- list(
     mapping = aes_(x = ~ m, y = ~ parameter),
     data = data,
-    size = 4,
+    size = point_size,
     shape = 21
   )
 
@@ -269,7 +292,8 @@ mcmc_areas <- function(x,
     x, pars, regex_pars, transformations,
     prob = prob, prob_outer = prob_outer,
     point_est = point_est, rhat = rhat,
-    bw = bw, adjust = adjust, kernel = kernel, n_dens = n_dens)
+    bw = bw, adjust = adjust, kernel = kernel, n_dens = n_dens
+  )
   datas <- split(data, data$interval)
 
   # Use a dummy empty dataframe if no point estimate
@@ -312,7 +336,11 @@ mcmc_areas <- function(x,
 
   datas$bottom <- datas$outer %>%
     group_by(!!! groups) %>%
-    summarise(ll = min(.data$x), hh = max(.data$x)) %>%
+    summarise(
+      ll = min(.data$x),
+      hh = max(.data$x),
+      .groups = "drop_last"
+    ) %>%
     ungroup()
 
   args_bottom <- list(
@@ -360,9 +388,16 @@ mcmc_areas <- function(x,
     args_outer$color <- get_color("dark")
   }
 
+  # An invisible layer that is 2.5% taller than the plotted one
+  args_outer2 <- args_outer
+  args_outer2$mapping <- args_outer2$mapping %>%
+    modify_aes_(scale = .925)
+  args_outer2$color <- NA
+
   layer_bottom <- do.call(geom_segment, args_bottom)
   layer_inner <- do.call(ggridges::geom_ridgeline, args_inner)
   layer_outer <- do.call(ggridges::geom_ridgeline, args_outer)
+  layer_outer2 <- do.call(ggridges::geom_ridgeline, args_outer2)
 
   point_geom <- if (no_point_est) {
     geom_ignore
@@ -386,12 +421,17 @@ mcmc_areas <- function(x,
     layer_inner +
     layer_point +
     layer_outer +
+    layer_outer2 +
     layer_bottom +
     scale_color +
     scale_fill +
     scale_y_discrete(
       limits = unique(rev(data$parameter)),
-      expand = expansion(add = c(0, .1), mult = c(.1, .3))) +
+      expand = expansion(
+        add = c(0, .5 + 1/(2 * nlevels(data$parameter))),
+        mult = c(.1, .1)
+      )
+    ) +
     xlim(x_lim) +
     bayesplot_theme_get() +
     legend_move(ifelse(color_by_rhat, "top", "none")) +