diff --git a/.gitignore b/.gitignore index 695e6e77..74301abb 100644 --- a/.gitignore +++ b/.gitignore @@ -16,3 +16,23 @@ release-prep.R # vscode/positron/etc settings .vscode/* + +# ignore draft vignette +vignettes/ppc_loo_pit.Rmd +vignettes/test_greyscale.qmd + +# internal knitting of vignettes +*_cache +*_files +vignettes/*.quarto +vignettes/*.html +vignettes/.gitignore + +# visual regression test +tests/vdiffr.Rout.fail + +# agent skills (Posit/Cursor skills installation) +.agent/ +.agents/ +skills/ +skills-lock.json \ No newline at end of file diff --git a/DESCRIPTION b/DESCRIPTION index 9f47b08b..26dd4dbf 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -43,6 +43,7 @@ Imports: tidyselect, utils Suggests: + brms (>= 2.23.0), ggdist, ggfortify, gridExtra (>= 2.2.1), diff --git a/R/helpers-ppc.R b/R/helpers-ppc.R index e268b315..cd263d09 100644 --- a/R/helpers-ppc.R +++ b/R/helpers-ppc.R @@ -595,3 +595,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 d86a8b4c..3920ca60 100644 --- a/R/ppc-distributions.R +++ b/R/ppc-distributions.R @@ -648,6 +648,39 @@ 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()` 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()`, 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()` when `method = "correlated"`, which +#' dependence-aware test to use: `"POT"`, `"PRIT"`, or `"PIET"`. +#' Defaults to `"POT"`. +#' @param gamma For `ppc_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 linewidth For `ppc_pit_ecdf()` when `method = "correlated"`, the line width of the ECDF +#' and highlighting points. Defaults to 0.3. +#' @param color For `ppc_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_pit_ecdf()` when `method = "correlated"`, a boolean +#' defining whether to add informative text to the plot. Defaults to `TRUE`. +#' @param pareto_pit For `ppc_pit_ecdf()`, 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. +#' When `TRUE`, this calls `posterior::pareto_pit()`, which is proposed in +#' \url{https://github.com/stan-dev/posterior/pull/435}. This code path requires that PR to be merged +#' and released in the posterior package. #' @rdname PPC-distributions #' ppc_pit_ecdf <- function(y, @@ -657,59 +690,310 @@ ppc_pit_ecdf <- function(y, K = NULL, prob = .99, plot_diff = FALSE, - interpolate_adj = NULL) { + interpolate_adj = NULL, + method = "independent", + test = NULL, + gamma = NULL, + linewidth = NULL, + color = NULL, + help_text = NULL, + pareto_pit = NULL + ) { + # Expected input combinations for ppc_pit_ecdf() based on method and test choices: + # | 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 | | + 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)) { + # --------------------------------------------------------------------------- + # Internal helpers + # --------------------------------------------------------------------------- + + .warn_ignored <- function(method_name, args) { + inform(paste0( + "As method = ", method_name, " specified; ignoring: ", + paste(args, collapse = ", "), "." + )) + } + + # --------------------------------------------------------------------------- + # Resolve and validate `method` + # --------------------------------------------------------------------------- + + method <- match.arg(method, choices = c("independent", "correlated")) + + + # --------------------------------------------------------------------------- + # Method-specific defaults, validation, and pareto_pit flag + # --------------------------------------------------------------------------- + + 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 applies for correlated only when the test is POT or PIET + # (not PRIT, which uses empirical PIT). + 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 applies for independent whenever pareto_pit = TRUE. + pareto_pit <- pareto_pit %||% FALSE + } + ) + + # --------------------------------------------------------------------------- + # Compute PIT values + # --------------------------------------------------------------------------- + + if (isTRUE(pareto_pit) && is.null(pit)) { + # --- Pareto-smoothed PIT --- + suggested_package("rstantools") + y <- validate_y(y) + yrep <- validate_predictions(yrep, length(y)) + + # TODO: posterior::pareto_pit() from https://github.com/stan-dev/posterior/pull/435 - requires that PR merged + 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) + + # Warn about any ignored arguments. + 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 --- 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) + } + + # --------------------------------------------------------------------------- + # Shared ECDF setup + # --------------------------------------------------------------------------- + + 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" + + # =========================================================================== + # Correlated method + # =========================================================================== + + if (method == "correlated") { + + # Compute the per-observation test statistics (sorted for Shapley values) + # and the combined Cauchy p-value. + # TODO: posterior::uniformity_test() from https://github.com/stan-dev/posterior/pull/435 - requires that PR merged + test_res <- posterior::uniformity_test(pit = pit, test = test) + p_value_CCT <- test_res$pvalue + pointwise_contrib <- test_res$pointwise + + # Validate gamma against computed Shapley values. + 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) + + # Build the main ECDF data frame over a dense grid that includes pit values + # so step discontinuities are rendered exactly. + 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 + ) + + # Sorted pit data frame (used for highlighting suspicious points). + pit_sorted <- sort(pit) + df_pit <- tibble::tibble( + pit = pit_sorted, + ecdf_val = ecdf_pit_fn(pit_sorted) - plot_diff * pit_sorted + ) + + # --- Base plot ----------------------------------------------------------- + + 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) + + # --- Highlight suspicious regions ---------------------------------------- + + 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, ] + + # Highlight contiguous groups as coloured step segments. + 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 + ) + } + + # Highlight isolated suspicious points as dots. + 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 + ) + } + } + } + + # --- Annotation and axis scaling ----------------------------------------- + + 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 - ) - 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" + + # =========================================================================== + # Independent method + # =========================================================================== + + gamma_indep <- adjust_gamma(N = n_obs, K = K, prob = prob, + interpolate_adj = interpolate_adj) + + lims <- ecdf_intervals(gamma = gamma_indep, N = n_obs, K = K) + + # `lims` contains K + 1 elements (including the boundary at 0); drop it so + # lengths match `unit_interval`. + 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 diff --git a/R/ppc-loo.R b/R/ppc-loo.R index 88a6692c..37427b5f 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 #' @@ -394,12 +396,39 @@ 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 not normally be modified by the user, +#' except for development purposes. When `TRUE`, this calls `posterior::pareto_pit()`, which is +#' proposed in \url{https://github.com/stan-dev/posterior/pull/435}. This code path requires +#' that PR to be merged and released in the posterior package. +#' @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, @@ -409,63 +438,326 @@ 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) { + # Expected input combinations for ppc_loo_pit_ecdf() based on method and test choices: + # | yrep | y | lw | psis_object | pit | method | test | pareto_pit | approach | + # |------|---|----|-------------|-----|-------------|------|------------|--------------------| + # | 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 | | + check_ignored_arguments(..., ok_args = list("moment_match")) - if (!is.null(pit)) { - inform("'pit' specified so ignoring 'y','yrep','lw' if specified.") + # --------------------------------------------------------------------------- + # Internal helpers + # --------------------------------------------------------------------------- + + .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.") + } + } + + # --------------------------------------------------------------------------- + # Method-specific defaults and validation + # --------------------------------------------------------------------------- + + 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 applies only when `pit` is not already supplied and the + # test is POT or PIET. + 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) + } + ) + + # --------------------------------------------------------------------------- + # Compute PIT values + # --------------------------------------------------------------------------- + + 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))) + + # TODO: posterior::pareto_pit() from https://github.com/stan-dev/posterior/pull/435 - requires that PR merged + 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) + + # Warn about any ignored arguments. + 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) + } + + # --------------------------------------------------------------------------- + # Shared ECDF setup + # --------------------------------------------------------------------------- + + 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" + + # =========================================================================== + # Correlated method + # =========================================================================== + + if (method == "correlated") { + + # Compute the per-observation test statistics (sorted for Shapley values) + # and the combined Cauchy p-value. + # TODO: posterior::uniformity_test() from https://github.com/stan-dev/posterior/pull/435 - requires that PR merged + test_res <- posterior::uniformity_test(pit = pit, test = test) + p_value_CCT <- test_res$pvalue + pointwise_contrib <- test_res$pointwise + + # Validate gamma against computed Shapley values. + 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 + )) + } + + # Build the main ECDF data frame over a dense grid that includes pit values + # so step discontinuities are rendered exactly. + 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 + ) + + # Sorted pit data frame (used for highlighting suspicious points). + pit_sorted <- sort(pit) + df_pit <- tibble::tibble( + pit = pit_sorted, + ecdf_val = ecdf_pit_fn(pit_sorted) - plot_diff * pit_sorted + ) + + # --- Base plot ----------------------------------------------------------- + + 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) + + # --- Highlight suspicious regions ---------------------------------------- + + 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, ] + + # Highlight contiguous groups as coloured step segments. + 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 + ) + } + + # Highlight isolated suspicious points as dots. + 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 + ) + } + } + } + + # --- Annotation and axis scaling ----------------------------------------- + + 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) } + + # =========================================================================== + # Independent method + # =========================================================================== - n_obs <- length(pit) - gamma <- 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" + gamma_indep <- adjust_gamma(N = n_obs, K = K, prob = prob, + interpolate_adj = interpolate_adj) + + lims <- ecdf_intervals(gamma = gamma_indep, N = n_obs, K = K) + + # `lims` contains K + 1 elements (including the boundary at 0); drop it so + # lengths match `unit_interval`. + 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 = "PIT", y = y_label) + yaxis_ticks(FALSE) + scale_color_ppc() + bayesplot_theme_get() + + return(p) } @@ -836,7 +1128,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 c31242e8..5cd1b4ee 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 = "independent", + test = NULL, + gamma = NULL, + linewidth = NULL, + color = NULL, + help_text = NULL, + pareto_pit = NULL ) ppc_pit_ecdf_grouped( @@ -239,11 +246,45 @@ 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()} 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()}, 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()} 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()} 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()} when \code{method = "correlated"}, the line width of the ECDF +and highlighting points. Defaults to 0.3.} + +\item{color}{For \code{ppc_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_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_pit_ecdf()}, 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 555898d8..52aa4010 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 not normally 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 @@ -230,6 +268,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-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..deefd121 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 +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..8f4da5ef 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 +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..294aac36 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 +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..8127fd63 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 +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 9062515f..d4cef0e9 100644 --- a/tests/testthat/test-ppc-distributions.R +++ b/tests/testthat/test-ppc-distributions.R @@ -104,15 +104,59 @@ 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")) + + # Specify 'pit' directly expect_message( ppc_pit_ecdf_grouped(pit = runif(length(group)), group = group, interpolate_adj = FALSE), "'pit' specified" ) }) +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_freqpoly_grouped returns a ggplot object", { expect_gg(ppc_freqpoly_grouped(y, yrep[1:4, ], group)) expect_gg(ppc_freqpoly_grouped(y, yrep[1:4, ], group, @@ -412,6 +456,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) @@ -421,4 +466,294 @@ 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) +}) + +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) { + is.call(e) && deparse(e[[1]]) == "if" && + grepl("pareto_pit", deparse(e[[2]])) + })) + ]] <- 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_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 9381692f..08464262 100644 --- a/tests/testthat/test-ppc-loo.R +++ b/tests/testthat/test-ppc-loo.R @@ -91,7 +91,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, "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)) @@ -103,6 +103,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") @@ -121,7 +213,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)) @@ -136,7 +228,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") @@ -199,6 +290,44 @@ test_that("error if subset is bigger than num obs", { ) }) +test_that("ppc_loo_pit_ecdf works with pareto_pit method", { + 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" + )) + # prit -> pareto_pit should not be default (doesn't matter whether y, yrep, pit is provided) + # y, yrep + pot, piet -> pareto_pit + # pit -> no additional pareto_pit + # add in the docs that the pareto_pit on/off should usually not be touched by the user. Default is okay + # secondary step: ppcheck in brms -> cdf based pit +}) + # Visual tests ------------------------------------------------------------ @@ -299,6 +428,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") @@ -338,4 +546,259 @@ 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) { + is.call(e) && deparse(e[[1]]) == "if" && + grepl("pareto_pit", deparse(e[[2]])) + })) + ]] <- 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" + ) })