diff --git a/DESCRIPTION b/DESCRIPTION index 8f4faf38..00752c14 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -35,7 +35,7 @@ Imports: ggridges (>= 0.5.5), glue, lifecycle, - posterior, + posterior (>= 1.7.0), reshape2, rlang (>= 1.0.0), stats, @@ -44,6 +44,7 @@ Imports: tidyselect, utils Suggests: + brms (>= 2.23.0), ggdist, ggfortify, gridExtra (>= 2.2.1), @@ -59,7 +60,8 @@ Suggests: shinystan (>= 2.3.0), survival, testthat (>= 3.0.0), - vdiffr (>= 1.0.2) + vdiffr (>= 1.0.2), + patchwork RoxygenNote: 7.3.3 VignetteBuilder: knitr Encoding: UTF-8 diff --git a/R/helpers-ppc.R b/R/helpers-ppc.R index 618b9fca..192ce565 100644 --- a/R/helpers-ppc.R +++ b/R/helpers-ppc.R @@ -590,3 +590,8 @@ create_rep_ids <- function(ids) paste('italic(y)[rep] (', ids, ")") y_label <- function() expression(italic(y)) yrep_label <- function() expression(italic(y)[rep]) ypred_label <- function() expression(italic(y)[pred]) + +# helper function for formatting p-value when displayed in a plot +fmt_p <- function(x) { + dplyr::if_else(x < 0.0005, "0.000", as.character(round(signif(x, 2) + 1e-10, 3))) +} diff --git a/R/ppc-distributions.R b/R/ppc-distributions.R index d32a2e1a..6db2ed8d 100644 --- a/R/ppc-distributions.R +++ b/R/ppc-distributions.R @@ -657,6 +657,44 @@ ppc_violin_grouped <- #' @param pit An optional vector of probability integral transformed values for #' which the ECDF is to be drawn. If NULL, PIT values are computed to `y` with #' respect to the corresponding values in `yrep`. +#' @param interpolate_adj For `ppc_pit_ecdf()` and `ppc_pit_ecdf_grouped()` +#' when `method = "independent"`, +#' a boolean defining if the simultaneous confidence bands should be +#' interpolated based on precomputed values rather than computed exactly. +#' Computing the bands may be computationally intensive and the approximation +#' gives a fast method for assessing the ECDF trajectory. The default is to use +#' interpolation if `K` is greater than 200. +#' @param method For `ppc_pit_ecdf()` and `ppc_pit_ecdf_grouped()`, the method +#' used to calculate the +#' uniformity test: +#' * `"independent"`: (default) Assumes independence (Säilynoja et al., 2022). +#' * `"correlated"`: Accounts for correlation (Tesso & Vehtari, 2026). +#' @param test For `ppc_pit_ecdf()` and `ppc_pit_ecdf_grouped()` when +#' `method = "correlated"`, which +#' dependence-aware test to use: `"POT"`, `"PRIT"`, or `"PIET"`. +#' Defaults to `"POT"`. +#' @param gamma For `ppc_pit_ecdf()` and `ppc_pit_ecdf_grouped()` when +#' `method = "correlated"`, tolerance +#' threshold controlling how strongly suspicious points are flagged. Larger +#' values (gamma > 0) emphasizes points with larger deviations. If `NULL`, automatically +#' determined based on p-value. +#' @param linewidth For `ppc_pit_ecdf()` and `ppc_pit_ecdf_grouped()` when +#' `method = "correlated"`, the line width of the ECDF and highlighting +#' points. Defaults to 0.3. +#' @param color For `ppc_pit_ecdf()` and `ppc_pit_ecdf_grouped()` when +#' `method = "correlated"`, a vector +#' with base color and highlight color for the ECDF plot. Defaults to +#' `c(ecdf = "grey60", highlight = "red")`. The first element is used for +#' the main ECDF line, the second for highlighted suspicious regions. +#' @param help_text For `ppc_pit_ecdf()` and `ppc_pit_ecdf_grouped()` when +#' `method = "correlated"`, a boolean defining whether to add informative +#' text to the plot. Defaults to `TRUE`. +#' @param pareto_pit For `ppc_pit_ecdf()` and `ppc_pit_ecdf_grouped()`, a +#' boolean defining whether to compute the PIT values using Pareto-smoothed +#' importance sampling (if `TRUE` and no pit values are provided). +#' Defaults to `TRUE` when `method = "correlated"` and `test` is `"POT"` or `"PIET"`. +#' Otherwise defaults to `FALSE`. If `TRUE` requires the specification of `lw` or `psis_object`. +#' The defaults should not be changed by the user, but the option is provided for developers. #' @rdname PPC-distributions #' ppc_pit_ecdf <- function(y, @@ -666,59 +704,258 @@ ppc_pit_ecdf <- function(y, K = NULL, prob = .99, plot_diff = FALSE, - interpolate_adj = NULL) { + interpolate_adj = NULL, + method = NULL, + test = NULL, + gamma = NULL, + linewidth = NULL, + color = NULL, + help_text = NULL, + pareto_pit = NULL + ) { check_ignored_arguments(..., - ok_args = c("K", "pit", "prob", "plot_diff", "interpolate_adj") + ok_args = c( + "K", "pareto_pit", "pit", "prob", "plot_diff", + "interpolate_adj", "method", "test", "gamma", "linewidth", + "color", "help_text" + ) ) - if (is.null(pit)) { - pit <- ppc_data(y, yrep) %>% - group_by(.data$y_id) %>% + .warn_ignored <- function(method_name, args) { + inform(paste0( + "As method = ", method_name, " specified; ignoring: ", + paste(args, collapse = ", "), "." + )) + } + + if (is.null(method)) { + inform(c( + "i" = paste( + "In the next major release, the default `method`", + "will change to 'correlated'." + ), + "*" = paste( + "To silence this message, explicitly set", + "`method = 'independent'` or `method = 'correlated'`." + ) + )) + method <- "independent" + } else { + method <- match.arg(method, choices = c("independent", "correlated")) + if (method == "independent") { + inform("The 'independent' method is superseded by the 'correlated' method.") + } + } + + switch(method, + "correlated" = { + if (!is.null(interpolate_adj)) .warn_ignored("'correlated'", "interpolate_adj") + + test <- match.arg(test %||% "POT", choices = c("POT", "PRIT", "PIET")) + alpha <- 1 - prob + gamma <- gamma %||% 0 + linewidth <- linewidth %||% 0.3 + color <- color %||% c(ecdf = "grey60", highlight = "red") + help_text <- help_text %||% TRUE + pareto_pit <- pareto_pit %||% is.null(pit) && test %in% c("POT", "PIET") + }, + "independent" = { + ignored <- c( + if (!is.null(test)) "test", + if (!is.null(gamma)) "gamma", + if (!is.null(help_text)) "help_text" + ) + if (length(ignored) > 0) .warn_ignored("'independent'", ignored) + pareto_pit <- pareto_pit %||% FALSE + } + ) + + if (isTRUE(pareto_pit) && is.null(pit)) { + # --- Pareto-smoothed PIT --- + suggested_package("rstantools") + y <- validate_y(y) + yrep <- validate_predictions(yrep, length(y)) + pit <- posterior::pareto_pit(x = yrep, y = y, weights = NULL, log = TRUE) + K <- K %||% length(pit) + + } else if (!is.null(pit)) { + # --- Pre-supplied PIT values --- + pit <- validate_pit(pit) + K <- K %||% length(pit) + ignored <- c( + if (!missing(y) && !is.null(y)) "y", + if (!missing(yrep) && !is.null(yrep)) "yrep" + ) + if (length(ignored) > 0) { + inform(paste0( + "As 'pit' specified; ignoring: ", + paste(ignored, collapse = ", "), "." + )) + } + } else { + # --- Empirical PIT ---' + y <- validate_y(y) + yrep <- validate_predictions(yrep, length(y)) + pit <- ppc_data(y, yrep) |> + group_by(.data$y_id) |> dplyr::group_map( ~ mean(.x$value[.x$is_y] > .x$value[!.x$is_y]) + runif(1, max = mean(.x$value[.x$is_y] == .x$value[!.x$is_y])) - ) %>% + ) |> unlist() - if (is.null(K)) { - K <- min(nrow(yrep) + 1, 1000) + K <- K %||% min(nrow(yrep) + 1, 1000) + } + + n_obs <- length(pit) + unit_interval <- seq(0, 1, length.out = K) + ecdf_pit_fn <- ecdf(pit) + y_label <- if (plot_diff) "ECDF difference" else "ECDF" + + if (method == "correlated") { + test_res <- posterior::uniformity_test(pit = pit, test = test) + p_value_CCT <- test_res$pvalue + pointwise_contrib <- test_res$pointwise + max_contrib <- max(pointwise_contrib) + if (gamma < 0 || gamma > max_contrib) { + stop(sprintf( + "gamma must be in [0, %.2f], but gamma = %s was provided.", + max_contrib, gamma + )) } - } else { - inform("'pit' specified so ignoring 'y', and 'yrep' if specified.") - pit <- validate_pit(pit) - if (is.null(K)) { - K <- length(pit) + x_combined <- sort(unique(c(unit_interval, pit))) + df_main <- tibble::tibble( + x = x_combined, + ecdf_val = ecdf_pit_fn(x_combined) - plot_diff * x_combined + ) + pit_sorted <- sort(pit) + df_pit <- tibble::tibble( + pit = pit_sorted, + ecdf_val = ecdf_pit_fn(pit_sorted) - plot_diff * pit_sorted + ) + + p <- ggplot() + + geom_step( + data = df_main, + mapping = aes(x = .data$x, y = .data$ecdf_val), + show.legend = FALSE, + linewidth = linewidth, + color = color["ecdf"] + ) + + geom_segment( + mapping = aes(x = 0, y = 0, xend = 1, yend = if (plot_diff) 0 else 1), + linetype = "dashed", + color = "darkgrey", + linewidth = 0.3 + ) + + labs(x = "PIT", y = y_label) + + if (p_value_CCT < alpha) { + red_idx <- which(pointwise_contrib > gamma) + + if (length(red_idx) > 0) { + df_red <- df_pit[red_idx, ] + df_red$segment <- cumsum(c(1, diff(red_idx) != 1)) + seg_sizes <- stats::ave(df_red$pit, df_red$segment, FUN = length) + df_isolated <- df_red[seg_sizes == 1, ] + df_grouped <- df_red[seg_sizes > 1, ] + + if (nrow(df_grouped) > 0) { + df_segments <- do.call(rbind, lapply( + split(df_grouped, df_grouped$segment), + function(grp) { + pit_idx <- match(grp$pit, x_combined) + idx_range <- seq(min(pit_idx), max(pit_idx)) + tibble::tibble( + x = df_main$x[idx_range], + ecdf_val = df_main$ecdf_val[idx_range], + segment = grp$segment[1L] + ) + } + )) + + p <- p + geom_step( + data = df_segments, + mapping = aes(x = .data$x, y = .data$ecdf_val, group = .data$segment), + color = color["highlight"], + linewidth = linewidth + 0.8 + ) + } + + if (nrow(df_isolated) > 0) { + p <- p + geom_point( + data = df_isolated, + aes(x = .data$pit, y = .data$ecdf_val), + color = color["highlight"], + size = linewidth + 1 + ) + } + } } + + if (isTRUE(help_text)) { + label_size <- 0.7 * bayesplot_theme_get()$text@size / ggplot2::.pt + p <- p + annotate( + "text", + x = -Inf, y = Inf, + label = sprintf("p[unif]^{%s} == '%s' ~ (alpha == '%.2f')", + test, fmt_p(p_value_CCT), alpha + ), + hjust = -0.05, + vjust = 1.5, + color = "black", + parse = TRUE, + size = label_size + ) + } + + if (plot_diff) { + epsilon = max( + sqrt(log(2 / (1 - prob)) / (2 * n_obs)), + max(abs(df_main$ecdf_val)) + ) + p <- p + scale_y_continuous(limits = c(-epsilon, epsilon)) + } + + p <- p + + yaxis_ticks(FALSE) + + scale_color_ppc() + + bayesplot_theme_get() + + return(p) } - N <- length(pit) - gamma <- adjust_gamma( - N = N, - K = K, - prob = prob, - interpolate_adj = interpolate_adj + + gamma_indep <- adjust_gamma( + N = n_obs, K = K, prob = prob, interpolate_adj = interpolate_adj ) - lims <- ecdf_intervals(gamma = gamma, N = N, K = K) - ggplot() + - aes( - x = seq(0,1,length.out = K), - y = ecdf(pit)(seq(0, 1, length.out = K)) - - (plot_diff == TRUE) * seq(0, 1, length.out = K), - color = "y" + lims <- ecdf_intervals(gamma = gamma_indep, N = n_obs, K = K) + lims_upper <- lims$upper[-1] / n_obs - plot_diff * unit_interval + lims_lower <- lims$lower[-1] / n_obs - plot_diff * unit_interval + ecdf_eval <- ecdf_pit_fn(unit_interval) - plot_diff * unit_interval + + p <- ggplot() + + geom_step( + mapping = aes(x = unit_interval, y = lims_upper, color = "yrep"), + linetype = "dashed", + linewidth = 0.3, + show.legend = FALSE + ) + + geom_step( + mapping = aes(x = unit_interval, y = lims_lower, color = "yrep"), + linetype = "dashed", + linewidth = 0.3, + show.legend = FALSE ) + - geom_step(show.legend = FALSE) + - geom_step(aes( - y = lims$upper[-1] / N - (plot_diff == TRUE) * seq(0, 1, length.out = K), - color = "yrep" - ), - linetype = 2, show.legend = FALSE) + - geom_step(aes( - y = lims$lower[-1] / N - (plot_diff == TRUE) * seq(0, 1, length.out = K), - color = "yrep" - ), - linetype = 2, show.legend = FALSE) + - labs(y = ifelse(plot_diff,"ECDF - difference","ECDF"), x = "PIT") + + geom_step( + mapping = aes(x = unit_interval, y = ecdf_eval, color = "y"), + linewidth = 0.5, + show.legend = FALSE + ) + + labs(x = "PIT", y = y_label) + yaxis_ticks(FALSE) + scale_color_ppc() + bayesplot_theme_get() + + return(p) } #' @export @@ -733,56 +970,300 @@ ppc_pit_ecdf_grouped <- pit = NULL, prob = .99, plot_diff = FALSE, - interpolate_adj = NULL) { + interpolate_adj = NULL, + method = NULL, + test = NULL, + gamma = NULL, + linewidth = NULL, + color = NULL, + help_text = NULL, + pareto_pit = NULL) { check_ignored_arguments(..., - ok_args = c("K", "pit", "prob", "plot_diff", "interpolate_adj") + ok_args = c("K", "pareto_pit", "pit", "prob", "plot_diff", + "interpolate_adj", "method", "test", "gamma", + "linewidth", "color", "help_text") ) - if (is.null(pit)) { + .warn_ignored <- function(method_name, args) { + inform(paste0( + "As method = ", method_name, " specified; ignoring: ", + paste(args, collapse = ", "), "." + )) + } + + # Resolve and validate `method` + if (is.null(method)) { + inform(c( + "i" = paste( + "In the next major release, the default `method`", + "will change to 'correlated'." + ), + "*" = paste( + "To silence this message, explicitly set", + "`method = 'independent'` or `method = 'correlated'`." + ) + )) + method <- "independent" + } else { + method <- match.arg(method, choices = c("independent", "correlated")) + if (method == "independent") { + inform("The 'independent' method is superseded by the 'correlated' method.") + } + } + + switch(method, + "correlated" = { + if (!is.null(interpolate_adj)) .warn_ignored("'correlated'", "interpolate_adj") + test <- match.arg(test %||% "POT", choices = c("POT", "PRIT", "PIET")) + alpha <- 1 - prob + gamma <- gamma %||% 0 + linewidth <- linewidth %||% 0.3 + color <- color %||% c(ecdf = "grey60", highlight = "red") + help_text <- help_text %||% TRUE + pareto_pit <- pareto_pit %||% is.null(pit) && test %in% c("POT", "PIET") + }, + "independent" = { + ignored <- c( + if (!is.null(test)) "test", + if (!is.null(gamma)) "gamma", + if (!is.null(help_text)) "help_text" + ) + if (length(ignored) > 0) .warn_ignored("'independent'", ignored) + pareto_pit <- pareto_pit %||% FALSE + } + ) + + if (isTRUE(pareto_pit) && is.null(pit)) { + suggested_package("rstantools") + y <- validate_y(y) + yrep <- validate_predictions(yrep, length(y)) + group <- validate_group(group, length(y)) + pit <- posterior::pareto_pit(x = yrep, y = y, weights = NULL, log = TRUE) + K <- K %||% length(pit) + } else if (!is.null(pit)) { + pit <- validate_pit(pit) + group <- validate_group(group, length(pit)) + K <- K %||% length(pit) + ignored <- c( + if (!missing(y) && !is.null(y)) "y", + if (!missing(yrep) && !is.null(yrep)) "yrep" + ) + if (length(ignored) > 0) { + inform(paste0( + "As 'pit' specified; ignoring: ", + paste(ignored, collapse = ", "), "." + )) + } + } else { + y <- validate_y(y) + yrep <- validate_predictions(yrep, length(y)) + group <- validate_group(group, length(y)) pit <- ppc_data(y, yrep, group) %>% group_by(.data$y_id) %>% dplyr::group_map( ~ mean(.x$value[.x$is_y] > .x$value[!.x$is_y]) + runif(1, max = mean(.x$value[.x$is_y] == .x$value[!.x$is_y])) - ) %>% + ) %>% unlist() - if (is.null(K)) { - K <- min(nrow(yrep) + 1, 1000) + K <- K %||% min(nrow(yrep) + 1, 1000) + } + + data <- data.frame(pit = pit, group = group, stringsAsFactors = FALSE) + group_levels <- unique(data$group) + + if (method == "correlated") { + data_cor <- dplyr::group_by(data, .data$group) %>% + dplyr::group_map(function(.x, .y) { + n_obs <- nrow(.x) + K_g <- K %||% n_obs + unit_interval <- seq(0, 1, length.out = K_g) + ecdf_pit_fn <- ecdf(.x$pit) + x_combined <- sort(unique(c(unit_interval, .x$pit))) + df_main <- data.frame( + x = x_combined, + ecdf_value = ecdf_pit_fn(x_combined) - plot_diff * x_combined, + group = .y[[1]], + stringsAsFactors = FALSE + ) + + test_res <- posterior::uniformity_test(pit = .x$pit, test = test) + p_value_CCT <- test_res$pvalue + pointwise_contrib <- test_res$pointwise + max_contrib <- max(pointwise_contrib) + if (gamma < 0 || gamma > max_contrib) { + stop(sprintf( + "gamma must be in [0, %.2f], but gamma = %s was provided.", + max_contrib, gamma + )) + } + + red <- NULL + red_points <- NULL + if (p_value_CCT < alpha) { + red_idx <- which(pointwise_contrib > gamma) + if (length(red_idx) > 0) { + pit_sorted <- sort(.x$pit) + df_pit <- data.frame( + pit = pit_sorted, + ecdf_value = ecdf_pit_fn(pit_sorted), + stringsAsFactors = FALSE + ) + df_red <- df_pit[red_idx, , drop = FALSE] + df_red$segment <- cumsum(c(1, diff(red_idx) != 1)) + seg_sizes <- stats::ave(df_red$pit, df_red$segment, FUN = length) + df_isolated <- df_red[seg_sizes == 1, , drop = FALSE] + df_grouped <- df_red[seg_sizes > 1, , drop = FALSE] + + if (nrow(df_grouped) > 0) { + red <- do.call(rbind, lapply( + split(df_grouped, df_grouped$segment), + function(grp) { + pit_idx <- match(grp$pit, x_combined) + idx_range <- seq(min(pit_idx), max(pit_idx)) + data.frame( + x = x_combined[idx_range], + ecdf_value = ecdf_pit_fn(x_combined[idx_range]) - + plot_diff * x_combined[idx_range], + segment = grp$segment[1], + group = .y[[1]], + stringsAsFactors = FALSE + ) + } + )) + } + + if (nrow(df_isolated) > 0) { + red_points <- data.frame( + x = df_isolated$pit, + ecdf_value = df_isolated$ecdf_value - plot_diff * df_isolated$pit, + group = .y[[1]], + stringsAsFactors = FALSE + ) + } + } + } + + ann <- NULL + if (isTRUE(help_text)) { + ann <- data.frame( + group = .y[[1]], + x = -Inf, + y = Inf, + label = sprintf( + "p[unif]^{%s} == '%s' ~ (alpha == '%.2f')", + test, fmt_p(p_value_CCT), alpha + ), + stringsAsFactors = FALSE + ) + } + + list(main = df_main, red = red, red_points = red_points, ann = ann) + }) + + main_df <- dplyr::bind_rows(lapply(data_cor, `[[`, "main")) + red_df <- dplyr::bind_rows(lapply(data_cor, `[[`, "red")) + red_points_df <- dplyr::bind_rows(lapply(data_cor, `[[`, "red_points")) + ann_df <- dplyr::bind_rows(lapply(data_cor, `[[`, "ann")) + ref_df <- data.frame( + group = group_levels, + x = 0, + y = 0, + xend = 1, + yend = if (plot_diff) 0 else 1, + stringsAsFactors = FALSE + ) + + p <- ggplot() + + geom_step( + data = main_df, + mapping = aes(x = .data$x, y = .data$ecdf_value, group = .data$group), + show.legend = FALSE, + linewidth = linewidth, + color = color["ecdf"] + ) + + geom_segment( + data = ref_df, + mapping = aes( + x = .data$x, + y = .data$y, + xend = .data$xend, + yend = .data$yend + ), + linetype = "dashed", + color = "darkgrey", + linewidth = 0.3 + ) + + if (nrow(red_df) > 0) { + p <- p + geom_step( + data = red_df, + mapping = aes(x = .data$x, y = .data$ecdf_value, group = interaction(.data$group, .data$segment)), + color = color["highlight"], + linewidth = linewidth + 0.8 + ) } - } else { - inform("'pit' specified so ignoring 'y' and 'yrep' if specified.") - pit <- validate_pit(pit) + + if (nrow(red_points_df) > 0) { + p <- p + geom_point( + data = red_points_df, + mapping = aes(x = .data$x, y = .data$ecdf_value), + color = color["highlight"], + size = linewidth + 1 + ) + } + + if (isTRUE(help_text) && nrow(ann_df) > 0) { + label_size <- 0.7 * bayesplot_theme_get()$text@size / ggplot2::.pt + p <- p + geom_text( + data = ann_df, + mapping = aes(x = .data$x, y = .data$y, label = .data$label), + hjust = -0.05, + vjust = 1.5, + color = "black", + parse = TRUE, + size = label_size + ) + } + + return( + p + + labs(y = if (plot_diff) "ECDF difference" else "ECDF", x = "PIT") + + yaxis_ticks(FALSE) + + bayesplot_theme_get() + + facet_wrap("group") + + scale_color_ppc() + + force_axes_in_facets() + ) } - N <- length(pit) - gammas <- lapply(unique(group), function(g) { - N_g <- sum(group == g) + gammas <- lapply(group_levels, function(g) { + N_g <- sum(data$group == g) adjust_gamma( N = N_g, - K = ifelse(is.null(K), N_g, K), + K = K %||% N_g, prob = prob, interpolate_adj = interpolate_adj ) }) - names(gammas) <- unique(group) + names(gammas) <- group_levels - data <- data.frame(pit = pit, group = group) %>% - group_by(group) %>% + data <- data %>% + dplyr::group_by(.data$group) %>% dplyr::group_map( ~ data.frame( - ecdf_value = ecdf(.x$pit)(seq(0, 1, length.out = ifelse(is.null(K), nrow(.x), K))), + ecdf_value = ecdf(.x$pit)(seq(0, 1, length.out = K %||% nrow(.x))), group = .y[1], lims_upper = ecdf_intervals( gamma = gammas[[unlist(.y[1])]], N = nrow(.x), - K = ifelse(is.null(K), nrow(.x), K) + K = K %||% nrow(.x) )$upper[-1] / nrow(.x), lims_lower = ecdf_intervals( gamma = gammas[[unlist(.y[1])]], N = nrow(.x), - K = ifelse(is.null(K), nrow(.x), K) + K = K %||% nrow(.x) )$lower[-1] / nrow(.x), - x = seq(0, 1, length.out = ifelse(is.null(K), nrow(.x), K)) + x = seq(0, 1, length.out = K %||% nrow(.x)) ) ) %>% dplyr::bind_rows() diff --git a/R/ppc-loo.R b/R/ppc-loo.R index 2f91c5b5..07cb9a76 100644 --- a/R/ppc-loo.R +++ b/R/ppc-loo.R @@ -21,7 +21,9 @@ #' [ggplot2::geom_density()], respectively. For `ppc_loo_intervals()`, `size` #' `linewidth` and `fatten` are passed to [ggplot2::geom_pointrange()]. For #' `ppc_loo_ribbon()`, `alpha` and `size` are passed to -#' [ggplot2::geom_ribbon()]. +#' [ggplot2::geom_ribbon()]. For `ppc_loo_pit_ecdf()`, linewidth for the ECDF plot. When +#' `method = "correlated"`, defaults to 0.3. When `method = "independent"`, +#' if `NULL` no linewidth is specified for the ECDF line. #' #' @template return-ggplot-or-data #' @@ -404,12 +406,37 @@ ppc_loo_pit_qq <- function(y, #' expectation for uniform PIT values rather than plotting the regular ECDF. #' The default is `FALSE`, but for large samples we recommend setting #' `plot_diff = TRUE` to better use the plot area. -#' @param interpolate_adj For `ppc_loo_pit_ecdf()`, a boolean defining if the -#' simultaneous confidence bands should be interpolated based on precomputed -#' values rather than computed exactly. Computing the bands may be -#' computationally intensive and the approximation gives a fast method for -#' assessing the ECDF trajectory. The default is to use interpolation if `K` -#' is greater than 200. +#' @param interpolate_adj For `ppc_loo_pit_ecdf()` when `method = "independent"`, +#' a boolean defining if the simultaneous confidence bands should be +#' interpolated based on precomputed values rather than computed exactly. +#' Computing the bands may be computationally intensive and the approximation +#' gives a fast method for assessing the ECDF trajectory. The default is to use +#' interpolation if `K` is greater than 200. +#' @param method For `ppc_loo_pit_ecdf()`, the method used to calculate the +#' uniformity test: +#' * `"independent"`: (Current default) Assumes independence (Säilynoja et al., 2022). +#' * `"correlated"`: (Recommended) Accounts for correlation (Tesso & Vehtari, 2026). +#' @param test For `ppc_loo_pit_ecdf()` when `method = "correlated"`, which +#' dependence-aware test to use: `"POT"`, `"PRIT"`, or `"PIET"`. +#' Defaults to `"POT"`. +#' @param gamma For `ppc_loo_pit_ecdf()` when `method = "correlated"`, tolerance +#' threshold controlling how strongly suspicious points are flagged. Larger +#' values (gamma > 0) emphasizes points with larger deviations. If `NULL`, automatically +#' determined based on p-value. +#' @param color For `ppc_loo_pit_ecdf()` when `method = "correlated"`, a vector +#' with base color and highlight color for the ECDF plot. Defaults to +#' `c(ecdf = "grey60", highlight = "red")`. The first element is used for +#' the main ECDF line, the second for highlighted suspicious regions. +#' @param help_text For `ppc_loo_pit_ecdf()` when `method = "correlated"`, a boolean +#' defining whether to add informative text to the plot. Defaults to `TRUE`. +#' @param pareto_pit For `ppc_loo_pit_ecdf()`. Computes PIT values using Pareto-PIT method. +#' Defaults to `TRUE` if `test` is either `"POT"` or `"PIET"` and no `pit` values are +#' provided otherwise `FALSE`. This argument should normally not be modified by the user, +#' except for development purposes. +#' @note +#' Note that the default "independent" method is **superseded** by +#' the "correlated" method (Tesso & Vehtari, 2026) which accounts for dependent +#' LOO-PIT values. ppc_loo_pit_ecdf <- function(y, yrep, lw = NULL, @@ -419,63 +446,252 @@ ppc_loo_pit_ecdf <- function(y, K = NULL, prob = .99, plot_diff = FALSE, - interpolate_adj = NULL) { + interpolate_adj = NULL, + method = NULL, + test = NULL, + gamma = NULL, + linewidth = NULL, + color = NULL, + help_text = NULL, + pareto_pit = NULL) { check_ignored_arguments(..., ok_args = list("moment_match")) - if (!is.null(pit)) { - inform("'pit' specified so ignoring 'y','yrep','lw' if specified.") + .warn_ignored <- function(method_name, args) { + inform(paste0( + "As method = ", method_name, " specified; ignoring: ", + paste(args, collapse = ", "), "." + )) + } + + if (is.null(method)) { + inform(c( + "i" = paste( + "In the next major release, the default `method`", + "will change to 'correlated'." + ), + "*" = paste( + "To silence this message, explicitly set", + "`method = 'independent'` or `method = 'correlated'`." + ) + )) + method <- "independent" + } else { + method <- match.arg(method, choices = c("independent", "correlated")) + if (method == "independent") { + inform("The 'independent' method is superseded by the 'correlated' method.") + } + } + + switch(method, + "correlated" = { + if (!is.null(interpolate_adj)) .warn_ignored("'correlated'", "interpolate_adj") + + test <- match.arg(test %||% "POT", choices = c("POT", "PRIT", "PIET")) + alpha <- 1 - prob + gamma <- gamma %||% 0 + linewidth <- linewidth %||% 0.3 + color <- color %||% c(ecdf = "grey60", highlight = "red") + help_text <- help_text %||% TRUE + pareto_pit <- pareto_pit %||% is.null(pit) && test %in% c("POT", "PIET") + }, + "independent" = { + # Collect args that are meaningless under the independent method. + ignored <- c( + if (!is.null(test)) "test", + if (!is.null(gamma)) "gamma", + if (!is.null(help_text)) "help_text" + ) + if (length(ignored) > 0) .warn_ignored("'independent'", ignored) + } + ) + + if (isTRUE(pareto_pit) && is.null(pit)) { + # --- Pareto-smoothed LOO PIT --- + suggested_package("rstantools") + y <- validate_y(y) + yrep <- validate_predictions(yrep, length(y)) + lw <- .get_lw(lw, psis_object) + stopifnot(identical(dim(yrep), dim(lw))) + pit <- posterior::pareto_pit(x = yrep, y = y, weights = lw, log = TRUE) + K <- K %||% length(pit) + } else if (!is.null(pit)) { + # --- Pre-supplied PIT values --- pit <- validate_pit(pit) - if (is.null(K)) { - K <- length(pit) + K <- K %||% length(pit) + ignored <- c( + if (!missing(y) && !is.null(y)) "y", + if (!missing(yrep) && !is.null(yrep)) "yrep", + if (!is.null(lw)) "lw" + ) + if (length(ignored) > 0) { + inform(paste0( + "As 'pit' specified; ignoring: ", + paste(ignored, collapse = ", "), "." + )) } } else { + # --- Standard LOO PIT --- suggested_package("rstantools") y <- validate_y(y) yrep <- validate_predictions(yrep, length(y)) lw <- .get_lw(lw, psis_object) stopifnot(identical(dim(yrep), dim(lw))) pit <- pmin(1, rstantools::loo_pit(object = yrep, y = y, lw = lw)) - if (is.null(K)) { - K <- min(nrow(yrep) + 1, 1000) - } + K <- K %||% min(nrow(yrep) + 1, 1000) } n_obs <- length(pit) - gamma <- adjust_gamma( - N = n_obs, - K = K, - prob = prob, - interpolate_adj = interpolate_adj + unit_interval <- seq(0, 1, length.out = K) + ecdf_pit_fn <- ecdf(pit) + y_label <- if (plot_diff) "ECDF difference" else "ECDF" + + if (method == "correlated") { + test_res <- posterior::uniformity_test(pit = pit, test = test) + p_value_CCT <- test_res$pvalue + pointwise_contrib <- test_res$pointwise + max_contrib <- max(pointwise_contrib) + if (gamma < 0 || gamma > max_contrib) { + stop(sprintf( + "gamma must be in [0, %.2f], but gamma = %s was provided.", + max_contrib, gamma + )) + } + + x_combined <- sort(unique(c(unit_interval, pit))) + df_main <- tibble::tibble( + x = x_combined, + ecdf_val = ecdf_pit_fn(x_combined) - plot_diff * x_combined + ) + pit_sorted <- sort(pit) + df_pit <- tibble::tibble( + pit = pit_sorted, + ecdf_val = ecdf_pit_fn(pit_sorted) - plot_diff * pit_sorted + ) + + p <- ggplot() + + geom_step( + data = df_main, + mapping = aes(x = .data$x, y = .data$ecdf_val), + show.legend = FALSE, + linewidth = linewidth, + color = color["ecdf"] + ) + + geom_segment( + mapping = aes(x = 0, y = 0, xend = 1, yend = if (plot_diff) 0 else 1), + linetype = "dashed", + color = "darkgrey", + linewidth = 0.3 + ) + + labs(x = "LOO-PIT", y = y_label) + + if (p_value_CCT < alpha) { + red_idx <- which(pointwise_contrib > gamma) + + if (length(red_idx) > 0) { + df_red <- df_pit[red_idx, ] + df_red$segment <- cumsum(c(1, diff(red_idx) != 1)) + seg_sizes <- stats::ave(df_red$pit, df_red$segment, FUN = length) + df_isolated <- df_red[seg_sizes == 1, ] + df_grouped <- df_red[seg_sizes > 1, ] + + if (nrow(df_grouped) > 0) { + df_segments <- do.call(rbind, lapply( + split(df_grouped, df_grouped$segment), + function(grp) { + pit_idx <- match(grp$pit, x_combined) + idx_range <- seq(min(pit_idx), max(pit_idx)) + tibble::tibble( + x = df_main$x[idx_range], + ecdf_val = df_main$ecdf_val[idx_range], + segment = grp$segment[1L] + ) + } + )) + + p <- p + geom_step( + data = df_segments, + mapping = aes(x = .data$x, y = .data$ecdf_val, group = .data$segment), + color = color["highlight"], + linewidth = linewidth + 0.8 + ) + } + + if (nrow(df_isolated) > 0) { + p <- p + geom_point( + data = df_isolated, + mapping = aes(x = .data$pit, y = .data$ecdf_val), + color = color["highlight"], + size = linewidth + 1 + ) + } + } + } + + if (isTRUE(help_text)) { + label_size <- 0.7 * bayesplot_theme_get()$text@size / ggplot2::.pt + p <- p + annotate( + "text", + x = -Inf, y = Inf, + label = sprintf( + "p[unif]^{%s} == '%s' ~ (alpha == '%.2f')", + test, fmt_p(p_value_CCT), alpha + ), + hjust = -0.05, + vjust = 1.5, + color = "black", + parse = TRUE, + size = label_size + ) + } + + if (plot_diff) { + epsilon <- max( + sqrt(log(2 / (1 - prob)) / (2 * n_obs)), + max(abs(df_main$ecdf_val)) + ) + p <- p + scale_y_continuous(limits = c(-epsilon, epsilon)) + } + + p <- p + + yaxis_ticks(FALSE) + + scale_color_ppc() + + bayesplot_theme_get() + + return(p) + } + + gamma_indep <- adjust_gamma( + N = n_obs, K = K, prob = prob, interpolate_adj = interpolate_adj ) - lims <- ecdf_intervals(gamma = gamma, N = n_obs, K = K) - ggplot() + - aes( - x = seq(0, 1, length.out = K), - y = ecdf(pit)(seq(0, 1, length.out = K)) - - (plot_diff == TRUE) * seq(0, 1, length.out = K), - color = "y" + lims <- ecdf_intervals(gamma = gamma_indep, N = n_obs, K = K) + lims_upper <- lims$upper[-1L] / n_obs - plot_diff * unit_interval + lims_lower <- lims$lower[-1L] / n_obs - plot_diff * unit_interval + ecdf_eval <- ecdf_pit_fn(unit_interval) - plot_diff * unit_interval + + p <- ggplot() + + geom_step( + mapping = aes(x = unit_interval, y = lims_upper, color = "yrep"), + linetype = "dashed", + linewidth = 0.3, + show.legend = FALSE ) + - geom_step(show.legend = FALSE) + geom_step( - aes( - y = lims$upper[-1] / n_obs - - (plot_diff == TRUE) * seq(0, 1, length.out = K), - color = "yrep" - ), - linetype = 2, show.legend = FALSE + mapping = aes(x = unit_interval, y = lims_lower, color = "yrep"), + linetype = "dashed", + linewidth = 0.3, + show.legend = FALSE ) + geom_step( - aes( - y = lims$lower[-1] / n_obs - - (plot_diff == TRUE) * seq(0, 1, length.out = K), - color = "yrep" - ), - linetype = 2, show.legend = FALSE + mapping = aes(x = unit_interval, y = ecdf_eval, color = "y"), + linewidth = 0.5, + show.legend = FALSE ) + - labs(y = ifelse(plot_diff, "ECDF difference", "ECDF"), x = "LOO PIT") + + labs(x = "LOO-PIT", y = y_label) + yaxis_ticks(FALSE) + scale_color_ppc() + bayesplot_theme_get() + + return(p) } @@ -846,7 +1062,7 @@ ppc_loo_ribbon <- bc_mat <- matrix(0, nrow(unifs), ncol(unifs)) # Generate boundary corrected reference values - for (i in 1:nrow(unifs)) { + for (i in seq_len(nrow(unifs))) { bc_list <- .kde_correction(unifs[i, ], bw = bw, grid_len = grid_len diff --git a/man-roxygen/args-methods-dots.R b/man-roxygen/args-methods-dots.R new file mode 100644 index 00000000..4e0e6c4a --- /dev/null +++ b/man-roxygen/args-methods-dots.R @@ -0,0 +1 @@ +#' @param ... Arguments passed to individual methods (if applicable). diff --git a/man/PPC-distributions.Rd b/man/PPC-distributions.Rd index b0099f2d..2a373b60 100644 --- a/man/PPC-distributions.Rd +++ b/man/PPC-distributions.Rd @@ -131,7 +131,14 @@ ppc_pit_ecdf( K = NULL, prob = 0.99, plot_diff = FALSE, - interpolate_adj = NULL + interpolate_adj = NULL, + method = NULL, + test = NULL, + gamma = NULL, + linewidth = NULL, + color = NULL, + help_text = NULL, + pareto_pit = NULL ) ppc_pit_ecdf_grouped( @@ -143,7 +150,14 @@ ppc_pit_ecdf_grouped( pit = NULL, prob = 0.99, plot_diff = FALSE, - interpolate_adj = NULL + interpolate_adj = NULL, + method = NULL, + test = NULL, + gamma = NULL, + linewidth = NULL, + color = NULL, + help_text = NULL, + pareto_pit = NULL ) } \arguments{ @@ -239,11 +253,53 @@ values rather than plotting the regular ECDF. The default is \code{FALSE}, but for large samples we recommend setting \code{plot_diff=TRUE} as the difference plot will visually show a more dynamic range.} -\item{interpolate_adj}{A boolean defining if the simultaneous confidence -bands should be interpolated based on precomputed values rather than -computed exactly. Computing the bands may be computationally intensive and -the approximation gives a fast method for assessing the ECDF trajectory. -The default is to use interpolation if \code{K} is greater than 200.} +\item{interpolate_adj}{For \code{ppc_pit_ecdf()} and \code{ppc_pit_ecdf_grouped()} +when \code{method = "independent"}, +a boolean defining if the simultaneous confidence bands should be +interpolated based on precomputed values rather than computed exactly. +Computing the bands may be computationally intensive and the approximation +gives a fast method for assessing the ECDF trajectory. The default is to use +interpolation if \code{K} is greater than 200.} + +\item{method}{For \code{ppc_pit_ecdf()} and \code{ppc_pit_ecdf_grouped()}, the method +used to calculate the +uniformity test: +\itemize{ +\item \code{"independent"}: (default) Assumes independence (Säilynoja et al., 2022). +\item \code{"correlated"}: Accounts for correlation (Tesso & Vehtari, 2026). +}} + +\item{test}{For \code{ppc_pit_ecdf()} and \code{ppc_pit_ecdf_grouped()} when +\code{method = "correlated"}, which +dependence-aware test to use: \code{"POT"}, \code{"PRIT"}, or \code{"PIET"}. +Defaults to \code{"POT"}.} + +\item{gamma}{For \code{ppc_pit_ecdf()} and \code{ppc_pit_ecdf_grouped()} when +\code{method = "correlated"}, tolerance +threshold controlling how strongly suspicious points are flagged. Larger +values (gamma > 0) emphasizes points with larger deviations. If \code{NULL}, automatically +determined based on p-value.} + +\item{linewidth}{For \code{ppc_pit_ecdf()} and \code{ppc_pit_ecdf_grouped()} when +\code{method = "correlated"}, the line width of the ECDF and highlighting +points. Defaults to 0.3.} + +\item{color}{For \code{ppc_pit_ecdf()} and \code{ppc_pit_ecdf_grouped()} when +\code{method = "correlated"}, a vector +with base color and highlight color for the ECDF plot. Defaults to +\code{c(ecdf = "grey60", highlight = "red")}. The first element is used for +the main ECDF line, the second for highlighted suspicious regions.} + +\item{help_text}{For \code{ppc_pit_ecdf()} and \code{ppc_pit_ecdf_grouped()} when +\code{method = "correlated"}, a boolean defining whether to add informative +text to the plot. Defaults to \code{TRUE}.} + +\item{pareto_pit}{For \code{ppc_pit_ecdf()} and \code{ppc_pit_ecdf_grouped()}, a +boolean defining whether to compute the PIT values using Pareto-smoothed +importance sampling (if \code{TRUE} and no pit values are provided). +Defaults to \code{TRUE} when \code{method = "correlated"} and \code{test} is \code{"POT"} or \code{"PIET"}. +Otherwise defaults to \code{FALSE}. If \code{TRUE} requires the specification of \code{lw} or \code{psis_object}. +The defaults should not be changed by the user, but the option is provided for developers.} } \value{ The plotting functions return a ggplot object that can be further diff --git a/man/PPC-loo.Rd b/man/PPC-loo.Rd index 7e9ebe02..0656c250 100644 --- a/man/PPC-loo.Rd +++ b/man/PPC-loo.Rd @@ -65,7 +65,14 @@ ppc_loo_pit_ecdf( K = NULL, prob = 0.99, plot_diff = FALSE, - interpolate_adj = NULL + interpolate_adj = NULL, + method = NULL, + test = NULL, + gamma = NULL, + linewidth = NULL, + color = NULL, + help_text = NULL, + pareto_pit = NULL ) ppc_loo_pit( @@ -148,7 +155,9 @@ and \code{alpha} are passed to \code{\link[ggplot2:geom_point]{ggplot2::geom_poi \code{\link[ggplot2:geom_density]{ggplot2::geom_density()}}, respectively. For \code{ppc_loo_intervals()}, \code{size} \code{linewidth} and \code{fatten} are passed to \code{\link[ggplot2:geom_linerange]{ggplot2::geom_pointrange()}}. For \code{ppc_loo_ribbon()}, \code{alpha} and \code{size} are passed to -\code{\link[ggplot2:geom_ribbon]{ggplot2::geom_ribbon()}}.} +\code{\link[ggplot2:geom_ribbon]{ggplot2::geom_ribbon()}}. For \code{ppc_loo_pit_ecdf()}, linewidth for the ECDF plot. When +\code{method = "correlated"}, defaults to 0.3. When \code{method = "independent"}, +if \code{NULL} no linewidth is specified for the ECDF line.} \item{boundary_correction}{For \code{ppc_loo_pit_overlay()}, when set to \code{TRUE} (the default) the function will compute boundary corrected density values @@ -192,12 +201,41 @@ expectation for uniform PIT values rather than plotting the regular ECDF. The default is \code{FALSE}, but for large samples we recommend setting \code{plot_diff = TRUE} to better use the plot area.} -\item{interpolate_adj}{For \code{ppc_loo_pit_ecdf()}, a boolean defining if the -simultaneous confidence bands should be interpolated based on precomputed -values rather than computed exactly. Computing the bands may be -computationally intensive and the approximation gives a fast method for -assessing the ECDF trajectory. The default is to use interpolation if \code{K} -is greater than 200.} +\item{interpolate_adj}{For \code{ppc_loo_pit_ecdf()} when \code{method = "independent"}, +a boolean defining if the simultaneous confidence bands should be +interpolated based on precomputed values rather than computed exactly. +Computing the bands may be computationally intensive and the approximation +gives a fast method for assessing the ECDF trajectory. The default is to use +interpolation if \code{K} is greater than 200.} + +\item{method}{For \code{ppc_loo_pit_ecdf()}, the method used to calculate the +uniformity test: +\itemize{ +\item \code{"independent"}: (Current default) Assumes independence (Säilynoja et al., 2022). +\item \code{"correlated"}: (Recommended) Accounts for correlation (Tesso & Vehtari, 2026). +}} + +\item{test}{For \code{ppc_loo_pit_ecdf()} when \code{method = "correlated"}, which +dependence-aware test to use: \code{"POT"}, \code{"PRIT"}, or \code{"PIET"}. +Defaults to \code{"POT"}.} + +\item{gamma}{For \code{ppc_loo_pit_ecdf()} when \code{method = "correlated"}, tolerance +threshold controlling how strongly suspicious points are flagged. Larger +values (gamma > 0) emphasizes points with larger deviations. If \code{NULL}, automatically +determined based on p-value.} + +\item{color}{For \code{ppc_loo_pit_ecdf()} when \code{method = "correlated"}, a vector +with base color and highlight color for the ECDF plot. Defaults to +\code{c(ecdf = "grey60", highlight = "red")}. The first element is used for +the main ECDF line, the second for highlighted suspicious regions.} + +\item{help_text}{For \code{ppc_loo_pit_ecdf()} when \code{method = "correlated"}, a boolean +defining whether to add informative text to the plot. Defaults to \code{TRUE}.} + +\item{pareto_pit}{For \code{ppc_loo_pit_ecdf()}. Computes PIT values using Pareto-PIT method. +Defaults to \code{TRUE} if \code{test} is either \code{"POT"} or \code{"PIET"} and no \code{pit} values are +provided otherwise \code{FALSE}. This argument should normally not be modified by the user, +except for development purposes.} \item{subset}{For \code{ppc_loo_intervals()} and \code{ppc_loo_ribbon()}, an optional integer vector indicating which observations in \code{y} (and \code{yrep}) to @@ -233,6 +271,11 @@ Leave-One-Out (LOO) predictive checks. See the \strong{Plot Descriptions} sectio below, and \href{https://github.com/jgabry/bayes-vis-paper#readme}{Gabry et al. (2019)} for details. } +\note{ +Note that the default "independent" method is \strong{superseded} by +the "correlated" method (Tesso & Vehtari, 2026) which accounts for dependent +LOO-PIT values. +} \section{Plot Descriptions}{ \describe{ diff --git a/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-color-change.svg b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-color-change.svg new file mode 100644 index 00000000..18b5da8f --- /dev/null +++ b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-color-change.svg @@ -0,0 +1,78 @@ + + + + + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +O +T += +0.000 + +( +α += +0.01 +) + + + +0.00 +0.25 +0.50 +0.75 +1.00 + + + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +PIT +ECDF +ppc_pit_ecdf (color change) + + diff --git a/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-correlated-diff.svg b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-correlated-diff.svg new file mode 100644 index 00000000..9b25f423 --- /dev/null +++ b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-correlated-diff.svg @@ -0,0 +1,71 @@ + + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +O +T += +0.000 + +( +α += +0.01 +) + + + +-0.1 +0.0 +0.1 + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +PIT +ECDF difference +ppc_pit_ecdf (correlated diff) + + diff --git a/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-correlated-piet-2.svg b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-correlated-piet-2.svg new file mode 100644 index 00000000..5be35332 --- /dev/null +++ b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-correlated-piet-2.svg @@ -0,0 +1,75 @@ + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +I +E +T += +0.19 + +( +α += +0.01 +) + + + +0.00 +0.25 +0.50 +0.75 +1.00 + + + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +PIT +ECDF +ppc_pit_ecdf (correlated piet 2) + + diff --git a/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-correlated-piet.svg b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-correlated-piet.svg new file mode 100644 index 00000000..3b47aa76 --- /dev/null +++ b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-correlated-piet.svg @@ -0,0 +1,76 @@ + + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +I +E +T += +0.000 + +( +α += +0.01 +) + + + +0.00 +0.25 +0.50 +0.75 +1.00 + + + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +PIT +ECDF +ppc_pit_ecdf (correlated PIET) + + diff --git a/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-correlated-pot.svg b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-correlated-pot.svg new file mode 100644 index 00000000..41a31cbe --- /dev/null +++ b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-correlated-pot.svg @@ -0,0 +1,78 @@ + + + + + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +O +T += +0.000 + +( +α += +0.01 +) + + + +0.00 +0.25 +0.50 +0.75 +1.00 + + + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +PIT +ECDF +ppc_pit_ecdf (correlated pot) + + diff --git a/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-correlated-prit-2.svg b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-correlated-prit-2.svg new file mode 100644 index 00000000..3e3bda99 --- /dev/null +++ b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-correlated-prit-2.svg @@ -0,0 +1,81 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +R +I +T += +0.001 + +( +α += +0.01 +) + + + +0.00 +0.25 +0.50 +0.75 +1.00 + + + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +PIT +ECDF +ppc_pit_ecdf (correlated prit 2) + + diff --git a/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-correlated-prit.svg b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-correlated-prit.svg new file mode 100644 index 00000000..12d47b39 --- /dev/null +++ b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-correlated-prit.svg @@ -0,0 +1,76 @@ + + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +R +I +T += +0.000 + +( +α += +0.01 +) + + + +0.00 +0.25 +0.50 +0.75 +1.00 + + + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +PIT +ECDF +ppc_pit_ecdf (correlated PRIT) + + diff --git a/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-correlated.svg b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-correlated.svg new file mode 100644 index 00000000..980e8bd7 --- /dev/null +++ b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-correlated.svg @@ -0,0 +1,75 @@ + + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +O +T += +0.000 + +( +α += +0.01 +) + + + +0.00 +0.25 +0.50 +0.75 +1.00 + + + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +PIT +ECDF +ppc_pit_ecdf (correlated) + + diff --git a/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-default.svg b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-default.svg index 8de07070..8a46a332 100644 --- a/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-default.svg +++ b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-default.svg @@ -25,9 +25,9 @@ - - - + + + diff --git a/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-diff-correlated-piet.svg b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-diff-correlated-piet.svg new file mode 100644 index 00000000..e7ed2ab4 --- /dev/null +++ b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-diff-correlated-piet.svg @@ -0,0 +1,75 @@ + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +I +E +T += +0.19 + +( +α += +0.01 +) + + + +-0.10 +-0.05 +0.00 +0.05 +0.10 + + + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +PIT +ECDF difference +ppc_pit_ecdf (diff, correlated piet) + + diff --git a/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-diff-correlated-pot.svg b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-diff-correlated-pot.svg new file mode 100644 index 00000000..0b7d33b0 --- /dev/null +++ b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-diff-correlated-pot.svg @@ -0,0 +1,78 @@ + + + + + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +O +T += +0.000 + +( +α += +0.01 +) + + + +-0.10 +-0.05 +0.00 +0.05 +0.10 + + + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +PIT +ECDF difference +ppc_pit_ecdf (diff, correlated pot) + + diff --git a/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-diff-correlated-prit.svg b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-diff-correlated-prit.svg new file mode 100644 index 00000000..4f44596c --- /dev/null +++ b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-diff-correlated-prit.svg @@ -0,0 +1,81 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +R +I +T += +0.001 + +( +α += +0.01 +) + + + +-0.10 +-0.05 +0.00 +0.05 +0.10 + + + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +PIT +ECDF difference +ppc_pit_ecdf (diff, correlated prit) + + diff --git a/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-diff.svg b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-diff.svg index 94693daf..36e743b5 100644 --- a/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-diff.svg +++ b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-diff.svg @@ -25,9 +25,9 @@ - - - + + + @@ -51,7 +51,7 @@ 0.75 1.00 PIT -ECDF - difference +ECDF difference ppc_pit_ecdf (diff) diff --git a/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-grouped-correlated-diff.svg b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-grouped-correlated-diff.svg new file mode 100644 index 00000000..7f81cfcd --- /dev/null +++ b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-grouped-correlated-diff.svg @@ -0,0 +1,228 @@ + + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +O +T += +0.004 + +( +α += +0.01 +) + + + + + + + + + + + + + + + +p +u +n +i +f +P +O +T += +0.002 + +( +α += +0.01 +) + + + + + + + + + + + + + +p +u +n +i +f +P +O +T += +0.34 + +( +α += +0.01 +) + + + + + + + + + + + + + +p +u +n +i +f +P +O +T += +0.34 + +( +α += +0.01 +) + + + + + + + + + + + +C + + + + + + + + + +D + + + + + + + + + +A + + + + + + + + + +B + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 + +-0.2 +-0.1 +0.0 +0.1 +0.2 + + + + + + +-0.2 +-0.1 +0.0 +0.1 +0.2 + + + + + +PIT +ECDF difference +ppc_pit_ecdf_grouped (correlated diff) + + diff --git a/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-grouped-correlated.svg b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-grouped-correlated.svg new file mode 100644 index 00000000..8796ffb4 --- /dev/null +++ b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-grouped-correlated.svg @@ -0,0 +1,228 @@ + + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +O +T += +0.004 + +( +α += +0.01 +) + + + + + + + + + + + + + + + +p +u +n +i +f +P +O +T += +0.002 + +( +α += +0.01 +) + + + + + + + + + + + + + +p +u +n +i +f +P +O +T += +0.34 + +( +α += +0.01 +) + + + + + + + + + + + + + +p +u +n +i +f +P +O +T += +0.34 + +( +α += +0.01 +) + + + + + + + + + + + +C + + + + + + + + + +D + + + + + + + + + +A + + + + + + + + + +B + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 + +0.00 +0.25 +0.50 +0.75 +1.00 + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 + + + + + +PIT +ECDF +ppc_pit_ecdf_grouped (correlated) + + diff --git a/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-linewidth-1.svg b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-linewidth-1.svg new file mode 100644 index 00000000..6452a003 --- /dev/null +++ b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-linewidth-1.svg @@ -0,0 +1,78 @@ + + + + + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +O +T += +0.000 + +( +α += +0.01 +) + + + +0.00 +0.25 +0.50 +0.75 +1.00 + + + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +PIT +ECDF +ppc_pit_ecdf (linewidth = 1) + + diff --git a/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-linewidth-2.svg b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-linewidth-2.svg new file mode 100644 index 00000000..f851d807 --- /dev/null +++ b/tests/testthat/_snaps/ppc-distributions/ppc-pit-ecdf-linewidth-2.svg @@ -0,0 +1,78 @@ + + + + + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +O +T += +0.000 + +( +α += +0.01 +) + + + +0.00 +0.25 +0.50 +0.75 +1.00 + + + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +PIT +ECDF +ppc_pit_ecdf (linewidth = 2) + + diff --git a/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-alpha-0-05.svg b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-alpha-0-05.svg new file mode 100644 index 00000000..82ae233c --- /dev/null +++ b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-alpha-0-05.svg @@ -0,0 +1,74 @@ + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +O +T += +0.37 + +( +α += +0.05 +) + + + +-0.10 +-0.05 +0.00 +0.05 +0.10 + + + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +LOO-PIT +ECDF difference +ppc_loo_pit_ecdf (alpha=0.05) + + diff --git a/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-changed-theme.svg b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-changed-theme.svg new file mode 100644 index 00000000..bf39eb07 --- /dev/null +++ b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-changed-theme.svg @@ -0,0 +1,70 @@ + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +O +T += +0.37 + +( +α += +0.01 +) + + + +-0.1 +0.0 +0.1 + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +LOO-PIT +ECDF difference +ppc_loo_pit_ecdf (changed theme) + + diff --git a/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-color-change.svg b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-color-change.svg new file mode 100644 index 00000000..49763437 --- /dev/null +++ b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-color-change.svg @@ -0,0 +1,78 @@ + + + + + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +O +T += +0.000 + +( +α += +0.01 +) + + + +0.00 +0.25 +0.50 +0.75 +1.00 + + + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +LOO-PIT +ECDF +ppc_loo_pit_ecdf (color change) + + diff --git a/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-correlated-piet.svg b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-correlated-piet.svg new file mode 100644 index 00000000..9eb513f8 --- /dev/null +++ b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-correlated-piet.svg @@ -0,0 +1,75 @@ + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +I +E +T += +0.19 + +( +α += +0.01 +) + + + +0.00 +0.25 +0.50 +0.75 +1.00 + + + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +LOO-PIT +ECDF +ppc_loo_pit_ecdf (correlated piet) + + diff --git a/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-correlated-pot.svg b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-correlated-pot.svg new file mode 100644 index 00000000..dedf7b1b --- /dev/null +++ b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-correlated-pot.svg @@ -0,0 +1,78 @@ + + + + + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +O +T += +0.000 + +( +α += +0.01 +) + + + +0.00 +0.25 +0.50 +0.75 +1.00 + + + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +LOO-PIT +ECDF +ppc_loo_pit_ecdf (correlated pot) + + diff --git a/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-correlated-prit.svg b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-correlated-prit.svg new file mode 100644 index 00000000..70141b36 --- /dev/null +++ b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-correlated-prit.svg @@ -0,0 +1,81 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +R +I +T += +0.001 + +( +α += +0.01 +) + + + +0.00 +0.25 +0.50 +0.75 +1.00 + + + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +LOO-PIT +ECDF +ppc_loo_pit_ecdf (correlated prit) + + diff --git a/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-default.svg b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-default.svg index 9bdd3960..83330ce3 100644 --- a/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-default.svg +++ b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-default.svg @@ -25,9 +25,9 @@ - - - + + + @@ -52,7 +52,7 @@ 0.50 0.75 1.00 -LOO PIT +LOO-PIT ECDF ppc_loo_pit_ecdf (default) diff --git a/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-diff-correlated-piet.svg b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-diff-correlated-piet.svg new file mode 100644 index 00000000..66599f42 --- /dev/null +++ b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-diff-correlated-piet.svg @@ -0,0 +1,75 @@ + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +I +E +T += +0.19 + +( +α += +0.01 +) + + + +-0.10 +-0.05 +0.00 +0.05 +0.10 + + + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +LOO-PIT +ECDF difference +ppc_loo_pit_ecdf (diff, correlated piet) + + diff --git a/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-diff-correlated-pot.svg b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-diff-correlated-pot.svg new file mode 100644 index 00000000..8de13006 --- /dev/null +++ b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-diff-correlated-pot.svg @@ -0,0 +1,78 @@ + + + + + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +O +T += +0.000 + +( +α += +0.01 +) + + + +-0.10 +-0.05 +0.00 +0.05 +0.10 + + + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +LOO-PIT +ECDF difference +ppc_loo_pit_ecdf (diff, correlated pot) + + diff --git a/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-diff-correlated-prit.svg b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-diff-correlated-prit.svg new file mode 100644 index 00000000..5a26e10d --- /dev/null +++ b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-diff-correlated-prit.svg @@ -0,0 +1,81 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +R +I +T += +0.001 + +( +α += +0.01 +) + + + +-0.10 +-0.05 +0.00 +0.05 +0.10 + + + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +LOO-PIT +ECDF difference +ppc_loo_pit_ecdf (diff, correlated prit) + + diff --git a/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-ecdf-difference.svg b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-ecdf-difference.svg index 5441468f..9e353cba 100644 --- a/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-ecdf-difference.svg +++ b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-ecdf-difference.svg @@ -25,9 +25,9 @@ - - - + + + @@ -48,7 +48,7 @@ 0.50 0.75 1.00 -LOO PIT +LOO-PIT ECDF difference ppc_loo_pit_ecdf (ecdf difference) diff --git a/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-k.svg b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-k.svg index 48f3cb24..25ab6c43 100644 --- a/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-k.svg +++ b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-k.svg @@ -25,9 +25,9 @@ - - - + + + @@ -52,7 +52,7 @@ 0.50 0.75 1.00 -LOO PIT +LOO-PIT ECDF ppc_loo_pit_ecdf (K) diff --git a/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-linewidth-1.svg b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-linewidth-1.svg new file mode 100644 index 00000000..3fa8cd92 --- /dev/null +++ b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-linewidth-1.svg @@ -0,0 +1,78 @@ + + + + + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +O +T += +0.000 + +( +α += +0.01 +) + + + +0.00 +0.25 +0.50 +0.75 +1.00 + + + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +LOO-PIT +ECDF +ppc_loo_pit_ecdf (linewidth = 1) + + diff --git a/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-linewidth-2.svg b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-linewidth-2.svg new file mode 100644 index 00000000..ff6e2f74 --- /dev/null +++ b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-linewidth-2.svg @@ -0,0 +1,78 @@ + + + + + + + + + + + + + + + + + + + + + + + + + +p +u +n +i +f +P +O +T += +0.000 + +( +α += +0.01 +) + + + +0.00 +0.25 +0.50 +0.75 +1.00 + + + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +LOO-PIT +ECDF +ppc_loo_pit_ecdf (linewidth = 2) + + diff --git a/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-no-help-text.svg b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-no-help-text.svg new file mode 100644 index 00000000..310a4987 --- /dev/null +++ b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-no-help-text.svg @@ -0,0 +1,58 @@ + + + + + + + + + + + + + + + + + + + + + + + + +-0.10 +-0.05 +0.00 +0.05 +0.10 + + + + + + + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 +LOO-PIT +ECDF difference +ppc_loo_pit_ecdf (no help_text) + + diff --git a/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-prob.svg b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-prob.svg index dadcb9e6..97b89103 100644 --- a/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-prob.svg +++ b/tests/testthat/_snaps/ppc-loo/ppc-loo-pit-ecdf-prob.svg @@ -25,9 +25,9 @@ - - - + + + @@ -52,7 +52,7 @@ 0.50 0.75 1.00 -LOO PIT +LOO-PIT ECDF ppc_loo_pit_ecdf (prob) diff --git a/tests/testthat/test-helpers-ppc.R b/tests/testthat/test-helpers-ppc.R index 36f921d3..043159a7 100644 --- a/tests/testthat/test-helpers-ppc.R +++ b/tests/testthat/test-helpers-ppc.R @@ -124,3 +124,11 @@ test_that("ecdf_intervals returns right dimensions and values", { expect_equal(min(lims$lower), 0) expect_equal(max(lims$lower), 100) }) + +# display p-values in plots ------------------------------------------------ +test_that("formatting of p-values works as expected", { + expect_equal(fmt_p(0.446), "0.45") + expect_equal(fmt_p(0.045), "0.045") + expect_equal(fmt_p(0.0045), "0.005") + expect_equal(fmt_p(0.00045), "0.000") +}) \ No newline at end of file diff --git a/tests/testthat/test-ppc-distributions.R b/tests/testthat/test-ppc-distributions.R index 34e1c82f..586016c2 100644 --- a/tests/testthat/test-ppc-distributions.R +++ b/tests/testthat/test-ppc-distributions.R @@ -104,13 +104,112 @@ test_that("ppc_dots returns a ggplot object", { }) test_that("ppc_pit_ecdf, ppc_pit_ecdf_grouped returns a ggplot object", { + # Independent method (default) expect_gg(ppc_pit_ecdf(y, yrep, interpolate_adj = FALSE)) + expect_gg(ppc_pit_ecdf(y, yrep, method = "independent", interpolate_adj = FALSE)) expect_gg(ppc_pit_ecdf_grouped(y, yrep, group = group, interpolate_adj = FALSE)) - expect_message(ppc_pit_ecdf(pit = runif(100)), "'pit' specified") + + # Correlated method + expect_gg(ppc_pit_ecdf(y, yrep, method = "correlated")) + expect_gg(ppc_pit_ecdf(y, yrep, method = "correlated", plot_diff = TRUE)) + expect_gg(ppc_pit_ecdf(y, yrep, method = "correlated", test = "PRIT")) + expect_gg(ppc_pit_ecdf(y, yrep, method = "correlated", test = "PIET")) + expect_gg(ppc_pit_ecdf_grouped(y, yrep, group = group, method = "correlated")) + expect_gg(ppc_pit_ecdf_grouped(y, yrep, group = group, method = "correlated", plot_diff = TRUE)) + expect_gg(ppc_pit_ecdf_grouped(y, yrep, group = group, method = "correlated", test = "PRIT")) + expect_gg(ppc_pit_ecdf_grouped(y, yrep, group = group, method = "correlated", test = "PIET")) + + # Specify 'pit' directly (with y/yrep still supplied) expect_message( - ppc_pit_ecdf_grouped(pit = runif(length(group)), group = group, interpolate_adj = FALSE), + ppc_pit_ecdf_grouped( + y = y, yrep = yrep, pit = runif(length(group)), + group = group, interpolate_adj = FALSE + ), "'pit' specified" ) + + # No y/yrep provided with pit -> no ignored-input message but "independent" method message + expect_message( + ppc_pit_ecdf_grouped(pit = runif(length(group)), group = group, + method = "independent", interpolate_adj = FALSE), + "The 'independent' method is superseded by the 'correlated' method." + ) +}) + +test_that("ppc_pit_ecdf method validation and ignored-argument warnings", { + # Invalid method + expect_error(ppc_pit_ecdf(y, yrep, method = "bogus")) + + # method = "correlated" warns about interpolate_adj + expect_message( + ppc_pit_ecdf(y, yrep, method = "correlated", interpolate_adj = TRUE), + "ignoring.*interpolate_adj" + ) + + # method = "independent" warns about test and gamma + expect_message( + ppc_pit_ecdf(y, yrep, method = "independent", test = "POT", + interpolate_adj = FALSE), + "ignoring.*test" + ) + expect_message( + ppc_pit_ecdf(y, yrep, method = "independent", test = "POT", gamma = 0.5, + interpolate_adj = FALSE), + "ignoring.*test, gamma" + ) + + # Invalid test type for correlated + expect_error( + ppc_pit_ecdf(y, yrep, method = "correlated", test = "INVALID") + ) +}) + +test_that("ppc_pit_ecdf correlated method validates gamma", { + expect_error( + ppc_pit_ecdf(y, yrep, method = "correlated", gamma = -1), + regexp = "gamma must be in" + ) +}) + +test_that("ppc_pit_ecdf_grouped method validation and ignored-argument warnings", { + # Invalid method + expect_error(ppc_pit_ecdf_grouped(y, yrep, group = group, method = "bogus")) + + # method = "correlated" warns about interpolate_adj + expect_message( + ppc_pit_ecdf_grouped( + y, yrep, group = group, method = "correlated", interpolate_adj = TRUE + ), + "ignoring.*interpolate_adj" + ) + + # method = "independent" warns about correlated-only args + expect_message( + ppc_pit_ecdf_grouped( + y, yrep, group = group, method = "independent", + test = "POT", interpolate_adj = FALSE + ), + "ignoring.*test" + ) + expect_message( + ppc_pit_ecdf_grouped( + y, yrep, group = group, method = "independent", + test = "POT", gamma = 0.5, interpolate_adj = FALSE + ), + "ignoring.*test, gamma" + ) + + # Invalid test type for correlated + expect_error( + ppc_pit_ecdf_grouped(y, yrep, group = group, method = "correlated", test = "INVALID") + ) +}) + +test_that("ppc_pit_ecdf_grouped correlated method validates gamma", { + expect_error( + ppc_pit_ecdf_grouped(y, yrep, group = group, method = "correlated", gamma = -1), + regexp = "gamma must be in" + ) }) test_that("ppc_freqpoly_grouped returns a ggplot object", { @@ -503,6 +602,7 @@ test_that("ppc_pit_ecdf, ppc_pit_ecdf_grouped renders correctly", { testthat::skip_if_not_installed("vdiffr") skip_on_r_oldrel() + # Independent method p_base <- ppc_pit_ecdf(y, yrep, interpolate_adj = FALSE) g_base <- ppc_pit_ecdf_grouped(y, yrep, group = group, interpolate_adj = FALSE) p_diff <- ppc_pit_ecdf(y, yrep, plot_diff = TRUE, interpolate_adj = FALSE) @@ -512,4 +612,304 @@ test_that("ppc_pit_ecdf, ppc_pit_ecdf_grouped renders correctly", { vdiffr::expect_doppelganger("ppc_pit_ecdf_grouped (default)", g_base) vdiffr::expect_doppelganger("ppc_pit_ecdf (diff)", p_diff) vdiffr::expect_doppelganger("ppc_pit_ecdf_grouped (diff)", g_diff) + + # Correlated method + p_corr <- ppc_pit_ecdf(y, yrep, method = "correlated") + vdiffr::expect_doppelganger("ppc_pit_ecdf (correlated)", p_corr) + + p_corr_diff <- ppc_pit_ecdf(y, yrep, method = "correlated", plot_diff = TRUE) + vdiffr::expect_doppelganger("ppc_pit_ecdf (correlated diff)", p_corr_diff) + + p_corr_prit <- ppc_pit_ecdf(y, yrep, method = "correlated", test = "PRIT") + vdiffr::expect_doppelganger("ppc_pit_ecdf (correlated PRIT)", p_corr_prit) + + p_corr_piet <- ppc_pit_ecdf(y, yrep, method = "correlated", test = "PIET") + vdiffr::expect_doppelganger("ppc_pit_ecdf (correlated PIET)", p_corr_piet) + + g_corr <- ppc_pit_ecdf_grouped(y, yrep, group = group, method = "correlated") + vdiffr::expect_doppelganger("ppc_pit_ecdf_grouped (correlated)", g_corr) + + g_corr_diff <- ppc_pit_ecdf_grouped( + y, yrep, group = group, method = "correlated", plot_diff = TRUE + ) + vdiffr::expect_doppelganger("ppc_pit_ecdf_grouped (correlated diff)", g_corr_diff) }) + +test_that("ppc_pit_ecdf with method correlated renders different tests correctly", { + set.seed(2025) + pit <- 1 - (1 - runif(300))^(1.2) + + p_cor_pot <- ppc_pit_ecdf( + pit = pit, + method = "correlated" + ) + vdiffr::expect_doppelganger("ppc_pit_ecdf (correlated pot)", p_cor_pot) + + p_cor_prit <- ppc_pit_ecdf( + pit = pit, + method = "correlated", + test = "PRIT" + ) + vdiffr::expect_doppelganger("ppc_pit_ecdf (correlated prit 2)", p_cor_prit) + + p_cor_piet <- ppc_pit_ecdf( + pit = pit, + method = "correlated", + test = "PIET" + ) + vdiffr::expect_doppelganger("ppc_pit_ecdf (correlated piet 2)", p_cor_piet) +}) + +test_that("ppc_pit_ecdf with plot_diff=TRUE and method correlated renders different tests correctly", { + set.seed(2025) + pit <- 1 - (1 - runif(300))^(1.2) + + p_cor_pot <- ppc_pit_ecdf( + pit = pit, + method = "correlated", + plot_diff = TRUE + ) + vdiffr::expect_doppelganger("ppc_pit_ecdf (diff, correlated pot)", p_cor_pot) + + p_cor_prit <- ppc_pit_ecdf( + pit = pit, + method = "correlated", + test = "PRIT", + plot_diff = TRUE + ) + vdiffr::expect_doppelganger("ppc_pit_ecdf (diff, correlated prit)", p_cor_prit) + + p_cor_piet <- ppc_pit_ecdf( + pit = pit, + method = "correlated", + test = "PIET", + plot_diff = TRUE + ) + vdiffr::expect_doppelganger("ppc_pit_ecdf (diff, correlated piet)", p_cor_piet) +}) + +test_that("ppc_pit_ecdf renders different linewidths and colors correctly", { + set.seed(2025) + pit <- 1 - (1 - runif(300))^(1.2) + + p_cor_lw1 <- ppc_pit_ecdf( + pit = pit, + method = "correlated", + linewidth = 1. + ) + vdiffr::expect_doppelganger("ppc_pit_ecdf (linewidth = 1)", p_cor_lw1) + + p_cor_lw2 <- ppc_pit_ecdf( + pit = pit, + method = "correlated", + linewidth = 2. + ) + vdiffr::expect_doppelganger("ppc_pit_ecdf (linewidth = 2)", p_cor_lw2) + + p_cor_col <- ppc_pit_ecdf( + pit = pit, + method = "correlated", + color = c(ecdf = "darkblue", highlight = "red") + ) + vdiffr::expect_doppelganger("ppc_pit_ecdf (color change)", p_cor_col) +}) + + +# Test PIT computation branches ------------------------------------------------ +# use monkey-patching to test whether the correct branch of the +# PIT computation is taken + +testthat::test_that("ppc_pit_ecdf takes correct PIT computation branch", { + skip_on_cran() + skip_if_not_installed("loo") + skip_on_r_oldrel() + skip_if(packageVersion("rstantools") <= "2.4.0") + + ppc_pit_ecdf_patched <- ppc_pit_ecdf + + body(ppc_pit_ecdf_patched)[[ + # Replace the PIT computation block (the large if/else if/else) + # with a version that emits diagnostics + which(sapply(as.list(body(ppc_pit_ecdf)), function(e) { + if (!is.call(e)) return(FALSE) + identical(e[[1]], as.name("if")) && + grepl("pareto_pit", paste(deparse(e[[2]]), collapse = " ")) + }))[1] + ]] <- quote({ + + if (isTRUE(pareto_pit) && is.null(pit)) { + message("[PIT BRANCH] Pareto-smoothed LOO PIT") + suggested_package("rstantools") + y <- validate_y(y) + yrep <- validate_predictions(yrep, length(y)) + + pit <- posterior::pareto_pit(x = yrep, y = y, weights = NULL, log = TRUE) + K <- K %||% length(pit) + + } else if (!is.null(pit)) { + message("[PIT BRANCH] Pre-supplied PIT") + pit <- validate_pit(pit) + K <- K %||% length(pit) + + ignored <- c( + if (!missing(y) && !is.null(y)) "y", + if (!missing(yrep) && !is.null(yrep)) "yrep" + ) + if (length(ignored) > 0) { + inform(paste0("As 'pit' specified; ignoring: ", + paste(ignored, collapse = ", "), ".")) + } + + } else { + message("[PIT BRANCH] Empirical PIT") + pit <- ppc_data(y, yrep) %>% + group_by(.data$y_id) %>% + dplyr::group_map( + ~ mean(.x$value[.x$is_y] > .x$value[!.x$is_y]) + + runif(1, max = mean(.x$value[.x$is_y] == .x$value[!.x$is_y])) + ) %>% + unlist() + K <- K %||% min(nrow(yrep) + 1, 1000) + } + }) + + # | yrep | y | pit | method | test | pareto_pit | approach | + # |------|---|-----|-------------|------|------------|--------------------| + # | x | x | | independent | NULL | FALSE | empirical pit | + # | | | x | independent | NULL | FALSE | | + # | x | x | | independent | NULL | TRUE | compute pareto-pit | + # | x | x | | correlated | POT | TRUE | compute pareto-pit | + # | | | x | correlated | POT | FALSE | | + # | x | x | | correlated | PIET | TRUE | compute pareto-pit | + # | | | x | correlated | PIET | FALSE | | + # | x | x | | correlated | PRIT | FALSE | empirical pit | + # | | | x | correlated | PRIT | FALSE | | + + pits <- rstantools::loo_pit(vdiff_loo_yrep, vdiff_loo_y, vdiff_loo_lw) + + # method = independent ------------------------------------------ + expect_message( + ppc_pit_ecdf_patched( + vdiff_loo_y, + vdiff_loo_yrep, + method = "independent" + ), + regexp = "\\[PIT BRANCH\\] Empirical PIT" + ) + + expect_message( + ppc_pit_ecdf_patched( + vdiff_loo_y, + vdiff_loo_yrep, + method = "independent", + pareto_pit = TRUE + ), + regexp = "\\[PIT BRANCH\\] Pareto-smoothed LOO PIT" + ) + + expect_message( + ppc_pit_ecdf_patched( + method = "independent", + pit = pits, + ), + regexp = "\\[PIT BRANCH\\] Pre-supplied PIT" + ) + + # method = correlated + POT test ------------------------------- + expect_message( + ppc_pit_ecdf_patched( + vdiff_loo_y, + vdiff_loo_yrep, + method = "correlated" + ), + regexp = "\\[PIT BRANCH\\] Pareto-smoothed LOO PIT" + ) + + expect_message( + ppc_pit_ecdf_patched( + vdiff_loo_y, + vdiff_loo_yrep, + method = "correlated", + pareto_pit = FALSE + ), + regexp = "\\[PIT BRANCH\\] Empirical PIT" + ) + + expect_message( + ppc_pit_ecdf_patched( + method = "correlated", + pit = pits, + ), + regexp = "\\[PIT BRANCH\\] Pre-supplied PIT" + ) + + # method = correlated + PIET test ------------------------------- + expect_message( + ppc_pit_ecdf_patched( + vdiff_loo_y, + vdiff_loo_yrep, + method = "correlated", + test = "PIET" + ), + regexp = "\\[PIT BRANCH\\] Pareto-smoothed LOO PIT" + ) + + expect_message( + ppc_pit_ecdf_patched( + method = "correlated", + test = "PIET", + pit = pits, + ), + regexp = "\\[PIT BRANCH\\] Pre-supplied PIT" + ) + + # method = correlated + PRIT test ------------------------------- + expect_message( + ppc_pit_ecdf_patched( + vdiff_loo_y, + vdiff_loo_yrep, + method = "correlated", + test = "PRIT" + ), + regexp = "\\[PIT BRANCH\\] Empirical PIT" + ) + + expect_message( + ppc_pit_ecdf_patched( + method = "correlated", + test = "PRIT", + pit = pits, + ), + regexp = "\\[PIT BRANCH\\] Pre-supplied PIT" + ) +}) + +test_that("ppc_pit_ecdf works with pareto_pit method", { + skip_on_cran() + skip_if_not_installed("brms") + skip_if_not_installed("rstanarm") + + data("roaches", package = "rstanarm") + roaches$sqrt_roach1 <- sqrt(roaches$roach1) + + fit_p <- brms::brm(y ~ sqrt_roach1 + treatment + senior + offset(log(exposure2)), + data = roaches, + family = poisson, + prior = brms::prior(normal(0, 1), class = b), + refresh = 0) + + fit_p <- brms::add_criterion(fit_p, criterion = "loo") + fit_p <- brms::add_criterion(fit_p, criterion = "loo", moment_match = TRUE, overwrite = TRUE) + fit_nb <- update(fit_p, family = brms::negbinomial) + + expect_gg(brms::pp_check(fit_nb, type = "pit_ecdf")) + + draws <- brms::posterior_predict(fit_nb) + y <- roaches$y + + expect_gg(ppc_pit_ecdf( + y = y, yrep = draws, method = "correlated" + )) + + expect_gg(brms::pp_check(fit_nb, type = "pit_ecdf", method = "correlated")) +}) \ No newline at end of file diff --git a/tests/testthat/test-ppc-loo.R b/tests/testthat/test-ppc-loo.R index a722ad46..b7163c2f 100644 --- a/tests/testthat/test-ppc-loo.R +++ b/tests/testthat/test-ppc-loo.R @@ -104,7 +104,7 @@ test_that("ppc_loo_pit_ecdf returns a ggplot object", { } else { ll1 <- p1$labels } - expect_equal(ll1$x, "LOO PIT") + expect_equal(ll1$x, "LOO-PIT") expect_equal(ll1$y, "ECDF") expect_equal(p1$data, p2$data) expect_gg(p3 <- ppc_loo_pit_ecdf(y, yrep, lw, plot_diff = TRUE)) @@ -116,6 +116,98 @@ test_that("ppc_loo_pit_ecdf returns a ggplot object", { expect_equal(ll3$y, "ECDF difference") }) +test_that("ppc_loo_pit_ecdf with method='correlated' validates input correctly", { + set.seed(2025) + pit <- 1 - (1 - runif(300))^(1.2) + y_mock <- 1:length(pit) + + expect_message( + ppc_loo_pit_ecdf(pit = pit, method = "correlated", interpolate_adj = FALSE), + "As method = 'correlated' specified; ignoring: interpolate_adj." + ) + expect_message( + ppc_loo_pit_ecdf(pit = pit, method = "independent", y = y_mock), + "As 'pit' specified; ignoring: y." + ) + expect_message( + ppc_loo_pit_ecdf(pit = pit, method = "independent", gamma = 1.0), + "As method = 'independent' specified; ignoring: gamma." + ) + expect_message( + ppc_loo_pit_ecdf(pit = pit, method = "independent", test = "POT"), + "As method = 'independent' specified; ignoring: test." + ) +}) + +test_that("ppc_loo_pit_ecdf with method='correlated' returns ggplot object", { + skip_if_not_installed("rstanarm") + skip_if_not_installed("loo") + + # Test with POT-C (default) + expect_gg(p1 <- ppc_loo_pit_ecdf(y, yrep, lw, method = "correlated")) + + # Test with PRIT-C + expect_gg(p2 <- ppc_loo_pit_ecdf(y, yrep, lw, method = "correlated", test = "PRIT")) + + # Test with PIET-C + expect_gg(p3 <- ppc_loo_pit_ecdf(y, yrep, lw, method = "correlated", test = "PIET")) + + # Test with plot_diff = TRUE + expect_gg(p4 <- ppc_loo_pit_ecdf(y, yrep, lw, method = "correlated", plot_diff = TRUE)) + + # Test with gamma specified + expect_gg(p5 <- ppc_loo_pit_ecdf(y, yrep, lw, method = "correlated", gamma = 0.1)) +}) + +test_that("ppc_loo_pit_ecdf method argument works correctly", { + skip_if_not_installed("rstanarm") + skip_if_not_installed("loo") + + # Test default (should inform about upcoming change) + expect_message( + p1 <- ppc_loo_pit_ecdf(y, yrep, lw), + "In the next major release" + ) + expect_gg(p1) + + # Test explicit independent method (should inform about supersession) + expect_message( + p2 <- ppc_loo_pit_ecdf(y, yrep, lw, method = "independent"), + "superseded by the 'correlated' method" + ) + expect_gg(p2) + + # Test correlated method (no message expected) + expect_gg(p3 <- ppc_loo_pit_ecdf(y, yrep, lw, method = "correlated")) + + # Test that independent and correlated produce different plots + expect_true(!identical(p2$data, p3$data) || !identical(p2$layers, p3$layers)) +}) + +test_that("ppc_loo_pit_ecdf correlated method handles edge cases", { + skip_if_not_installed("rstanarm") + skip_if_not_installed("loo") + + set.seed(2026) + + # Test with small sample + small_pit <- runif(10) + expect_gg(p1 <- ppc_loo_pit_ecdf(pit = small_pit, method = "correlated")) + + # Test with perfect uniform + uniform_pit <- seq(0, 1, length.out = 100) + expect_gg(p2 <- ppc_loo_pit_ecdf(pit = uniform_pit, method = "correlated")) + + # Test with extreme values + extreme_pit <- c(rep(0, 10), rep(1, 10), runif(80)) + expect_gg(p3 <- ppc_loo_pit_ecdf(pit = extreme_pit, method = "correlated")) + + # Test with single value (edge case) + single_pit <- 0.5 + expect_error(ppc_loo_pit_ecdf(pit = single_pit, method = "correlated")) + expect_gg(p5 <- ppc_loo_pit_ecdf(pit = single_pit, method = "correlated", test = "PIET")) +}) + test_that("ppc_loo_pit functions work when pit specified instead of y, yrep, and lw", { skip_if_not_installed("rstanarm") skip_if_not_installed("loo") @@ -134,7 +226,7 @@ test_that("ppc_loo_pit functions work when pit specified instead of y, yrep, and expect_gg(ppc_loo_pit_ecdf(pit = rep(pits, 4))) expect_message( p1 <- ppc_loo_pit_ecdf(y = y, yrep = yrep, lw = lw, pit = rep(pits, 4)), - "'pit' specified so ignoring 'y','yrep','lw' if specified" + "As 'pit' specified; ignoring: y, yrep, lw." ) expect_message( p2 <- ppc_loo_pit_ecdf(pit = rep(pits, 4)) @@ -149,7 +241,6 @@ test_that("ppc_loo_pit functions work when pit specified instead of y, yrep, and ) }) - test_that("ppc_loo_intervals returns ggplot object", { skip_if_not_installed("rstanarm") skip_if_not_installed("loo") @@ -212,6 +303,40 @@ test_that("error if subset is bigger than num obs", { ) }) +test_that("ppc_loo_pit_ecdf works with pareto_pit method", { + skip_on_cran() + skip_if_not_installed("brms") + skip_if_not_installed("rstanarm") + + data("roaches", package = "rstanarm") + roaches$sqrt_roach1 <- sqrt(roaches$roach1) + + fit_zinb <- + brms::brm(brms::bf(y ~ sqrt_roach1 + treatment + senior + offset(log(exposure2)), + zi ~ sqrt_roach1 + treatment + senior + offset(log(exposure2))), + family = brms::zero_inflated_negbinomial(), data = roaches, + prior = c(brms::prior(normal(0, 1), class = "b"), + brms::prior(normal(0, 1), class = "b", dpar = "zi"), + brms::prior(normal(0, 1), class = "Intercept", dpar = "zi")), + seed = 1704009, refresh = 1000) + + fit_zinb <- brms::add_criterion(fit_zinb, criterion = "loo", save_psis = TRUE) + fit_zinb <- brms::add_criterion(fit_zinb, criterion = "loo", save_psis = TRUE, + moment_match = TRUE, overwrite = TRUE) + + draws <- brms::posterior_predict(fit_zinb) + psis_object <- brms::loo(fit_zinb, save_psis = TRUE)$psis_object + y <- roaches$y + + expect_gg(ppc_loo_pit_ecdf( + y = y, yrep = draws, psis_object = psis_object, method = "correlated" + )) + + expect_gg(brms::pp_check( + fit_zinb, type = "loo_pit_ecdf", moment_match = TRUE, method = "correlated" + )) +}) + # Visual tests ------------------------------------------------------------ @@ -312,6 +437,85 @@ test_that("ppc_loo_ribbon renders correctly", { vdiffr::expect_doppelganger("ppc_loo_ribbon (subset)", p_custom) }) +test_that("ppc_loo_pit_ecdf with method correlated renders different tests correctly", { + set.seed(2025) + pit <- 1 - (1 - runif(300))^(1.2) + + p_cor_pot <- ppc_loo_pit_ecdf( + pit = pit, + method = "correlated" + ) + vdiffr::expect_doppelganger("ppc_loo_pit_ecdf (correlated pot)", p_cor_pot) + + p_cor_prit <- ppc_loo_pit_ecdf( + pit = pit, + method = "correlated", + test = "PRIT" + ) + vdiffr::expect_doppelganger("ppc_loo_pit_ecdf (correlated prit)", p_cor_prit) + + p_cor_piet <- ppc_loo_pit_ecdf( + pit = pit, + method = "correlated", + test = "PIET" + ) + vdiffr::expect_doppelganger("ppc_loo_pit_ecdf (correlated piet)", p_cor_piet) +}) + +test_that("ppc_loo_pit_ecdf with plot_diff=TRUE and method correlated renders different tests correctly", { + set.seed(2025) + pit <- 1 - (1 - runif(300))^(1.2) + + p_cor_pot <- ppc_loo_pit_ecdf( + pit = pit, + method = "correlated", + plot_diff = TRUE + ) + vdiffr::expect_doppelganger("ppc_loo_pit_ecdf (diff, correlated pot)", p_cor_pot) + + p_cor_prit <- ppc_loo_pit_ecdf( + pit = pit, + method = "correlated", + test = "PRIT", + plot_diff = TRUE + ) + vdiffr::expect_doppelganger("ppc_loo_pit_ecdf (diff, correlated prit)", p_cor_prit) + + p_cor_piet <- ppc_loo_pit_ecdf( + pit = pit, + method = "correlated", + test = "PIET", + plot_diff = TRUE + ) + vdiffr::expect_doppelganger("ppc_loo_pit_ecdf (diff, correlated piet)", p_cor_piet) +}) + +test_that("ppc_loo_pit_ecdf renders different linewidths and colors correctly", { + set.seed(2025) + pit <- 1 - (1 - runif(300))^(1.2) + + p_cor_lw1 <- ppc_loo_pit_ecdf( + pit = pit, + method = "correlated", + linewidth = 1. + ) + vdiffr::expect_doppelganger("ppc_loo_pit_ecdf (linewidth = 1)", p_cor_lw1) + + p_cor_lw2 <- ppc_loo_pit_ecdf( + pit = pit, + method = "correlated", + linewidth = 2. + ) + vdiffr::expect_doppelganger("ppc_loo_pit_ecdf (linewidth = 2)", p_cor_lw2) + + p_cor_col <- ppc_loo_pit_ecdf( + pit = pit, + method = "correlated", + color = c(ecdf = "darkblue", highlight = "red") + ) + vdiffr::expect_doppelganger("ppc_loo_pit_ecdf (color change)", p_cor_col) +}) + test_that("ppc_loo_pit_ecdf renders correctly", { skip_on_cran() skip_if_not_installed("vdiffr") @@ -351,6 +555,262 @@ test_that("ppc_loo_pit_ecdf renders correctly", { K = 100 ) vdiffr::expect_doppelganger("ppc_loo_pit_ecdf (ecdf difference)", p_custom) + + p_custom <- ppc_loo_pit_ecdf( + vdiff_loo_y, + vdiff_loo_yrep, + psis_object = psis_object, + method = "correlated", + plot_diff = TRUE, + prob = 0.95 + ) + vdiffr::expect_doppelganger("ppc_loo_pit_ecdf (alpha=0.05)", p_custom) + + p_custom <- ppc_loo_pit_ecdf( + vdiff_loo_y, + vdiff_loo_yrep, + psis_object = psis_object, + method = "correlated", + plot_diff = TRUE, + prob = 0.95, + help_text = FALSE + ) + vdiffr::expect_doppelganger("ppc_loo_pit_ecdf (no help_text)", p_custom) + + + theme_set(bayesplot::theme_default(base_family = "sans", base_size = 12)) + p_custom <- ppc_loo_pit_ecdf( + vdiff_loo_y, + vdiff_loo_yrep, + psis_object = psis_object, + method = "correlated", + plot_diff = TRUE + ) + vdiffr::expect_doppelganger("ppc_loo_pit_ecdf (changed theme)", p_custom) +}) + +# Test PIT computation branches ------------------------------------------------ +# use monkey-patching to test whether the correct branch of the +# PIT computation is taken + +testthat::test_that("ppc_loo_pit_ecdf takes correct PIT computation branch", { + skip_on_cran() + skip_if_not_installed("loo") + skip_on_r_oldrel() + skip_if(packageVersion("rstantools") <= "2.4.0") + + ppc_loo_pit_ecdf_patched <- ppc_loo_pit_ecdf + + body(ppc_loo_pit_ecdf_patched)[[ + # Replace the PIT computation block (the large if/else if/else) + # with a version that emits diagnostics + which(sapply(as.list(body(ppc_loo_pit_ecdf)), function(e) { + if (!is.call(e)) return(FALSE) + identical(e[[1]], as.name("if")) && + grepl("pareto_pit", paste(deparse(e[[2]]), collapse = " ")) + }))[1] + ]] <- quote({ + + if (isTRUE(pareto_pit) && is.null(pit)) { + message("[PIT BRANCH] Pareto-smoothed LOO PIT") + suggested_package("rstantools") + y <- validate_y(y) + yrep <- validate_predictions(yrep, length(y)) + lw <- .get_lw(lw, psis_object) + stopifnot(identical(dim(yrep), dim(lw))) + pit <- posterior::pareto_pit(x = yrep, y = y, weights = lw, log = TRUE) + K <- K %||% length(pit) + + } else if (!is.null(pit)) { + message("[PIT BRANCH] Pre-supplied PIT") + pit <- validate_pit(pit) + K <- K %||% length(pit) + + ignored <- c( + if (!missing(y) && !is.null(y)) "y", + if (!missing(yrep) && !is.null(yrep)) "yrep", + if (!is.null(lw)) "lw" + ) + if (length(ignored) > 0) { + inform(paste0("As 'pit' specified; ignoring: ", + paste(ignored, collapse = ", "), ".")) + } + + } else { + message("[PIT BRANCH] Standard LOO PIT") + suggested_package("rstantools") + y <- validate_y(y) + yrep <- validate_predictions(yrep, length(y)) + lw <- .get_lw(lw, psis_object) + stopifnot(identical(dim(yrep), dim(lw))) + pit <- pmin(1, rstantools::loo_pit(object = yrep, y = y, lw = lw)) + K <- K %||% min(nrow(yrep) + 1, 1000) + } + }) + + # | yrep | y | lw | psis_object | pit | method | test | pareto_pit | approach | + # |------|---|----|-------------|-----|-------------|------|------------|--------------------| + # | x | x | x | | | independent | NULL | FALSE (D) | compute loo-pit | + # | x | x | | x | | independent | NULL | FALSE (D) | compute loo-pit | + # | x | x | x | | | independent | NULL | TRUE | compute pareto-pit | + # | x | x | | x | | independent | NULL | TRUE | compute pareto-pit | + # | | | | | x | independent | NULL | FALSE | | + # | x | x | x | | | correlated | POT | TRUE | compute pareto-pit | + # | x | x | | x | | correlated | POT | TRUE | compute pareto-pit | + # | | | | | x | correlated | POT | FALSE | | + # | x | x | x | | | correlated | PIET | TRUE | compute pareto-pit | + # | x | x | | x | | correlated | PIET | TRUE | compute pareto-pit | + # | | | | | x | correlated | PIET | FALSE | | + # | x | x | x | | | correlated | PRIT | FALSE | compute loo-pit | + # | x | x | | x | | correlated | PRIT | FALSE | compute loo-pit | + # | | | | | x | correlated | PRIT | FALSE | | + + psis_object <- suppressWarnings(loo::psis(-vdiff_loo_lw)) + pits <- rstantools::loo_pit(vdiff_loo_yrep, vdiff_loo_y, vdiff_loo_lw) + + # method = independent ------------------------------------------ + expect_message( + ppc_loo_pit_ecdf_patched( + vdiff_loo_y, + vdiff_loo_yrep, + method = "independent", + lw = vdiff_loo_lw + ), + regexp = "\\[PIT BRANCH\\] Standard LOO PIT" + ) + + expect_message( + ppc_loo_pit_ecdf_patched( + vdiff_loo_y, + vdiff_loo_yrep, + method = "independent", + psis_object = psis_object, + ), + regexp = "\\[PIT BRANCH\\] Standard LOO PIT" + ) + + expect_message( + ppc_loo_pit_ecdf_patched( + vdiff_loo_y, + vdiff_loo_yrep, + method = "independent", + psis_object = psis_object, + pareto_pit = TRUE + ), + regexp = "\\[PIT BRANCH\\] Pareto-smoothed LOO PIT" + ) + + expect_message( + ppc_loo_pit_ecdf_patched( + method = "independent", + pit = pits, + ), + regexp = "\\[PIT BRANCH\\] Pre-supplied PIT" + ) + + # method = correlated + POT test ------------------------------- + expect_message( + ppc_loo_pit_ecdf_patched( + vdiff_loo_y, + vdiff_loo_yrep, + method = "correlated", + lw = vdiff_loo_lw + ), + regexp = "\\[PIT BRANCH\\] Pareto-smoothed LOO PIT" + ) + + expect_message( + ppc_loo_pit_ecdf_patched( + vdiff_loo_y, + vdiff_loo_yrep, + method = "correlated", + lw = vdiff_loo_lw, + pareto_pit = FALSE + ), + regexp = "\\[PIT BRANCH\\] Standard LOO PIT" + ) + + expect_message( + ppc_loo_pit_ecdf_patched( + vdiff_loo_y, + vdiff_loo_yrep, + method = "correlated", + psis_object = psis_object, + ), + regexp = "\\[PIT BRANCH\\] Pareto-smoothed LOO PIT" + ) + + expect_message( + ppc_loo_pit_ecdf_patched( + method = "correlated", + pit = pits, + ), + regexp = "\\[PIT BRANCH\\] Pre-supplied PIT" + ) + + # method = correlated + PIET test ------------------------------- + expect_message( + ppc_loo_pit_ecdf_patched( + vdiff_loo_y, + vdiff_loo_yrep, + method = "correlated", + test = "PIET", + lw = vdiff_loo_lw + ), + regexp = "\\[PIT BRANCH\\] Pareto-smoothed LOO PIT" + ) + + expect_message( + ppc_loo_pit_ecdf_patched( + vdiff_loo_y, + vdiff_loo_yrep, + method = "correlated", + test = "PIET", + psis_object = psis_object, + ), + regexp = "\\[PIT BRANCH\\] Pareto-smoothed LOO PIT" + ) + + expect_message( + ppc_loo_pit_ecdf_patched( + method = "correlated", + test = "PIET", + pit = pits, + ), + regexp = "\\[PIT BRANCH\\] Pre-supplied PIT" + ) + + # method = correlated + PRIT test ------------------------------- + expect_message( + ppc_loo_pit_ecdf_patched( + vdiff_loo_y, + vdiff_loo_yrep, + method = "correlated", + test = "PRIT", + lw = vdiff_loo_lw + ), + regexp = "\\[PIT BRANCH\\] Standard LOO PIT" + ) + + expect_message( + ppc_loo_pit_ecdf_patched( + vdiff_loo_y, + vdiff_loo_yrep, + method = "correlated", + test = "PRIT", + psis_object = psis_object, + ), + regexp = "\\[PIT BRANCH\\] Standard LOO PIT" + ) + + expect_message( + ppc_loo_pit_ecdf_patched( + method = "correlated", + test = "PRIT", + pit = pits, + ), + regexp = "\\[PIT BRANCH\\] Pre-supplied PIT" + ) }) diff --git a/vignettes/loo-pit-correlated-tests.Rmd b/vignettes/loo-pit-correlated-tests.Rmd new file mode 100644 index 00000000..a1898b68 --- /dev/null +++ b/vignettes/loo-pit-correlated-tests.Rmd @@ -0,0 +1,448 @@ +--- +title: "Model checking using ppc_pit_ecdf() and ppc_loo_pit_ecdf()" +author: "bayesplot team" +date: "`r Sys.Date()`" +output: + rmarkdown::html_vignette: + toc: true + toc_depth: 3 +params: + EVAL: !r identical(Sys.getenv("NOT_CRAN"), "true") +vignette: > + %\VignetteIndexEntry{Model checking using ppc_pit_ecdf() and ppc_loo_pit_ecdf()} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, child="children/SETTINGS-knitr.txt"} +``` + +```{r pkgs, include=FALSE} +library(ggplot2) +library(patchwork) +library(loo) +library(brms) + +bayesplot_theme_set() + +set.seed(2026) +``` + +```{r, eval=FALSE} +library(bayesplot) +library(ggplot2) +library(patchwork) +library(loo) +library(brms) + +set.seed(2026) +``` + +## Introduction +This vignette focuses on pointwise predictive checks based on the probability integral transform (PIT). The PIT is defined as the model’s univariate predictive cumulative distribution function (CDF) evaluated given the observation. For a well-calibrated model, these transformations result in values that should be approximately uniformly distributed between 0 and 1 (Tesso & Vehtari, [2026](#Tesso2026); Säilynoja et al., [2022](#Säilynoja2022)). + +In `bayesplot`, PIT-based checks are implemented in the `ppc_pit_ecdf()` and `ppc_loo_pit_ecdf()` functions, which plot the empirical CDF of the PIT values against the expected uniform distribution. The `ppc_pit_ecdf_grouped()` function extends this functionality to allow for grouping by covariates. + +**New Correlation-Aware Method** + +Following the work of Tesso & Vehtari ([2026](#Tesso2026)), `bayesplot` now offers a new method for evaluating the uniformity of PIT values that accounts for their correlation structure. This method is implemented via the argument `method = "correlated"` and is recommended for most applications. While the previous approach is retained via `method = "independent"` for backward compatibility and is the current default. + +This vignette focuses specifically on the changes introduced by this new correlation-aware method. For background information on graphical uniformity tests using PIT, see Säilynoja et al. ([2022](#Säilynoja2022)). For a more general discussion on the use of Leave-One-Out Cross-Validation (LOO-CV), see Vehtari et al. ([2017](#Vehtari2017), [2024](#Vehtari2024)), among others. + +## The `method` argument +### `method = "independent"` (superseded) +When `method = "independent"` is selected, simultaneous confidence bands for the ECDF are constructed under the assumption that the PIT values are both independent and uniform (Säilynoja et al., [2022](#Säilynoja2022)). However, if this independence assumption is violated, the resulting bands can be too wide, which reduces the test's sensitivity to actual miscalibration (Tesso & Vehtari, [2026](#Tesso2026)). + +**Deprecation and Compatibility** + +As of `bayesplot vX.X.X`, the `"independent"` method is officially superseded. To maintain backward compatibility, `"independent"` remains the current default; however, using it will now trigger a message informing the user: +``` +"The 'independent' method is superseded by the 'correlated' method." +``` +This is intended to encourage a transition to the `"correlated"` method, which will become the default in a future release. + +### `method = "correlated"` (new, recommended) +This method employes one of three dependence-aware uniformity tests (selected via the `test` argument) to compute a global p-value for the null hypothesis of uniformity. Unlike the independent method, it accounts for the correlation among PIT values (Tesso & Vehtari, [2026](#Tesso2026)). + +Instead of drawing traditional confidence bands, the plot highlights ECDF regions in red where the pointwise contribution to the test statistic is largest. This visualization makes it easier to diagnose the *type* and *location* of miscalibration. + +## Reading the plots for different (mis)calibration scenarios +The shape of the ECDF curve provides direct insight into *how* a predictive distribution is miscalibrated. To illustrate this, the following examples utilize simulated scenarios where "observed" values (`y`) are drawn from a `normal(0, sd)` distribution, while "replicated" values (`yrep`) are generated from a non-central t-distribution. By varying the degrees of freedom (`df`) and non-centrality parameter (`ncp`), we can simulate and visualize several distinct types of miscalibration. + +```{r, echo=FALSE} +plot_scenarios <- function( + df = NULL, ncp = NULL, sd = 1, n = 300, S = 100, seed = 2026, + xlim = c(-5, 5) + ) { + set.seed(seed) + y <- rnorm(n, mean = 0, sd = sd) + yrep <- matrix(rt(n * S, df = df, ncp = ncp), nrow = S, ncol = n) + + yrep_df <- data.frame( + value = as.vector(t(yrep)), + draw = rep(seq_len(nrow(yrep)), each = ncol(yrep)) + ) + + p1 <- ggplot(yrep_df, aes(x = .data$value, group = .data$draw, color = "yrep")) + + geom_density(linewidth = 0.5, alpha = 0.3) + + geom_density( + data = data.frame(value = y), + aes(x = .data$value, color = "y"), + linewidth = 1, + inherit.aes = FALSE + ) + + scale_color_manual( + values = c(yrep = "lightblue", y = "darkblue"), + breaks = c("yrep", "y"), + name = NULL + ) + + xlim(c(xlim[1], xlim[2])) + + yaxis_ticks(FALSE) + + xlab("y / yrep") + + bayesplot_theme_get() + + theme( + legend.position = "top", + legend.direction = "horizontal" + ) + p2 <- ppc_pit_ecdf(y = y, yrep = yrep, method = "correlated", prob = 0.95) + + labs(subtitle = "method = 'correlated'") + p3 <- ppc_pit_ecdf(y = y, yrep = yrep, method = "independent", prob = 0.95) + + labs(subtitle = "method = 'independent'") + + bayesplot_grid( + p1, p2, p3, + grid_args = list(ncol = 3, widths = c(1.1, 1.1, 1.1)) + ) +} + +``` + +### Calibrated model +We begin by simulating a well-calibrated scenario where the predictive distribution (`ypred`) closely matches the data-generating process (`y`). In this case, the PIT values should be approximately uniform, and the ECDF should closely follow the identity (diagonal) line. + +$$ +\begin{aligned} +y &\sim \text{Normal}(0, 1) \\ +y_{rep} &\sim \text{Student}_{\nu = 100}(0, 1) \quad (\approx \text{Normal}(0, 1)) +\end{aligned} +$$ +```{r calibrated, fig.width=7, fig.height=3, fig.asp=NULL, out.width="80%", fig.align="center", warning=FALSE, message=FALSE} +plot_scenarios(df = 100, ncp = 0) +``` + +### Predictive distribution too wide (overdispersed) +When the predictive distribution is too wide, the model overestimates the probability mass in the tails. As a result, observed values tend to fall near the *centre* of the predictive distribution more often than expected. This leads to an `S`-shaped ECDF curve that typically falls below the diagonal line in the lower tail and rises above it in the upper tail. + +$$ +\begin{aligned} +y &\sim \text{Normal}(0, 1) \\ +y_{rep} &\sim \text{Student}_{\nu = 1}(0, 1) \quad (= \text{Cauchy}(0, 1)) +\end{aligned} +$$ + +```{r too-wide, fig.width=7, fig.height=3, fig.asp=NULL, out.width="80%", fig.align="center", warning=FALSE, message=FALSE} +plot_scenarios(df = 1, ncp = 0, xlim = c(-7, 7)) +``` + +### Predictive distribution too narrow (underdispersed) +When the predictive distribution is too narrow, the model underestimates the potential for extreme values, causing observations to fall frequently in the tails. This results in a reverse `S`-shaped ECDF curve that typically lies above the diagonal line in the lower tail and below it in the upper tail. + +$$ +\begin{aligned} +y &\sim \text{Normal}(0, 0.5) \\ +y_{rep} &\sim \text{Student}_{\nu = 100}(0, 1) \quad (\approx \text{Normal}(0, 1)) +\end{aligned} +$$ + +```{r too-narrow, fig.width=7, fig.height=3, fig.asp=NULL, out.width="80%", fig.align="center", warning=FALSE, message=FALSE} +plot_scenarios(df = 100, ncp = 0, sd = 2) +``` + +### Systematically biased predictions (location shift) +When the model's predictions are systematically shifted relative to the data, observations tend to fall consistently in one tail of the predictive distribution. This results in an ECDF curve that "bows" aways from the identity line: + ++ Underestimation (Positive Bias): If the model predicts values that are systematically too low, the PIT values will cluster near 1. The ECDF curve will stay below the diagonal line. ++ Overestimation (Negative Bias): If the model predicts values that are systematically too high, the PIT values will cluster near 0. The ECDF curve will stay above the diagonal line. + +$$ +\begin{aligned} +y &\sim \text{Normal}(1, 1) \\ +y_{rep} &\sim \text{Student}_{\nu = 100}(0, 1) \quad (\approx \text{Normal}(0, 1)) +\end{aligned} +$$ + +```{r biased, fig.width=7, fig.height=3, fig.asp=NULL, out.width="80%", fig.align="center", warning=FALSE, message=FALSE} +plot_scenarios(df = 100, ncp = 1) +``` + +## The three uniformity-tests within `method = "correlated"` +When using `method = "correlated"`, you select the testing procedure via the `test` argument. The three options correspond to the three procedures proposed by Tesso & Vehtari ([2026](#Tesso2026)): + +### `test = "POT"` +Pointwise Order Tests Combination (`POT`) is based on order statistics of the PIT values and is the **recommended default** for most workflows. It should ideally use continuous LOO-PITs and works for PSIS-LOO weighted rank-based PITs. + +### `test = "PRIT"` +Pointwise Rank-based Individual Tests Combination (`PRIT`) is designed for **rank-based or discrete PIT values**. It is most compatible with PITs computed as normalized ranks and is weaker than "POT" in testing weighted rank based LOO-PIT values for the same reason previously outlined. + +### `test = "PIET"` +Pointwise Inverse-CDF Evaluation Tests Combination (`PIET`) applies an inverse-CDF transformation before testing and is **tail-sensitive**. It is exclusively recommended for detecting deviations in tails (or outliers) and requires continuous PIT values for reliable behavior. + +## Showcase functions for PIT-based checks using nabiximols data +To illustrate the functions for PIT-based checks, we use data from the study "Nabiximols for the Treatment of Cannabis Dependence: A Randomized Clinical Trial" by Lintzeris et al. ([2019](#Lintzeris2019)). The dataset includes 128 participants (`id`) in two groups (`group`): Placebo and Nabiximols. + +Participants received a 12-week treatment involving weekly clinical reviews, structured counseling, and flexible medication doses, dispensed weekly, of up to 32 sprays daily (86.4 mg tetrahydrocannabinol and 80 mg cannabidiol). The number of cannabis use days (`cu`) during the previous 28 days (`set`) was recorded at 0, 4, 8, and 12 weeks (`week`). + +We fit a simple linear mixed model using the `brms` package to demonstrate the use of `ppc_pit_ecdf()`, `ppc_loo_pit_ecdf()`, and `ppc_pit_ecdf_grouped()`. + +For a more detailed analysis of this data within the Bayesian workflow, see the [Nabiximols Case Study](https://users.aalto.fi/~ave/casestudies/Nabiximols/nabiximols.html#ref-Lintzeris+etal:2019:Nabiximols) by Aki Vehtari. + +```{r naboximols-data, warning=FALSE, message=FALSE} +id <- factor(c(1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4, 5, 5, 6, 6, 6, 6, 7, 7, 8, 8, 8, 8, 9, 9, 9, 10, 10, 10, 10, 11, 11, 11, 12, 12, 13, 14, 15, 16, 16, 17, 18, 18, 18, 18, 19, 20, 20, 20, 20, 21, 21, 21, 21, 22, 22, 23, 23, 23, 24, 24, 24, 24, 25, 25, 25, 25, 26, 27, 27, 28, 28, 28, 28, 29, 30, 30, 30, 30, 31, 31, 32, 32, 32, 32, 33, 33, 33, 34, 34, 34, 35, 35, 36, 36, 37, 37, 37, 37, 38, 39, 39, 39, 39, 40, 40, 40, 41, 42, 42, 42, 42, 43, 43, 43, 43, 44, 44, 45, 45, 46, 46, 46, 46, 47, 47, 47, 47, 48, 48, 49, 49, 49, 50, 50, 50, 50, 51, 51, 51, 52, 52, 52, 52, 53, 53, 53, 53, 54, 54, 55, 55, 55, 55, 56, 57, 57, 57, 57, 58, 58, 58, 58, 59, 59, 59, 59, 60, 60, 60, 60, 61, 61, 61, 62, 63, 63, 64, 64, 64, 65, 65, 65, 65, 66, 66, 66, 66, 67, 67, 67, 67, 68, 68, 68, 69, 69, 69, 69, 70, 70, 70, 70, 71, 71, 71, 71, 72, 73, 73, 73, 73, 74, 74, 74, 75, 76, 76, 76, 76, 77, 77, 77, 77, 78, 78, 78, 79, 79, 79, 79, 80, 80, 80, 80, 81, 81, 81, 81, 82, 82, 83, 83, 84, 84, 84, 85, 85, 85, 86, 86, 86, 86, 87, 87, 87, 87, 88, 88, 88, 88, 89, 89, 89, 89, 90, 90, 90, 90, 91, 91, 91, 91, 92, 92, 92, 92, 93, 93, 93, 93, 94, 94, 94, 94, 95, 95, 95, 95, 96, 96, 96, 96, 97, 97, 97, 98, 98, 98, 98, 99, 99, 99, 99, 100, 101, 101, 101, 102, 102, 102, 102, 103, 103, 103, 103, 104, 104, 105, 105, 105, 105, 106, 106, 106, 106, 107, 107, 107, 107, 108, 108, 108, 108, 109, 109, 109, 109, 110, 110, 111, 111, 112, 112, 112, 112, 113, 113, 113, 113, 114, 115, 115, 115, 115, 116, 116, 116, 116, 117, 117, 117, 117, 118, 118, 119, 119, 119, 119, 120, 120, 120, 120, 121, 121, 121, 122, 123, 123, 123, 123, 124, 124, 124, 125, 125, 125, 125, 126, 126, 126, 126, 127, 127, 128)) +group <- factor(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0), levels = 0:1, labels = c("placebo", "nabiximols")) +week <- factor(c(0, 4, 8, 12, 0, 4, 8, 0, 4, 8, 0, 4, 8, 12, 0, 4, 0, 4, 8, 12, 0, 4, 0, 4, 8, 12, 0, 4, 8, 0, 4, 8, 12, 0, 4, 8, 0, 4, 0, 0, 0, 0, 4, 0, 0, 4, 8, 12, 0, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 0, 4, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 0, 4, 0, 4, 8, 12, 0, 0, 4, 8, 12, 0, 4, 0, 4, 8, 12, 0, 4, 8, 0, 4, 12, 0, 8, 0, 4, 0, 4, 8, 12, 0, 0, 4, 8, 12, 0, 4, 8, 0, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 0, 4, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 0, 4, 12, 0, 4, 8, 12, 0, 4, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 0, 4, 8, 12, 0, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 12, 0, 0, 4, 0, 4, 8, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 0, 4, 8, 12, 0, 4, 8, 0, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 0, 4, 0, 4, 8, 0, 4, 8, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 0, 4, 8, 12, 0, 4, 8, 12, 0, 0, 4, 8, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 0, 4, 0, 4, 8, 12, 0, 4, 8, 12, 0, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 12, 0, 0, 4, 8, 12, 0, 4, 12, 0, 4, 8, 12, 0, 4, 8, 12, 0, 4, 0)) +cu <- c(13, 12, 12, 12, 28, 0, NA, 16, 9, 2, 28, 28, 28, 28, 28, NA, 28, 28, 17, 28, 28, NA, 16, 0, 0, NA, 28, 28, 28, 28, 17, 0, NA, 28, 27, 28, 28, 26, 24, 28, 28, 28, 25, 28, 26, 28, 18, 16, 28, 28, 7, 0, 2, 28, 2, 4, 1, 28, 28, 16, 28, 28, 24, 26, 15, 28, 25, 17, 1, 8, 28, 24, 27, 28, 28, 28, 28, 28, 27, 28, 28, 28, 28, 20, 28, 28, 28, 28, 12, 28, NA, 17, 15, 14, 28, 0, 28, 28, 28, 0, 0, 0, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 21, 24, 28, 27, 28, 28, 26, NA, 28, NA, 20, 2, 3, 7, 28, 1, 19, 8, 21, 7, 28, 28, 20, 28, 28, 28, 24, 20, 17, 11, 25, 25, 28, 26, 28, 24, 17, 16, 27, 14, 28, 28, 28, 28, 28, 28, 14, 13, 4, 24, 28, 28, 28, 21, 28, 21, 26, 28, 28, 0, 0, 28, 23, 20, 28, 20, 16, 28, 28, 28, 10, 1, 1, 2, 28, 28, 28, 28, 18, 22, 9, 15, 28, 9, 1, 20, 18, 20, 24, 28, 28, 28, 19, 28, 28, 28, 28, 28, 28, 28, 28, 28, 4, 14, 20, 28, 28, 0, 0, 0, 28, 20, 9, 24, 28, 28, 28, 28, 28, 21, 28, 28, 14, 24, 28, 23, 0, 0, 0, 28, NA, 28, NA, 28, 15, NA, 12, 25, NA, 28, 2, 0, 0, 28, 10, 0, 0, 28, 0, 0, 0, 23, 0, 0, 0, 28, 0, 0, 0, 28, 0, 0, 0, 28, 2, 1, 0, 21, 14, 7, 8, 28, 28, 28, 0, 28, 28, 20, 18, 24, 0, 0, 0, 28, 15, NA, 28, 1, 1, 2, 28, 1, 0, 0, 28, 28, 14, 21, 25, 19, 16, 13, 28, 28, 28, 28, 28, 28, 28, 27, 19, 21, 18, 1, 0, 0, 28, 28, 28, 28, 28, 24, 27, 28, 18, 0, 3, 8, 28, 28, 28, 9, 20, 25, 20, 12, 19, 0, 0, 0, 27, 28, 0, 0, 0, 20, 17, 16, 14, 28, 7, 0, 1, 28, 24, 28, 25, 23, 20, 28, 14, 16, 7, 28, 28, 26, 28, 28, 26, 28, 28, 28, 24, 20, 28, 28, 28, 28, 28, 8, 6, 4, 28, 20, 28) +set <- rep(28, length(cu)) +cu_df <- data.frame(id, group, week, cu, set) + +# remove the rows with NA's as they are not used for posteriors anyway +cu_df <- cu_df |> + tidyr::drop_na(cu) + +fit_normal <- brms::brm(formula = cu ~ group*week + (1 | id), + data = cu_df, + family = gaussian(), + prior = c(brms::prior(normal(14, 1.5), class = Intercept), + brms::prior(normal(0, 11), class = b), + brms::prior(cauchy(1,2), class = sd)), + save_pars = brms::save_pars(all=TRUE), + seed = 1234, + refresh = 0) +fit_normal <- brms::add_criterion(fit_normal, criterion="loo", save_psis=TRUE) +``` + +### `ppc_pit_ecdf()` + +```{r ppc-pit-example-1, fig.width=7, fig.height=3, fig.asp=NULL, out.width="80%", warning=FALSE, message=FALSE} +p1 <- ppc_pit_ecdf( + y = cu_df$cu, + yrep = brms::posterior_predict(fit_normal), + method = "independent" +) + +p2 <- ppc_pit_ecdf( + y = cu_df$cu, + yrep = brms::posterior_predict(fit_normal), + method = "correlated" +) + +p1 | p2 +``` + +```{r ppc-pit-example-2, fig.width=7, fig.height=3, fig.asp=NULL, out.width="80%", warning=FALSE, message=FALSE} +p1 <- ppc_pit_ecdf( + y = cu_df$cu, + yrep = brms::posterior_predict(fit_normal), + method = "independent", + plot_diff = TRUE +) + +p2 <- ppc_pit_ecdf( + y = cu_df$cu, + yrep = brms::posterior_predict(fit_normal), + method = "correlated", + plot_diff = TRUE +) + +p1 | p2 +``` + +### `ppc_loo_pit_ecdf()` + +```{r loo-pit-example-1, fig.width=7, fig.height=3, fig.asp=NULL, out.width="80%", warning=FALSE, message=FALSE} +p1 <- ppc_loo_pit_ecdf( + y = cu_df$cu, + yrep = brms::posterior_predict(fit_normal), + method = "independent", + lw = weights(loo(fit_normal)$psis_object) +) + +p2 <- ppc_loo_pit_ecdf( + y = cu_df$cu, + yrep = brms::posterior_predict(fit_normal), + method = "correlated", + lw = weights(loo(fit_normal)$psis_object) +) + +p1 | p2 +``` + +```{r loo-pit-example-2, fig.width=7, fig.height=3, fig.asp=NULL, out.width="80%", warning=FALSE, message=FALSE} +p1 <- ppc_loo_pit_ecdf( + y = cu_df$cu, + yrep = brms::posterior_predict(fit_normal), + method = "independent", + lw = weights(loo(fit_normal)$psis_object), + plot_diff = TRUE +) + +p2 <- ppc_loo_pit_ecdf( + y = cu_df$cu, + yrep = brms::posterior_predict(fit_normal), + method = "correlated", + lw = weights(loo(fit_normal)$psis_object), + plot_diff = TRUE +) + +p1 | p2 +``` + +### `ppc_pit_ecdf_grouped()` + +```{r grouped-example, fig.width=7, fig.height=3, fig.asp=NULL, out.width="80%", warning=FALSE, message=FALSE} +ppc_pit_ecdf_grouped( + y = cu_df$cu, + yrep = brms::posterior_predict(fit_normal), + lw = weights(loo(fit_normal)$psis_object), + group = cu_df$group, + method = "correlated" + ) +``` + +## Using `brms::pp_check()` +It is also possible to use `brms::pp_check()` with `type = "loo_pit_ecdf"` to perform the same testing and plotting procedure as `ppc_loo_pit_ecdf()`. The following example demonstrates this using the same fitted model as above with `method = "correlated"`. +```{r pp-check-example-1, fig.show="hold", out.width="80%", warning=FALSE, message=FALSE} +brms::pp_check( + fit_normal, + type = "loo_pit_ecdf", + method = "correlated" +) +``` + +Or for `method = "independent"`: + +```{r pp-check-example-2, fig.show="hold", out.width="80%", warning=FALSE, message=FALSE} +brms::pp_check(fit_normal, method = "independent", type = "pit_ecdf") +``` + +## Additional arguments +With the introduction of the `method = "correlated"` option, the three functions now have additional arguments that control the appearance and behavior of the plot when using correlated testing procedures. These arguments are: + +### `gamma` +When `method = "correlated"` and the global test rejects uniformity, the plot uses color coding to highlight influential regions. The `gamma` argument is a non-negative numeric threshold: a point is highlighted if its pointwise contribution to the test statistic exceeds `gamma`. The default is `0`, which highlights all points that contribute at all when the test rejects. Increasing `gamma` raises the bar, flagging only the most strongly deviant points. + +```{r gamma-example, fig.width=7, fig.height=3, fig.asp=NULL, out.width="80%", fig.align="center"} +set.seed(2026) +pit <- rbeta(300, 1, 1.2) + +# Highlight only the most strongly influential regions +p1 <- ppc_loo_pit_ecdf( + pit = pit, + method = "correlated", + plot_diff = TRUE, + gamma = 0 + ) + + labs(subtitle = "gamma = 0") + +# Highlight more broadly +p2 <- ppc_loo_pit_ecdf( + pit = pit, method = "correlated", + plot_diff = TRUE, + gamma = 80 + ) + + labs(subtitle = "gamma = 80") + +p1 | p2 +``` + +### `linewidth` and `color` +`linewidth` controls the width of the ECDF line (default `0.3`). `color` is a named character vector with two elements controlling the appearance of the plot: `c(ecdf = "grey60", highlight = "red")`. Whereby `ecdf` is the color of the main ECDF line; `highlight` is the color used for the +flagged regions when the global test rejects. Both can be overridden: + +```{r color-example, fig.width=7, fig.height=3, fig.asp=NULL, out.width="80%", fig.align="center"} +pit <- rbeta(300, 2, 2) + +p1 <- ppc_loo_pit_ecdf( + pit = pit, + method = "correlated" +) + + labs(subtitle = "default style") + +p2 <- ppc_loo_pit_ecdf( + pit = pit, + method = "correlated", + color = c(ecdf = "steelblue", highlight = "orange"), + linewidth = 0.5 +) + + labs(subtitle = "custom color / linewidth") + +p1 | p2 +``` + +### `help_text` +Setting `help_text = TRUE` adds an annotation to the plot showing the `test` name, (1-`prob`) level ($\alpha$), and global p-value. This is the **default** (`help_text = TRUE`). Set `help_text = FALSE` to suppress the annotation, for example when embedding the plot in a figure with a caption. + +```{r help-text-example, fig.width=7, fig.height=3, fig.asp=NULL, out.width="80%", fig.align="center"} +pit <- rbeta(300, 2, 2) + +p1 <- ppc_loo_pit_ecdf( + pit = pit, + method = "correlated", + plot_diff = TRUE +) + + labs(subtitle = "help_text = TRUE") + +p2 <- ppc_loo_pit_ecdf( + pit = pit, + method = "correlated", + plot_diff = TRUE, + help_text = FALSE +) + + labs(subtitle = "help_text = FALSE") + +p1 | p2 +``` + +### `pareto_pit` +The `pareto_pit` argument controls how PIT values are computed internally when precomputed `pit` values are not supplied. + +**Default behavior (recommended):** leave `pareto_pit` at its default value. The function chooses the appropriate PIT computation automatically based on `method` and `test`: + +| method | test | `pareto_pit` default | +|:---|:------:|:------| +| `"independent"` | -- | `FALSE` (empirical/rank-based PIT) | +| `"correlated"` | `"PRIT"` | `FALSE` (empirical/rank-based PIT) | +| `"correlated"` | `"POT"` | `TRUE` (Pareto-smoothed PIT) | +| `"correlated"` | `"PIET"` | `TRUE` (Pareto-smoothed PIT) | + +When `pareto_pit = TRUE`, `ppc_loo_pit_ecdf()` requires LOO importance weights to be supplied (via `lw` or `psis_object`), since the LOO-PIT computation is weighted. `ppc_pit_ecdf()` and `ppc_pit_ecdf_grouped()` do not require weights; they call `posterior::pareto_pit()` with `weights = NULL`. + +**When to override:** You can set `pareto_pit` manually if you need to force a specific PIT computation path for diagnostic comparison. In practice this is rarely necessary. + +## Summary of new arguments + +| Argument | Functions | Values | Default | +|---|---|---|---| +| `method` | all three | `"independent"`, `"correlated"` | `NULL` → falls back to `"independent"` with a deprecation notice | +| `test` | all three | `"POT"`, `"PRIT"`, `"PIET"` | `"POT"` | +| `gamma` | all three | numeric $\ge 0$ | `0` | +| `linewidth` | all three | numeric | `0.3` | +| `color` | all three | named character vector | `c(ecdf = "grey60", highlight = "red")` | +| `help_text` | all three | `TRUE`, `FALSE` | `TRUE` | +| `pareto_pit` | all three | `TRUE`, `FALSE` | `TRUE` if `test` is `"POT"` or `"PIET"` and no `pit` supplied; `FALSE` otherwise | + +## References + +Tesso, H., & Vehtari, A. (2026). LOO-PIT predictive model checking. http://arxiv.org/abs/2603.02928 + + +Säilynoja, T., Bürkner, P.-C., and Vehtari, A. (2022). Graphical test for discrete uniformity and its applications in goodness-of-fit evaluation and multiple sample comparison. *Statistics and Computing*, 32, 32. https://link.springer.com/10.1007/s11222-022-10090-6 + + +Vehtari, A., Gelman, A., and Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. *Statistics and Computing*, 27(5), 1413–1432. https://link.springer.com/article/10.1007/S11222-016-9696-4 + + +Vehtari, A., Simpson, D., Gelman, A., Yao, Y., and Gabry, J. (2024). Pareto smoothed importance sampling. *Journal of Machine Learning Research*, 25(72), 1–58. https://www.jmlr.org/papers/v25/19-556.html + + +Lintzeris, Nicholas, Anjali Bhardwaj, Llewellyn Mills, Adrian Dunlop, Jan Copeland, Iain McGregor, Raimondo Bruno, et al. 2019. “Nabiximols for the Treatment of Cannabis Dependence: A Randomized Clinical Trial.” JAMA Internal Medicine 179 (9): 1242–53. https://jamanetwork.com/journals/jamainternalmedicine/fullarticle/2737918 + \ No newline at end of file