Merge pull request #154 from mcol/fix_152

Add prob_outer option to ppc_ribbon() and ppc_intervals()

Author: jgabry

NEWS.md

@@ -27,6 +27,13 @@ now has an argument `grid_args` to provide a way of passing optional arguments t
 gains an argument `discrete`, which is `FALSE` by default, but can be used to make the 
 Geom more appropriate for discrete data. (#145)
 
+* [PPC intervals plots](http://mc-stan.org/bayesplot/reference/PPC-intervals.html)
+and [LOO predictive checks](http://mc-stan.org/bayesplot/reference/PPC-loo.html)
+now draw both an outer and an inner probability interval, which can be
+controlled through the new argument `prob_outer` and the already existing
+`prob`. This is consistent with what is produced by `mcmc_intervals()`.
+(#152, #154, @mcol)
+
 
 # bayesplot 1.5.0
 

R/ppc-intervals.R

@@ -7,13 +7,12 @@
 #' @family PPCs
 #'
 #' @template args-y-yrep
+#' @template args-prob-prob_outer
 #' @param x A numeric vector the same length as \code{y} to use as the x-axis
 #'   variable. For example, \code{x} could be a predictor variable from a
 #'   regression model, a time variable for time-series models, etc. If \code{x}
 #'   is missing or NULL, then \code{1:length(y)} is used for the x-axis.
 #' @param ... Currently unused.
-#' @param prob A value between 0 and 1 indicating the desired probability mass
-#'   to include in the \code{yrep} intervals. The default is 0.9.
 #' @param alpha,size,fatten Arguments passed to geoms. For ribbon plots
 #'   \code{alpha} and \code{size} are passed to
 #'   \code{\link[ggplot2]{geom_ribbon}}. For interval plots \code{size} and
@@ -109,16 +108,17 @@ NULL
 
 #' @rdname PPC-intervals
 #' @export
-ppc_intervals <- function(y, yrep, x = NULL, ..., prob = 0.9, size = 1,
-                          fatten = 3) {
+ppc_intervals <- function(y, yrep, x = NULL, ..., prob = 0.5, prob_outer = 0.9,
+                          size = 1, fatten = 3) {
   check_ignored_arguments(...)
 
   data <- ppc_intervals_data(
       y = y,
       yrep = yrep,
       x = x,
       group = NULL,
-      prob = prob
+      prob = prob,
+      prob_outer = prob_outer
   )
 
   .ppc_intervals(
@@ -143,7 +143,8 @@ ppc_intervals_grouped <- function(y,
                                   group,
                                   facet_args = list(),
                                   ...,
-                                  prob = 0.9,
+                                  prob = 0.5,
+                                  prob_outer = 0.9,
                                   size = 1,
                                   fatten = 3) {
   check_ignored_arguments(...)
@@ -157,7 +158,8 @@ ppc_intervals_grouped <- function(y,
     yrep = yrep,
     x = x,
     group = group,
-    prob = prob
+    prob = prob,
+    prob_outer = prob_outer
   )
 
   .ppc_intervals(
@@ -178,7 +180,8 @@ ppc_ribbon <- function(y,
                        yrep,
                        x = NULL,
                        ...,
-                       prob = 0.9,
+                       prob = 0.5,
+                       prob_outer = 0.9,
                        alpha = 0.33,
                        size = 0.25) {
   check_ignored_arguments(...)
@@ -210,7 +213,8 @@ ppc_ribbon_grouped <- function(y,
                                group,
                                facet_args = list(),
                                ...,
-                               prob = 0.9,
+                               prob = 0.5,
+                               prob_outer = 0.9,
                                alpha = 0.33,
                                size = 0.25) {
   check_ignored_arguments(...)
@@ -224,7 +228,8 @@ ppc_ribbon_grouped <- function(y,
     yrep = yrep,
     x = x,
     group = group,
-    prob = prob
+    prob = prob,
+    prob_outer = prob_outer
   )
 
   .ppc_intervals(
@@ -241,9 +246,11 @@ ppc_ribbon_grouped <- function(y,
 
 #' @rdname PPC-intervals
 #' @export
-ppc_intervals_data <- function(y, yrep, x = NULL, group = NULL, prob = 0.9, ...) {
+ppc_intervals_data <- function(y, yrep, x = NULL, group = NULL,
+                               prob = 0.5, prob_outer = 0.9, ...) {
   check_ignored_arguments(...)
-  .ppc_intervals_data(y = y, yrep = yrep, x = x, group = group, prob = prob)
+  .ppc_intervals_data(y = y, yrep = yrep, x = x, group = group,
+                      prob = prob, prob_outer = prob_outer)
 }
 
 
@@ -259,9 +266,14 @@ label_x <- function(x) {
   if (missing(x)) "Index" else NULL
 }
 
-.ppc_intervals_data <- function(y, yrep, x = NULL, group = NULL, prob = 0.8) {
+.ppc_intervals_data <- function(y, yrep, x = NULL, group = NULL,
+                                prob = 0.5, prob_outer = 0.9) {
   grouped <- !is.null(group)
   stopifnot(prob > 0 && prob < 1)
+  stopifnot(prob_outer > 0 && prob_outer <= 1)
+  probs <- sort(c(prob, prob_outer))
+  prob <- probs[1]
+  prob_outer <- probs[2]
 
   y <- validate_y(y)
   yrep <- validate_yrep(yrep, y)
@@ -283,16 +295,19 @@ label_x <- function(x) {
   }
 
   grouped_d <- dplyr::group_by(molten_reps, !!! group_vars)
-  alpha <- (1 - prob) / 2
+  alpha <- (1 - probs) / 2
   probs <- sort(c(alpha, 0.5, 1 - alpha))
 
   val_col <- sym("value")
   dplyr::ungroup(dplyr::summarise(
     grouped_d,
-    prob = prob,
-    lo = quantile(!! val_col, prob = probs[1]),
-    mid = quantile(!! val_col, prob = probs[2]),
-    hi = quantile(!! val_col, prob = probs[3])
+    outer_width = prob_outer,
+    inner_width = prob,
+    ll = quantile(!! val_col, prob = probs[1]),
+    l  = quantile(!! val_col, prob = probs[2]),
+    m  = quantile(!! val_col, prob = probs[3]),
+    h  = quantile(!! val_col, prob = probs[4]),
+    hh = quantile(!! val_col, prob = probs[5])
   ))
 }
 
@@ -317,14 +332,19 @@ label_x <- function(x) {
     data = data,
     mapping = aes_(
       x = ~ x,
-      y = ~ mid,
-      ymin = ~ lo,
-      ymax = ~ hi
+      y = ~ m,
+      ymin = ~ l,
+      ymax = ~ h
     )
   )
 
   if (style == "ribbon") {
     graph <- graph +
+      geom_ribbon(
+        aes_(color = "yrep", fill = "yrep", ymin = ~ ll, ymax = ~ hh),
+        alpha = alpha,
+        size = size
+      ) +
       geom_ribbon(
         aes_(color = "yrep", fill = "yrep"),
         alpha = alpha,
@@ -341,6 +361,13 @@ label_x <- function(x) {
       )
   } else {
     graph <- graph +
+      geom_pointrange(
+        mapping = aes_(color = "yrep", fill = "yrep", ymin = ~ ll, ymax = ~ hh),
+        shape = 21,
+        alpha = alpha,
+        size = size,
+        fatten = fatten
+      ) +
       geom_pointrange(
         mapping = aes_(color = "yrep", fill = "yrep"),
         shape = 21,

R/ppc-loo.R

@@ -90,11 +90,11 @@
 #'
 #' # loo predictive intervals vs observations
 #' keep_obs <- 1:50
-#' ppc_loo_intervals(y, yrep, psis_object = psis1, subset = keep_obs, prob = 0.9)
+#' ppc_loo_intervals(y, yrep, psis_object = psis1, subset = keep_obs)
 #'
 #' color_scheme_set("gray")
 #' ppc_loo_intervals(y, yrep, psis_object = psis1, subset = keep_obs,
-#'                   order = "median", prob = 0.9)
+#'                   order = "median")
 #' }
 #'
 NULL
@@ -266,20 +266,19 @@ ppc_loo_pit <-
 
 #' @rdname PPC-loo
 #' @export
+#' @template args-prob-prob_outer
 #' @param psis_object If using \pkg{loo} version \code{2.0.0} or greater, an
 #'   object returned by the \code{psis} function (or by the \code{loo} function
 #'   with argument \code{save_psis} set to \code{TRUE}).
-#' @param prob A value between 0 and 1 indicating the desired probability mass
-#'   to include in the intervals. The default is 0.9.
 #' @param intervals For \code{ppc_loo_intervals} and \code{ppc_loo_ribbon},
-#'   optionally a matrix of precomputed LOO predictive intervals intervals with
+#'   optionally a matrix of precomputed LOO predictive intervals
 #'   that can be specified instead of \code{yrep} and \code{lw} (these are both
 #'   ignored if \code{intervals} is specified). If not specified the intervals
 #'   are computed internally before plotting. If specified, \code{intervals}
 #'   must be a matrix with number of rows equal to the number of data points and
-#'   three columns in the following order: the first for the lower bound of the
-#'   interval, the second for median (50\%), and the third for the interval
-#'   upper bound (column names are ignored).
+#'   five columns in the following order: lower outer interval, lower inner
+#'   interval, median (50\%), upper inner interval and upper outer interval
+#'   (column names are ignored).
 #' @param order For \code{ppc_loo_intervals}, a string indicating how to arrange
 #'   the plotted intervals. The default (\code{"index"}) is to plot them in the
 #'   order of the observations. The alternative (\code{"median"}) arranges them
@@ -301,7 +300,8 @@ ppc_loo_intervals <-
            subset = NULL,
            intervals = NULL,
            ...,
-           prob = 0.9,
+           prob = 0.5,
+           prob_outer = 0.9,
            size = 1,
            fatten = 3,
            order = c("index", "median")) {
@@ -310,11 +310,14 @@ ppc_loo_intervals <-
     y <- validate_y(y)
     order_by_median <- match.arg(order) == "median"
     if (!is.null(intervals)) {
-      stopifnot(is.matrix(intervals), ncol(intervals) == 3)
+      stopifnot(is.matrix(intervals), ncol(intervals) %in% c(3, 5))
       message(
         "'intervals' specified so ignoring ",
         "'yrep', 'psis_object', 'subset', if specified."
       )
+      if (ncol(intervals) == 3) {
+          intervals <- cbind(intervals[, 1], intervals, intervals[, 3])
+      }
     } else {
       suggested_package("loo", min_version = "2.0.0")
       yrep <- validate_yrep(yrep, y)
@@ -324,7 +327,8 @@ ppc_loo_intervals <-
         yrep <- yrep[, subset, drop=FALSE]
         psis_object <- .psis_subset(psis_object, subset)
       }
-      a <- (1 - prob) / 2
+      probs <- sort(c(prob, prob_outer))
+      a <- (1 - probs) / 2
       stopifnot(identical(dim(psis_object), dim(yrep)))
       intervals <- suppressWarnings(t(loo::E_loo(
         x = yrep,
@@ -368,17 +372,21 @@ ppc_loo_ribbon <-
            subset = NULL,
            intervals = NULL,
            ...,
-           prob = 0.9,
+           prob = 0.5,
+           prob_outer = 0.9,
            alpha = 0.33,
            size = 0.25) {
     check_ignored_arguments(...)
     y <- validate_y(y)
     if (!is.null(intervals)) {
-      stopifnot(is.matrix(intervals), ncol(intervals) == 3)
+      stopifnot(is.matrix(intervals), ncol(intervals) %in% c(3, 5))
       message(
         "'intervals' specified so ignoring ",
         "'yrep', 'psis_object', 'subset', if specified."
       )
+      if (ncol(intervals) == 3) {
+          intervals <- cbind(intervals[, 1], intervals, intervals[, 3])
+      }
     } else {
       suggested_package("loo", min_version = "2.0.0")
       yrep <- validate_yrep(yrep, y)
@@ -388,7 +396,8 @@ ppc_loo_ribbon <-
         yrep <- yrep[, subset, drop=FALSE]
         psis_object <- .psis_subset(psis_object, subset)
       }
-      a <- (1 - prob) / 2
+      probs <- sort(c(prob, prob_outer))
+      a <- (1 - probs) / 2
       stopifnot(identical(dim(psis_object), dim(yrep)))
       intervals <- suppressWarnings(t(loo::E_loo(
         x = yrep,
@@ -417,9 +426,11 @@ ppc_loo_ribbon <-
     y_id = seq_along(y),
     y_obs = y,
     x = x,
-    lo = intervals[, 1],
-    mid = intervals[, 2],
-    hi = intervals[, 3])
+    ll = intervals[, 1],
+    l  = intervals[, 2],
+    m  = intervals[, 3],
+    h  = intervals[, 4],
+    hh = intervals[, 5])
 }
 
 # subset a psis_object without breaking it

man-roxygen/args-prob-prob_outer.R

@@ -0,0 +1,3 @@
+#' @param prob,prob_outer Values between 0 and 1 indicating the desired
+#'   probability mass to include in the inner and outer intervals. The defaults
+#'   are \code{prob=0.5} and \code{prob_outer=0.9}.