From 7235fa9d3ec2adc328eeb3fef5030853dd3856fc Mon Sep 17 00:00:00 2001 From: Tristan Jay Mahr Date: Mon, 7 Aug 2017 11:10:52 -0500 Subject: [PATCH 1/8] import rlang package --- DESCRIPTION | 3 ++- NAMESPACE | 1 + R/bayesplot-package.R | 2 +- 3 files changed, 4 insertions(+), 2 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 0a54b030..fcd44f59 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -22,7 +22,8 @@ Imports: ggplot2 (>= 2.2.1), reshape2, stats, - utils + utils, + rlang Suggests: arm, gridExtra (>= 2.2.1), diff --git a/NAMESPACE b/NAMESPACE index bd5eafe9..eeb877be 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -114,6 +114,7 @@ export(yaxis_text) export(yaxis_ticks) export(yaxis_title) import(ggplot2) +import(rlang) import(stats) importFrom(dplyr,"%>%") importFrom(dplyr,arrange_) diff --git a/R/bayesplot-package.R b/R/bayesplot-package.R index e43e63ec..5faac322 100644 --- a/R/bayesplot-package.R +++ b/R/bayesplot-package.R @@ -4,7 +4,7 @@ #' @name bayesplot-package #' @aliases bayesplot #' -#' @import ggplot2 stats +#' @import ggplot2 stats rlang #' #' @description #' \if{html}{ From 3ee739154a14ca9cb2153f949e292c1e35e9c0fc Mon Sep 17 00:00:00 2001 From: Tristan Jay Mahr Date: Mon, 7 Aug 2017 11:25:18 -0500 Subject: [PATCH 2/8] don't use is.missing() in validate_x(). use {}'s --- R/helpers-ppc.R | 81 ++++++++++++++++++++++++++++++++++--------------- 1 file changed, 57 insertions(+), 24 deletions(-) diff --git a/R/helpers-ppc.R b/R/helpers-ppc.R index ef57e47a..0ed54d36 100644 --- a/R/helpers-ppc.R +++ b/R/helpers-ppc.R @@ -1,10 +1,13 @@ # Check if an object is a vector (but not list) or a 1-D array is_vector_or_1Darray <- function(x) { - if (is.vector(x) && !is.list(x)) + if (is.vector(x) && !is.list(x)) { return(TRUE) + } + isTRUE(is.array(x) && length(dim(x)) == 1) } + # Validate y # # Checks that y is numeric, doesn't have any NAs, and is either a vector, 1-D @@ -17,13 +20,15 @@ validate_y <- function(y) { stopifnot(is.numeric(y)) if (!(inherits(y, "ts") && is.null(dim(y)))) { - if (!is_vector_or_1Darray(y)) + if (!is_vector_or_1Darray(y)) { stop("'y' must be a vector or 1D array.") + } y <- as.vector(y) } - if (anyNA(y)) + if (anyNA(y)) { stop("NAs not allowed in 'y'.") + } unname(y) } @@ -40,15 +45,21 @@ validate_y <- function(y) { validate_yrep <- function(yrep, y) { stopifnot(is.matrix(yrep), is.numeric(yrep)) if (is.integer(yrep)) { - if (nrow(yrep) == 1) + if (nrow(yrep) == 1) { yrep[1, ] <- as.numeric(yrep[1,, drop = FALSE]) - else + } + else { yrep <- apply(yrep, 2, as.numeric) + } } - if (anyNA(yrep)) + + if (anyNA(yrep)) { stop("NAs not allowed in 'yrep'.") - if (ncol(yrep) != length(y)) + } + + if (ncol(yrep) != length(y)) { stop("ncol(yrep) must be equal to length(y).") + } unclass(unname(yrep)) } @@ -64,18 +75,27 @@ validate_yrep <- function(yrep, y) { # validate_group <- function(group, y) { stopifnot(is.vector(group) || is.factor(group)) - if (!is.factor(group)) + + if (!is.factor(group)) { group <- as.factor(group) - if (anyNA(group)) + } + + if (anyNA(group)) { stop("NAs not allowed in 'group'.") - if (length(group) != length(y)) + } + + if (length(group) != length(y)) { stop("length(group) must be equal to length(y).") - if (length(unique(group)) == 1) + } + + if (length(unique(group)) == 1) { stop("'group' must have more than one unique value.") + } unname(group) } + # Validate x # # Checks that x is a numeric vector, doesn't have any NAs, and has the @@ -84,25 +104,33 @@ validate_group <- function(group, y) { # @param x,y The user's x vector and the y object returned by validate_y. # @param unique_x T/F indicating whether to require all unique values in x. # @return Either throws an error or returns a numeric vector. -# -validate_x <- function(x, y, unique_x = FALSE) { - if (missing(x)) { +validate_x <- function(x = NULL, y, unique_x = FALSE) { + if (is.null(x)) { if (inherits(y, "ts") && is.null(dim(y))) { - return(stats::time(y)) - } else return(1:length(y)) + x <- stats::time(y) + } else { + x <- seq_along(y) + } } stopifnot(is.numeric(x)) - if (!is_vector_or_1Darray(x)) + + if (!is_vector_or_1Darray(x)) { stop("'x' must be a vector or 1D array.") + } x <- as.vector(x) - if (length(x) != length(y)) + if (length(x) != length(y)) { stop("length(x) must be equal to length(y).") - if (anyNA(x)) + } + + if (anyNA(x)) { stop("NAs not allowed in 'x'.") - if (unique_x) + } + + if (unique_x) { stopifnot(identical(length(x), length(unique(x)))) + } unname(x) } @@ -129,6 +157,7 @@ melt_yrep <- function(yrep, label = TRUE) { out } + # Stack y below melted yrep data # # @param y Validated y input. @@ -152,6 +181,7 @@ melt_and_stack <- function(y, yrep, label = TRUE) { }) } + # Prepare data for use in PPCs by group # # @param y,yrep,group Validated y, yrep, and group objects from the user. @@ -169,21 +199,24 @@ ppc_group_data <- function(y, yrep, group, stat = NULL) { colnames(d) <- gsub(".", "_", colnames(d), fixed = TRUE) molten_d <- reshape2::melt(d, id.vars = "group") molten_d <- dplyr::group_by_(molten_d, .dots = list(~group, ~variable)) - if (is.null(stat)) + if (is.null(stat)) { return(molten_d) + } - if (!is.function(stat)) + if (!is.function(stat)) { stat <- match.fun(stat) + } dplyr::summarise_(molten_d, value = ~stat(value)) } # set mapping depending on freq argument set_hist_aes <- function(freq = TRUE, ...) { - if (freq) + if (freq) { aes_(x = ~ value, ...) - else + } else { aes_(x = ~ value, y = ~ ..density.., ...) + } } # check if x consists of whole numbers (very close to integers) From f43a65d3fc2fe52bd36f24d46f44a6d94ce07edc Mon Sep 17 00:00:00 2001 From: Tristan Jay Mahr Date: Mon, 7 Aug 2017 11:26:46 -0500 Subject: [PATCH 3/8] refactor ppc_intervals_data() and plotting functions to remove `is_y` column. --- R/ppc-intervals.R | 282 +++++++++++++++------------- R/ppc-loo.R | 29 +-- tests/testthat/test-ppc-intervals.R | 6 +- tests/testthat/test-ppc-loo.R | 1 + 4 files changed, 170 insertions(+), 148 deletions(-) diff --git a/R/ppc-intervals.R b/R/ppc-intervals.R index dfd86f3f..c1f20998 100644 --- a/R/ppc-intervals.R +++ b/R/ppc-intervals.R @@ -10,7 +10,7 @@ #' @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 then \code{1:length(y)} is used for the x-axis. +#' 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. @@ -104,21 +104,22 @@ NULL #' @export ppc_intervals <- function(y, yrep, - x, + x = NULL, ..., prob = 0.9, size = 1, fatten = 3) { check_ignored_arguments(...) - y <- validate_y(y) - .ppc_intervals( - data = .ppc_intervals_data( + data <- .ppc_intervals_data( y = y, - yrep = validate_yrep(yrep, y), - x = validate_x(x, y), + yrep = yrep, + x = x, group = NULL, prob = prob - ), + ) + + .ppc_intervals( + data = data, size = size, fatten = fatten, grouped = FALSE, @@ -135,7 +136,7 @@ ppc_intervals <- function(y, #' ppc_intervals_grouped <- function(y, yrep, - x, + x = NULL, group, facet_args = list(), ..., @@ -143,18 +144,21 @@ ppc_intervals_grouped <- function(y, size = 1, fatten = 3) { check_ignored_arguments(...) - y <- validate_y(y) - if (is.null(facet_args[["scales"]])) + + if (is.null(facet_args[["scales"]])) { facet_args[["scales"]] <- "free" + } + + data <- .ppc_intervals_data( + y = y, + yrep = yrep, + x = x, + group = group, + prob = prob + ) .ppc_intervals( - data = .ppc_intervals_data( - y = y, - yrep = validate_yrep(yrep, y), - x = validate_x(x, y), - group = validate_group(group, y), - prob = prob - ), + data = data, facet_args = facet_args, size = size, fatten = fatten, @@ -168,22 +172,22 @@ ppc_intervals_grouped <- function(y, #' @export ppc_ribbon <- function(y, yrep, - x, + x = NULL, ..., prob = 0.9, alpha = 0.33, size = 0.25) { check_ignored_arguments(...) + data <- .ppc_intervals_data( + y = y, + yrep = yrep, + x = x, + group = NULL, + prob = prob + ) - y <- validate_y(y) .ppc_intervals( - data = .ppc_intervals_data( - y = y, - yrep = validate_yrep(yrep, y), - x = validate_x(x, y), - group = NULL, - prob = prob - ), + data = data, alpha = alpha, size = size, grouped = FALSE, @@ -197,36 +201,38 @@ ppc_ribbon <- function(y, ppc_ribbon_grouped <- function(y, yrep, - x, + x = NULL, group, facet_args = list(), ..., prob = 0.9, alpha = 0.33, size = 0.25) { - check_ignored_arguments(...) - - y <- validate_y(y) - if (is.null(facet_args[["scales"]])) - facet_args[["scales"]] <- "free" + check_ignored_arguments(...) - .ppc_intervals( - data = .ppc_intervals_data( - y = y, - yrep = validate_yrep(yrep, y), - x = validate_x(x, y), - group = validate_group(group, y), - prob = prob - ), - facet_args = facet_args, - alpha = alpha, - size = size, - grouped = TRUE, - style = "ribbon", - x_lab = label_x(x) - ) + if (is.null(facet_args[["scales"]])) { + facet_args[["scales"]] <- "free" } + data <- .ppc_intervals_data( + y = y, + yrep = yrep, + x = x, + group = group, + prob = prob + ) + + .ppc_intervals( + data = data, + facet_args = facet_args, + alpha = alpha, + size = size, + grouped = TRUE, + style = "ribbon", + x_lab = label_x(x) + ) +} + # internal ---------------------------------------------------------------- label_x <- function(x) { @@ -236,36 +242,47 @@ label_x <- function(x) { .ppc_intervals_data <- function(y, yrep, - x, + x = NULL, group = NULL, prob = 0.8) { - grouped <- !is.null(group) - stopifnot(prob > 0 && prob < 1) + grouped <- !is.null(group) + stopifnot(prob > 0 && prob < 1) - molten_d <- dplyr::rename_(.data = melt_and_stack(y, yrep), - .dots = setNames(list(~ y_id), "x")) - molten_d <- dplyr::arrange_(.data = molten_d, - .dots = list(~ rep_id, ~ x)) - molten_d$x <- x - if (grouped) { - molten_d$group <- group - dots <- list(~ x, ~ group, ~ is_y) - } else { - dots <- list(~ x, ~ is_y) - } - grouped_d <- dplyr::group_by_(molten_d, .dots = dots) - alpha <- (1 - prob) / 2 - probs <- sort(c(alpha, 0.5, 1 - alpha)) - dplyr::summarise_( - .data = grouped_d, - .dots = list( - lo = ~ quantile(value, prob = probs[1]), - mid = ~ quantile(value, prob = probs[2]), - hi = ~ quantile(value, prob = probs[3]) - ) - ) + y <- validate_y(y) + yrep <- validate_yrep(yrep, y) + x <- validate_x(x, y) + + long_d <- melt_and_stack(y, yrep, label = FALSE) + long_d$x <- x[long_d$y_id] + long_d$y_obs <- y[long_d$y_id] + + molten_reps <- long_d[!as.logical(long_d[["is_y"]]), , drop = FALSE] + molten_reps$is_y <- NULL + + if (grouped) { + group <- validate_group(group, y) + molten_reps$group <- group[molten_reps$y_id] + group_vars <- quos(y_id, y_obs, group, x) + } else { + group_vars <- quos(y_id, y_obs, x) } + grouped_d <- dplyr::group_by(molten_reps, !!! group_vars) + alpha <- (1 - prob) / 2 + probs <- sort(c(alpha, 0.5, 1 - alpha)) + + val_col <- quo(value) + dplyr::ungroup(dplyr::summarise( + grouped_d, + # lo_prob = probs[1], + # mid_prob = probs[2], + # hi_prob = probs[3], + lo = quantile(!! val_col, prob = probs[1]), + mid = quantile(!! val_col, prob = probs[2]), + hi = quantile(!! val_col, prob = probs[3]) + )) +} + # Make intervals or ribbon plot # # @param data The object returned by .ppc_intervals_data @@ -281,75 +298,72 @@ label_x <- function(x) { style = c("intervals", "ribbon"), x_lab = NULL) { - style <- match.arg(style) - data[["is_y"]] <- as.logical(data[["is_y"]]) - yrep_data <- dplyr::filter_(data, ~ !is_y) - y_data <- dplyr::filter_(data, ~ is_y) + style <- match.arg(style) - graph <- ggplot( - data = yrep_data, - mapping = aes_( - x = ~ x, - y = ~ mid, - ymin = ~ lo, - ymax = ~ hi - ) + graph <- ggplot( + data = data, + mapping = aes_( + x = ~ x, + y = ~ mid, + ymin = ~ lo, + ymax = ~ hi ) + ) - if (style == "ribbon") { - graph <- graph + - geom_ribbon( - aes_(color = "yrep", fill = "yrep"), - alpha = alpha, - size = size - ) + - geom_blank(aes_(fill = "y")) + - geom_line( - aes_(color = "yrep"), - size = size/2 - ) + - geom_line( - aes_(color = "y"), - data = y_data, - size = 0.5 - ) - } else { - graph <- graph + - geom_pointrange( - mapping = aes_(color = "yrep", fill = "yrep"), - shape = 21, - size = size, - fatten = fatten - ) + - geom_point( - data = y_data, - mapping = aes_(color = "y", fill = "y"), - shape = 21, - size = 1.5 - ) - } + if (style == "ribbon") { + graph <- graph + + geom_ribbon( + aes_(color = "yrep", fill = "yrep"), + alpha = alpha, + size = size + ) + + geom_line( + aes_(color = "yrep"), + size = size/2 + ) + + geom_blank(aes_(fill = "y")) + + geom_line( + aes_(y = ~ y_obs, color = "y"), + size = 0.5 + ) + } else { graph <- graph + - scale_color_manual( - name = "", - values = setNames(get_color(c("lh", "dh")), c("yrep", "y")), - labels = c(yrep = yrep_label(), y = y_label()) + geom_pointrange( + mapping = aes_(color = "yrep", fill = "yrep"), + shape = 21, + size = size, + fatten = fatten ) + - scale_fill_manual( - name = "", - values = c(yrep = get_color("l"), - y = if (style == "ribbon") NA else get_color("d")), - labels = c(yrep = yrep_label(), y = y_label()) + geom_point( + mapping = aes_(y = ~ y_obs, color = "y", fill = "y"), + shape = 21, + size = 1.5 ) + } + + graph <- graph + + scale_color_manual( + name = "", + values = setNames(get_color(c("lh", "dh")), c("yrep", "y")), + labels = c(yrep = yrep_label(), y = y_label()) + ) + + scale_fill_manual( + name = "", + values = c(yrep = get_color("l"), + y = if (style == "ribbon") NA else get_color("d")), + labels = c(yrep = yrep_label(), y = y_label()) + ) - if (grouped) { - facet_args[["facets"]] <- "group" - if (is.null(facet_args[["scales"]])) - facet_args[["scales"]] <- "free" - graph <- graph + - do.call("facet_wrap", facet_args) + if (grouped) { + facet_args[["facets"]] <- "group" + if (is.null(facet_args[["scales"]])) { + facet_args[["scales"]] <- "free" } - graph + - labs(y = NULL, x = x_lab %||% expression(italic(x))) + graph <- graph + do.call("facet_wrap", facet_args) } + graph + + labs(y = NULL, x = x_lab %||% expression(italic(x))) +} + diff --git a/R/ppc-loo.R b/R/ppc-loo.R index e9ada85b..5ca5fc5a 100644 --- a/R/ppc-loo.R +++ b/R/ppc-loo.R @@ -172,16 +172,17 @@ ppc_loo_intervals <- function(y, yrep, lw, - intervals, + intervals = NULL, ..., prob = 0.9, size = 1, fatten = 3, order = c("index", "median")) { + check_ignored_arguments(...) y <- validate_y(y) order_by_median <- match.arg(order) == "median" - if (!missing(intervals)) { + if (!is.null(intervals)) { stopifnot(is.matrix(intervals), ncol(intervals) == 3) message("'intervals' specified so ignoring 'yrep' and 'lw' if specified.") } else { @@ -198,8 +199,9 @@ ppc_loo_intervals <- } x <- seq_along(y) - if (order_by_median) + if (order_by_median) { x <- reorder(x, intervals[, 2]) + } graph <- .ppc_intervals( data = .loo_intervals_data(y, x, intervals), @@ -209,8 +211,10 @@ ppc_loo_intervals <- fatten = fatten, x_lab = "Data point (index)" ) - if (!order_by_median) + + if (!order_by_median) { return(graph) + } graph + xlab("Ordered by median") + @@ -224,14 +228,14 @@ ppc_loo_ribbon <- function(y, yrep, lw, - intervals, + intervals = NULL, ..., prob = 0.9, alpha = 0.33, size = 0.25) { check_ignored_arguments(...) y <- validate_y(y) - if (!missing(intervals)) { + if (!is.null(intervals)) { stopifnot(is.matrix(intervals), ncol(intervals) == 3) message("'intervals' specified so ignoring 'yrep' and 'lw' if specified.") } else { @@ -260,11 +264,14 @@ ppc_loo_ribbon <- # internal ---------------------------------------------------------------- .loo_intervals_data <- function(y, x, intervals) { - colnames(intervals) <- c("lo", "mid", "hi") stopifnot(length(y) == nrow(intervals), length(x) == length(y)) - dplyr::bind_rows( - data.frame(x, is_y = TRUE, lo = y, mid = y, hi = y), - data.frame(x, is_y = FALSE, intervals) - ) + + data.frame( + y_id = seq_along(y), + y_obs = y, + x = x, + lo = intervals[, 1], + mid = intervals[, 2], + hi = intervals[, 3]) } diff --git a/tests/testthat/test-ppc-intervals.R b/tests/testthat/test-ppc-intervals.R index 7467b642..b648b74b 100644 --- a/tests/testthat/test-ppc-intervals.R +++ b/tests/testthat/test-ppc-intervals.R @@ -1,7 +1,7 @@ library(bayesplot) context("PPC: intervals & ribbon") -source("data-for-ppc-tests.R") +source(test_path("data-for-ppc-tests.R")) test_that("ppc_intervals returns ggplot object", { expect_gg(ppc_intervals(y, yrep)) @@ -31,8 +31,8 @@ test_that("ppc_ribbon_grouped returns ggplot object", { test_that(".ppc_intervals_data returns correct structure", { d <- .ppc_intervals_data(y, yrep, x = 1:length(y)) d_group <- .ppc_intervals_data(y, yrep, x, group) - expect_named(d, c("x", "is_y", "lo", "mid", "hi")) - expect_named(d_group, c("x", "group", "is_y","lo", "mid", "hi")) + expect_named(d, c("y_id", "y_obs", "x", "lo", "mid", "hi")) + expect_named(d_group, c("y_id", "y_obs", "group", "x", "lo", "mid", "hi")) expect_error(.ppc_intervals_data(y, yrep, x = 1:length(y), prob = 0), "prob") expect_error(.ppc_intervals_data(y, yrep, x = 1:length(y), prob = 1.01), "prob") diff --git a/tests/testthat/test-ppc-loo.R b/tests/testthat/test-ppc-loo.R index 06707651..43a2e3e6 100644 --- a/tests/testthat/test-ppc-loo.R +++ b/tests/testthat/test-ppc-loo.R @@ -44,6 +44,7 @@ test_that("ppc_loo_intervals returns ggplot object", { expect_s3_class(g$data$x, "factor") expect_equal(nlevels(g$data$x), length(g$data$x)) }) + test_that("ppc_loo_ribbon returns ggplot object", { expect_gg(ppc_loo_ribbon(y, yrep, lw, prob = 0.7, alpha = 0.1)) }) From e457ebe91c9cfbf26103086a0318831113291e6c Mon Sep 17 00:00:00 2001 From: Tristan Jay Mahr Date: Mon, 7 Aug 2017 14:01:40 -0500 Subject: [PATCH 4/8] document and export ppc_intervals_data and ppc_ribbon_data --- NAMESPACE | 2 ++ R/ppc-intervals.R | 54 +++++++++++++++++++---------- man-roxygen/return-ggplot-or-data.R | 3 ++ man/PPC-intervals.Rd | 25 +++++++++---- man/PPC-loo.Rd | 8 ++--- 5 files changed, 63 insertions(+), 29 deletions(-) create mode 100644 man-roxygen/return-ggplot-or-data.R diff --git a/NAMESPACE b/NAMESPACE index eeb877be..12ef385b 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -88,11 +88,13 @@ export(ppc_freqpoly) export(ppc_freqpoly_grouped) export(ppc_hist) export(ppc_intervals) +export(ppc_intervals_data) export(ppc_intervals_grouped) export(ppc_loo_intervals) export(ppc_loo_pit) export(ppc_loo_ribbon) export(ppc_ribbon) +export(ppc_ribbon_data) export(ppc_ribbon_grouped) export(ppc_rootogram) export(ppc_scatter) diff --git a/R/ppc-intervals.R b/R/ppc-intervals.R index c1f20998..5ab9a1e8 100644 --- a/R/ppc-intervals.R +++ b/R/ppc-intervals.R @@ -19,7 +19,7 @@ #' \code{\link[ggplot2]{geom_ribbon}}. For interval plots \code{size} and #' \code{fatten} are passed to \code{\link[ggplot2]{geom_pointrange}}. #' -#' @template return-ggplot +#' @template return-ggplot-or-data #' #' @templateVar bdaRef (Ch. 6) #' @template reference-bda @@ -73,6 +73,9 @@ #' xaxis_ticks(FALSE) + #' panel_bg(fill = "gray20") #' +#' ppc_data <- ppc_intervals_data(y, yrep, x = year, prob = 0.5) +#' ppc_group_data <- ppc_intervals_data(y, yrep, x = year, group, prob = 0.5) +#' #' \dontrun{ #' library("rstanarm") #' fit <- stan_glmer(mpg ~ wt + (1|cyl), data = mtcars) @@ -110,7 +113,8 @@ ppc_intervals <- function(y, size = 1, fatten = 3) { check_ignored_arguments(...) - data <- .ppc_intervals_data( + + data <- ppc_intervals_data( y = y, yrep = yrep, x = x, @@ -149,7 +153,7 @@ ppc_intervals_grouped <- function(y, facet_args[["scales"]] <- "free" } - data <- .ppc_intervals_data( + data <- ppc_intervals_data( y = y, yrep = yrep, x = x, @@ -168,6 +172,7 @@ ppc_intervals_grouped <- function(y, ) } + #' @rdname PPC-intervals #' @export ppc_ribbon <- function(y, @@ -178,7 +183,8 @@ ppc_ribbon <- function(y, alpha = 0.33, size = 0.25) { check_ignored_arguments(...) - data <- .ppc_intervals_data( + + data <- ppc_intervals_data( y = y, yrep = yrep, x = x, @@ -196,25 +202,25 @@ ppc_ribbon <- function(y, ) } + #' @export #' @rdname PPC-intervals -ppc_ribbon_grouped <- - function(y, - yrep, - x = NULL, - group, - facet_args = list(), - ..., - prob = 0.9, - alpha = 0.33, - size = 0.25) { +ppc_ribbon_grouped <- function(y, + yrep, + x = NULL, + group, + facet_args = list(), + ..., + prob = 0.9, + alpha = 0.33, + size = 0.25) { check_ignored_arguments(...) if (is.null(facet_args[["scales"]])) { facet_args[["scales"]] <- "free" } - data <- .ppc_intervals_data( + data <- ppc_intervals_data( y = y, yrep = yrep, x = x, @@ -234,6 +240,21 @@ ppc_ribbon_grouped <- } +#' @rdname PPC-intervals +#' @export +ppc_intervals_data <- function(y, yrep, x = NULL, group = NULL, prob = 0.9, ...) { + check_ignored_arguments(...) + .ppc_intervals_data(y = y, yrep = yrep, x = x, group = group, prob = prob) +} + + +#' @rdname PPC-intervals +#' @export +ppc_ribbon_data <- ppc_intervals_data + + + + # internal ---------------------------------------------------------------- label_x <- function(x) { if (missing(x)) "Index" else NULL @@ -274,9 +295,6 @@ label_x <- function(x) { val_col <- quo(value) dplyr::ungroup(dplyr::summarise( grouped_d, - # lo_prob = probs[1], - # mid_prob = probs[2], - # hi_prob = probs[3], lo = quantile(!! val_col, prob = probs[1]), mid = quantile(!! val_col, prob = probs[2]), hi = quantile(!! val_col, prob = probs[3]) diff --git a/man-roxygen/return-ggplot-or-data.R b/man-roxygen/return-ggplot-or-data.R new file mode 100644 index 00000000..6909c6e3 --- /dev/null +++ b/man-roxygen/return-ggplot-or-data.R @@ -0,0 +1,3 @@ +#' @return A ggplot object that can be further customized using the +#' \pkg{ggplot2} package. The \code{_data} functions return the data that +#' would have be drawn by the plotting function. diff --git a/man/PPC-intervals.Rd b/man/PPC-intervals.Rd index 2648c8ba..6d7c7a8d 100644 --- a/man/PPC-intervals.Rd +++ b/man/PPC-intervals.Rd @@ -6,17 +6,24 @@ \alias{ppc_intervals_grouped} \alias{ppc_ribbon} \alias{ppc_ribbon_grouped} +\alias{ppc_intervals_data} +\alias{ppc_ribbon_data} \title{PPC intervals} \usage{ -ppc_intervals(y, yrep, x, ..., prob = 0.9, size = 1, fatten = 3) +ppc_intervals(y, yrep, x = NULL, ..., prob = 0.9, size = 1, fatten = 3) -ppc_intervals_grouped(y, yrep, x, group, facet_args = list(), ..., +ppc_intervals_grouped(y, yrep, x = NULL, group, facet_args = list(), ..., prob = 0.9, size = 1, fatten = 3) -ppc_ribbon(y, yrep, x, ..., prob = 0.9, alpha = 0.33, size = 0.25) +ppc_ribbon(y, yrep, x = NULL, ..., prob = 0.9, alpha = 0.33, + size = 0.25) -ppc_ribbon_grouped(y, yrep, x, group, facet_args = list(), ..., prob = 0.9, - alpha = 0.33, size = 0.25) +ppc_ribbon_grouped(y, yrep, x = NULL, group, facet_args = list(), ..., + prob = 0.9, alpha = 0.33, size = 0.25) + +ppc_intervals_data(y, yrep, x = NULL, group = NULL, prob = 0.9, ...) + +ppc_ribbon_data(y, yrep, x = NULL, group = NULL, prob = 0.9, ...) } \arguments{ \item{y}{A vector of observations. See \strong{Details}.} @@ -32,7 +39,7 @@ instructions.} \item{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 then \code{1:length(y)} is used for the x-axis.} +is missing or NULL, then \code{1:length(y)} is used for the x-axis.} \item{...}{Currently unused.} @@ -53,7 +60,8 @@ passed to \code{\link[ggplot2]{facet_wrap}} to control faceting.} } \value{ A ggplot object that can be further customized using the - \pkg{ggplot2} package. + \pkg{ggplot2} package. The \code{_data} functions return the data that + would have be drawn by the plotting function. } \description{ Medians and central interval estimates of \code{yrep} with \code{y} overlaid. @@ -110,6 +118,9 @@ ppc_ribbon_grouped( xaxis_ticks(FALSE) + panel_bg(fill = "gray20") +ppc_data <- ppc_intervals_data(y, yrep, x = year, prob = 0.5) +ppc_group_data <- ppc_intervals_data(y, yrep, x = year, group, prob = 0.5) + \dontrun{ library("rstanarm") fit <- stan_glmer(mpg ~ wt + (1|cyl), data = mtcars) diff --git a/man/PPC-loo.Rd b/man/PPC-loo.Rd index 871a9d3e..fcae88e3 100644 --- a/man/PPC-loo.Rd +++ b/man/PPC-loo.Rd @@ -10,11 +10,11 @@ ppc_loo_pit(y, yrep, lw, pit, compare = c("uniform", "normal"), ..., size = 2, alpha = 1) -ppc_loo_intervals(y, yrep, lw, intervals, ..., prob = 0.9, size = 1, - fatten = 3, order = c("index", "median")) +ppc_loo_intervals(y, yrep, lw, intervals = NULL, ..., prob = 0.9, + size = 1, fatten = 3, order = c("index", "median")) -ppc_loo_ribbon(y, yrep, lw, intervals, ..., prob = 0.9, alpha = 0.33, - size = 0.25) +ppc_loo_ribbon(y, yrep, lw, intervals = NULL, ..., prob = 0.9, + alpha = 0.33, size = 0.25) } \arguments{ \item{y}{A vector of observations. See \strong{Details}.} From a0bcab377af25a9c372f8f8c8c8941cb536603b4 Mon Sep 17 00:00:00 2001 From: Tristan Jay Mahr Date: Mon, 7 Aug 2017 14:47:51 -0500 Subject: [PATCH 5/8] fix global variable warnings from R CMD check --- R/ppc-intervals.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/ppc-intervals.R b/R/ppc-intervals.R index 5ab9a1e8..5eb78fc8 100644 --- a/R/ppc-intervals.R +++ b/R/ppc-intervals.R @@ -283,16 +283,16 @@ label_x <- function(x) { if (grouped) { group <- validate_group(group, y) molten_reps$group <- group[molten_reps$y_id] - group_vars <- quos(y_id, y_obs, group, x) + group_vars <- quos(!!! syms(c("y_id", "y_obs", "group", "x"))) } else { - group_vars <- quos(y_id, y_obs, x) + group_vars <- quos(!!! syms(c("y_id", "y_obs", "x"))) } grouped_d <- dplyr::group_by(molten_reps, !!! group_vars) alpha <- (1 - prob) / 2 probs <- sort(c(alpha, 0.5, 1 - alpha)) - val_col <- quo(value) + val_col <- quo(!! sym("value")) dplyr::ungroup(dplyr::summarise( grouped_d, lo = quantile(!! val_col, prob = probs[1]), From 1b1485d06b4901fd42b51778c4fb61e32db681a7 Mon Sep 17 00:00:00 2001 From: Tristan Jay Mahr Date: Mon, 7 Aug 2017 15:12:37 -0500 Subject: [PATCH 6/8] simplify standard evaluation code --- R/ppc-intervals.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/ppc-intervals.R b/R/ppc-intervals.R index 5eb78fc8..14a05824 100644 --- a/R/ppc-intervals.R +++ b/R/ppc-intervals.R @@ -283,16 +283,16 @@ label_x <- function(x) { if (grouped) { group <- validate_group(group, y) molten_reps$group <- group[molten_reps$y_id] - group_vars <- quos(!!! syms(c("y_id", "y_obs", "group", "x"))) + group_vars <- syms(c("y_id", "y_obs", "group", "x")) } else { - group_vars <- quos(!!! syms(c("y_id", "y_obs", "x"))) + group_vars <- syms(c("y_id", "y_obs", "x")) } grouped_d <- dplyr::group_by(molten_reps, !!! group_vars) alpha <- (1 - prob) / 2 probs <- sort(c(alpha, 0.5, 1 - alpha)) - val_col <- quo(!! sym("value")) + val_col <- sym("value") dplyr::ungroup(dplyr::summarise( grouped_d, lo = quantile(!! val_col, prob = probs[1]), From 271a668952125b5589c9ded7c256c0550084cc1c Mon Sep 17 00:00:00 2001 From: Tristan Jay Mahr Date: Tue, 8 Aug 2017 08:07:59 -0500 Subject: [PATCH 7/8] require newer dplyr version --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 82e53edc..3b140216 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -19,7 +19,7 @@ SystemRequirements: pandoc Depends: R (>= 3.1.0) Imports: - dplyr (>= 0.4.3), + dplyr (>= 0.7.1), ggplot2 (>= 2.2.1), reshape2, stats, From 972e29098f0993e4d990ca2358a465e17b8f544e Mon Sep 17 00:00:00 2001 From: Tristan Jay Mahr Date: Tue, 8 Aug 2017 08:30:08 -0500 Subject: [PATCH 8/8] have ppc_intervals_data return prob value. improve tests. --- R/ppc-intervals.R | 1 + tests/testthat/test-ppc-intervals.R | 39 +++++++++++++++++++++++------ 2 files changed, 32 insertions(+), 8 deletions(-) diff --git a/R/ppc-intervals.R b/R/ppc-intervals.R index 14a05824..6a1b3da0 100644 --- a/R/ppc-intervals.R +++ b/R/ppc-intervals.R @@ -295,6 +295,7 @@ label_x <- function(x) { 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]) diff --git a/tests/testthat/test-ppc-intervals.R b/tests/testthat/test-ppc-intervals.R index b648b74b..fbf353c0 100644 --- a/tests/testthat/test-ppc-intervals.R +++ b/tests/testthat/test-ppc-intervals.R @@ -28,13 +28,36 @@ test_that("ppc_ribbon_grouped returns ggplot object", { expect_gg(ppc_ribbon_grouped(y, yrep, x, group, facet_args = list(scales = "fixed"))) }) -test_that(".ppc_intervals_data returns correct structure", { - d <- .ppc_intervals_data(y, yrep, x = 1:length(y)) - d_group <- .ppc_intervals_data(y, yrep, x, group) - expect_named(d, c("y_id", "y_obs", "x", "lo", "mid", "hi")) - expect_named(d_group, c("y_id", "y_obs", "group", "x", "lo", "mid", "hi")) - - expect_error(.ppc_intervals_data(y, yrep, x = 1:length(y), prob = 0), "prob") - expect_error(.ppc_intervals_data(y, yrep, x = 1:length(y), prob = 1.01), "prob") +test_that("ppc_intervals_data returns correct structure", { + d <- ppc_intervals_data(y, yrep, x = 1:length(y), prob = .9) + d_group <- ppc_intervals_data(y, yrep, x, group) + expect_named(d, c("y_id", "y_obs", "x", + "prob", "lo", "mid", "hi")) + expect_named(d_group, c("y_id", "y_obs", "group", "x", + "prob", "lo", "mid", "hi")) + + expect_error(ppc_intervals_data(y, yrep, x = 1:length(y), prob = 0), "prob") + expect_error(ppc_intervals_data(y, yrep, x = 1:length(y), prob = 1.01), "prob") +}) + +test_that("ppc_intervals_data does math correctly", { + d <- ppc_intervals_data(y, yrep, prob = .9) + qs <- unname(quantile(yrep[, 1], c(.05, .5, .95))) + expect_equal(d$lo[1], qs[1]) + expect_equal(d$mid[1], qs[2]) + expect_equal(d$hi[1], qs[3]) + + # Testing groups and known quantiles + y <- rep(10, 4) + group <- c("a", "a", "b", "b") + yrep_g1 <- matrix(rep((0:20), 2), ncol = 2) + yrep_g2 <- yrep_g1 - 10 + yrep <- cbind(yrep_g1, yrep_g2) + + d <- ppc_intervals_data(y, yrep, group = group, prob = .9) + expect_equal(unique(d$prob), .9) + expect_equal(d$lo, c( 1, 1, -9, -9)) + expect_equal(d$mid, c(10, 10, 0, 0)) + expect_equal(d$hi, c(19, 19, 9, 9)) })