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..d790e60c --- /dev/null +++ b/vignettes/loo-pit-correlated-tests.Rmd @@ -0,0 +1,598 @@ +--- +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) +``` + +## Setup +```{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. + +## 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) + p3 <- ppc_pit_ecdf(y = y, yrep = yrep, method = "correlated", + plot_diff = TRUE, prob = 0.95) + + 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}(0, 1) \\ +y_{rep} &\sim \text{Student}_{\nu = 100}(-1, 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 `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. + +```{r comparison, fig.width=7, fig.height=3, fig.asp=NULL, out.width="80%", fig.align="center"} +set.seed(2026) +pit <- rbeta(300, 1, 1.2) + +p1 <- ppc_loo_pit_ecdf( + pit = pit, + method = "independent" + ) + + labs(subtitle = "method = 'independent'") + +p2 <- ppc_loo_pit_ecdf( + pit = pit, method = "correlated" + ) + + labs(subtitle = "method = 'correlated'") + +p1 | p2 +``` + +## 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, 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 = "correlated" +) + +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 ppc-loo-pit-example, 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 = "correlated", + 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), + 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, fig.show="hold", out.width="80%", warning=FALSE, message=FALSE} +brms::pp_check( + fit_normal, + type = "loo_pit_ecdf", + method = "correlated" +) +``` + +## 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. It follows the standard PIT +computation procedure, except that in the tails it substitutes the empirical CDF +(ECDF) with the CDF of a fitted generalized Pareto distribution. This reduces +variability and yields non-zero, non-unity PIT values even outside the range of +the reference draws — making it particularly useful for stabilizing PIT +uniformity checks, where a single 0 or 1 can otherwise introduce large +variation. + +The following figure compares simulated PIT values from `rstantools::loo_pit()` +(x-axis) with those from `posterior::pareto_pit()` (y-axis). Both axes are +logit-transformed, so PIT values of exactly 0 or 1 map to `-Inf` or `Inf`, +respectively; these are shown as red bars on the left (for 0) and right (for 1) +edge of the y-axis. The plot illustrates how `pareto_pit` replaces these extreme + values with finite, less extreme alternatives, thereby reducing variability. + +```{r pareto_pit-example, out.width="80%", fig.align="center", warning=FALSE, echo=FALSE} +set.seed(122) +n <- 300 +S <- 100 + +# True data-generating process +y <- rnorm(n, mean = 0, sd = 1) +yrep <- matrix(rt(n * S, df = 100), nrow = S, ncol = n) + +# Log-likelihood: depends only on y (not the draw), so compute once and broadcast +log_lik_vec <- dt(y, df = 100, log = TRUE) +log_lik <- matrix(rep(log_lik_vec, each = S), nrow = S, ncol = n) + +# LOO log importance weights via PSIS +psis_result <- loo::psis(-log_lik) +lw <- weights(psis_result, log = TRUE) + +# Compute and sort PITs once; reuse sorted vectors for both tail subsets +pit <- rstantools::loo_pit(object = yrep, y = y, lw = lw) +pit2 <- posterior::pareto_pit(x = yrep, y = y, weights = lw, log = TRUE) + +pit_sorted <- sort(pit) +pit2_sorted <- sort(pit2) + +# Lower- and upper-tail subsets (M points each) +df <- data.frame(pit = pit_sorted, pit2 = pit2_sorted) +df$clipped_0 <- df$pit == 0 | df$pit2 == 0 +df$clipped_1 <- df$pit == 1 | df$pit2 == 1 + +df_normal <- df[c(!df$clipped_0, !df$clipped_1), ] +df_clipped_0 <- df[df$clipped_0, ] +df_clipped_1 <- df[df$clipped_1, ] + +breaks <- c(0.001, 0.01, 0.1, 0.5, 0.9, 0.99, 0.999) + +ggplot(mapping = aes(x = pit, y = pit2)) + + geom_point(data = df_normal, colour = "black", size = 1) + + geom_rug(data = df_clipped_0, colour = "red", linewidth = 0.7, + length = unit(0.03, "npc"), sides = "l") + + geom_rug(data = df_clipped_1, colour = "red", linewidth = 0.7, + length = unit(0.03, "npc"), sides = "r") + + geom_abline(slope = 1, intercept = 0, linetype = "dashed") + + scale_x_continuous(transform = "logit", breaks = breaks, labels = breaks) + + scale_y_continuous(transform = "logit", breaks = breaks, labels = breaks) + + labs(x = "pit() [logit transform]", y = "pareto_pit() [logit transform]", + title = "PIT values from pit() vs. pareto_pit()") + + bayesplot_theme_get() +``` + +**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