diff --git a/.gitignore b/.gitignore
index 695e6e77..1e896959 100644
--- a/.gitignore
+++ b/.gitignore
@@ -16,3 +16,17 @@ 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
\ 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/NAMESPACE b/NAMESPACE
index b8cf818b..abc7365e 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -27,6 +27,9 @@ S3method(nuts_params,stanreg)
S3method(parameter_names,array)
S3method(parameter_names,default)
S3method(parameter_names,matrix)
+S3method(pareto_pit,default)
+S3method(pareto_pit,draws_matrix)
+S3method(pareto_pit,rvar)
S3method(plot,bayesplot_grid)
S3method(plot,bayesplot_scheme)
S3method(pp_check,default)
@@ -55,6 +58,7 @@ export(example_y_data)
export(example_yrep_draws)
export(facet_bg)
export(facet_text)
+export(gpdfit)
export(grid_lines)
export(hline_0)
export(hline_at)
@@ -115,6 +119,8 @@ export(panel_bg)
export(param_glue)
export(param_range)
export(parcoord_style_np)
+export(pareto_pit)
+export(pgeneralized_pareto)
export(plot_bg)
export(pp_check)
export(ppc_bars)
@@ -190,6 +196,7 @@ export(ppd_stat_data)
export(ppd_stat_freqpoly)
export(ppd_stat_freqpoly_grouped)
export(ppd_stat_grouped)
+export(qgeneralized_pareto)
export(rhat)
export(scatter_style_np)
export(theme_default)
diff --git a/R/helpers-ppc.R b/R/helpers-ppc.R
index e268b315..9f00fb27 100644
--- a/R/helpers-ppc.R
+++ b/R/helpers-ppc.R
@@ -595,3 +595,598 @@ 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])
+
+# Dependence-aware uniformity tests -------------------------------------
+
+#' Compute Shapley values
+#'
+#' Calculates the average marginal contribution of players across
+#' all random arrival orders in a cooperative game.
+#' Used to provide a principled approach for quantifying
+#' point-specific influences in a way that reflects local miscalibration.
+#'
+#' @param x Numeric vector of Cauchy-transformed PIT values.
+#' @return Numeric vector of Shapley values with the same length as `x`.
+#' @noRd
+.compute_shapley_values <- function(x) {
+ n <- length(x)
+ if (n == 0) {
+ return(numeric(0))
+ }
+ if (n == 1) {
+ return(0)
+ }
+
+ # Harmonic number
+ # H_n = sum(1/i) for i = 1 to n
+ harmonic_number <- sum(1 / seq_len(n))
+
+ shapley_values <- numeric(n)
+ for (i in seq_len(n)) {
+ mean_others <- sum(x[-i]) / (n - 1)
+ # Applies the closed-form formula to assign player i their fair share.
+ shapley_values[i] <- (1 / n) * x[i] + ((harmonic_number - 1) / n) * (x[i] - mean_others)
+ }
+
+ return(shapley_values)
+}
+
+#' Pointwise Inverse-CDF Evaluation Tests Combination (PIET)
+#'
+#' Uniformity test with respect to any continuous distribution.
+#' H0: The value obtained via the inverse CDF transformation F^(-1)(x_i)
+#' follows the distribution of X under uniformity.
+#' HA: The p-value p_(i) provides evidence against uniformity.
+#'
+#' @param x Numeric vector of PIT values in `[0, 1]`.
+#' @return Numeric vector of p-values.
+#' @noRd
+.piet_test <- function(x) {
+ cdf_exp <- pexp(-log(x), rate = 1) # same as 1-x but numerically more stable
+ p_values <- 2 * pmin(cdf_exp, 1 - cdf_exp)
+
+ return(p_values)
+}
+
+#' Pointwise Order Tests Combination (POT)
+#'
+#' Uniformity test based on a beta distribution.
+#' H0: The i-th order statistic u_(i) follows a beta distribution
+#' under uniformity.
+#' HA: The p-value p_(i) provides evidence against uniformity
+#' at the i-th order statistic u_(i).
+#'
+#' @param x Numeric vector of PIT values in `[0, 1]`.
+#' @return Numeric vector of p-values.
+#' @noRd
+.pot_test <- function(x) {
+ n <- length(x)
+ # keep NA values instead of silent recycling via sort()
+ # TODO: Can PIT values be NAN at this point?
+ cdf_beta <- pbeta(sort(x, na.last = TRUE), 1:n, seq(n, 1, by = -1))
+ p_values <- 2 * pmin(cdf_beta, 1 - cdf_beta)
+
+ return(p_values)
+}
+
+#' Pointwise Rank-based Individual Tests Combination (PRIT)
+#'
+#' Uniformity test based on a binomial distribution.
+#' H0: The number of observations falling at or below x_i
+#' follows a binomial distribution under uniformity.
+#' HA: The p-value p_i provides evidence against uniformity.
+#'
+#' @param x Numeric vector of PIT values in `[0, 1]`.
+#' @return Numeric vector of p-values.
+#' @noRd
+.prit_test <- function(x) {
+ n <- length(x)
+ scaled_ecdf <- n * ecdf(x)(x)
+ probs1 <- pbinom(scaled_ecdf - 1, n, x)
+ probs2 <- pbinom(scaled_ecdf, n, x)
+ p_values <- 2 * pmin(1 - probs1, probs2)
+
+ return(p_values)
+}
+
+#' Truncated Cauchy combination test
+#'
+#' Combines dependent p-values using the Cauchy combination method.
+#' If truncate, only p-values less than 0.5 are included.
+#'
+#' @param x Numeric vector of p-values transformed to follow a standard
+#' Cauchy distribution.
+#' @param truncate Boolean; If TRUE only p-values less than 0.5 are
+#' included.
+#' @return p-value of the Cauchy combination method.
+#' @noRd
+.cauchy_combination_test <- function(x, truncate = NULL) {
+ if (truncate) {
+ idx <- which(x < 0.5)
+ if (length(idx) == 0) {
+ stop("Cannot compute truncated Cauchy combination test. ",
+ "No p-values below 0.5 found.")
+ }
+ 1 - pcauchy(mean(-qcauchy(x[idx])))
+ } else {
+ 1 - pcauchy(mean(-qcauchy(x)))
+ }
+}
+
+#' Compute Cauchy transformation
+#'
+#' Transforms PIT values to follow a standard Cauchy distribution.
+#'
+#' @param x Numeric vector of PIT values in `[0, 1]`.
+#' @return Numeric vector of Cauchy-transformed values.
+#' @noRd
+.compute_cauchy <- function(x) {
+ tan((0.5 - x) * pi)
+}
+
+# Support pareto_pit in ppc_loo_pit_ecdf and ppc_pit_ecdf --------------
+# TODO: Once this function is available via posterior, remove this internal implementation
+# and import from posterior.
+
+#' Pareto-smoothed probability integral transform
+#'
+#' Compute PIT values using the empirical CDF, then refine values in
+#' the tails by fitting a generalized Pareto distribution (GPD) to
+#' the tail draws. This gives smoother, more accurate PIT values in
+#' the tails where the ECDF is coarse, and avoids PIT values of 0 and 1.
+#' Due to use of generalized Pareto distribution CDF in tails, the
+#' PIT values are not anymore rank based and continuous uniformity
+#' test is appropriate.
+#'
+#' @name pareto_pit
+#'
+#' @param x (draws) A `posterior::draws_matrix` object or one coercible to a
+#' `posterior::draws_matrix` object, or an `posterior::rvar` object.
+#'
+#' @param y (observations) A 1D vector, or an array of dim(x), if x is `posterior::rvar`.
+#' Each element of `y` corresponds to a variable in `x`.
+#'
+#' @param weights A matrix of weights for each draw and variable. `weights`
+#' should have one column per variable in `x`, and `posterior::ndraws(x)` rows.
+#'
+#' @param log (logical) Are the weights passed already on the log scale? The
+#' default is `FALSE`, that is, expecting `weights` to be on the standard
+#' (non-log) scale.
+#'
+#' @param ndraws_tail (integer) Number of tail draws to use for GPD
+#' fitting. If `NULL` (the default), computed using `posterior::ps_tail_length()`.
+#'
+#' @template args-methods-dots
+#'
+#' @details The function first computes raw PIT values identically to
+#' `posterior::pit()` (including support for weighted draws). It then fits a
+#' GPD to both tails of the draws (using the same approach as
+#' `posterior::pareto_smooth()`) and replaces PIT values for observations falling in
+#' the tail regions:
+#'
+#' For a right-tail observation \eqn{y_i > c_R} (where \eqn{c_R} is
+#' the right-tail cutoff):
+#'
+#' \deqn{PIT(y_i) = 1 - p_{tail}(1 - F_{GPD}(y_i; c_R, \sigma_R, k_R))}
+#'
+#' For a left-tail observation \eqn{y_i < c_L}:
+#'
+#' \deqn{PIT(y_i) = p_{tail}(1 - F_{GPD}(-y_i; -c_L, \sigma_L, k_L))}
+#'
+#' where \eqn{p_{tail}} is the proportion of (weighted) mass in the
+#' tail.
+#'
+#' When (log-)weights in `weights` are provided, they are used for
+#' the raw PIT computation (as in `posterior::pit()`) and for GPD fit.
+#'
+#' @return A numeric vector of length `length(y)` containing the PIT values, or
+#' an array of shape `dim(y)`, if `x` is an `rvar`.
+#'
+#' @seealso `posterior::pit()` for the unsmoothed version, `posterior::pareto_smooth()` for
+#' Pareto smoothing of draws.
+#'
+#' @examples
+#' x <- posterior::example_draws()
+#' y <- rnorm(posterior::nvariables(x), 5, 5)
+#' pareto_pit(x, y)
+#'
+NULL
+
+#' @rdname pareto_pit
+#' @export
+pareto_pit <- function(x, y, ...) UseMethod("pareto_pit")
+
+#' @rdname pareto_pit
+#' @export
+pareto_pit.default <- function(x, y, weights = NULL, log = FALSE,
+ ndraws_tail = NULL, ...) {
+ x <- posterior::as_draws_matrix(x)
+ if (!is.null(weights)) {
+ weights <- posterior::as_draws_matrix(weights)
+ }
+ pareto_pit(x, y, weights = weights, log = log,
+ ndraws_tail = ndraws_tail, ...)
+}
+
+#' @rdname pareto_pit
+#' @export
+pareto_pit.draws_matrix <- function(x, y, weights = NULL, log = FALSE,
+ ndraws_tail = NULL, ...) {
+ y <- pareto_pit_validate_y(y, x)
+
+ # validate and normalize weights to log scale (same as pit.draws_matrix)
+ if (!is.null(weights)) {
+ weights <- sapply(seq_len(posterior::nvariables(x)), function(var_idx) {
+ posterior:::validate_weights(weights[, var_idx], x[, var_idx], log)
+ })
+ weights <- normalize_log_weights(weights)
+ }
+
+ ndraws <- posterior:::ndraws(x)
+
+ if (is.null(ndraws_tail)) {
+ ndraws_tail <- posterior::ps_tail_length(ndraws, 1)
+ } else {
+ ndraws_tail <- posterior:::as_one_integer(ndraws_tail)
+ }
+
+ # validate ndraws_tail once for all variables
+ gpd_ok <- !is.na(ndraws_tail) && ndraws_tail >= 5
+ if (gpd_ok && ndraws_tail > ndraws / 2) {
+ ndraws_tail <- floor(ndraws / 2)
+ }
+ if (gpd_ok && ndraws_tail >= ndraws) {
+ gpd_ok <- FALSE
+ }
+
+ # precompute tail indices (shared across all variables)
+ if (gpd_ok) {
+ tail_ids <- seq(ndraws - ndraws_tail + 1, ndraws)
+ }
+
+ pit_values <- vapply(seq_len(ncol(x)), function(j) {
+ draws <- x[, j]
+
+ # --- raw PIT (same logic as pit.draws_matrix) ---
+ sel_min <- draws < y[j]
+ if (!any(sel_min)) {
+ raw_pit <- 0
+ } else {
+ if (is.null(weights)) {
+ raw_pit <- mean(sel_min)
+ } else {
+ raw_pit <- exp(posterior:::log_sum_exp(weights[sel_min, j]))
+ }
+ }
+
+ sel_sup <- draws == y[j]
+ if (any(sel_sup)) {
+ if (is.null(weights)) {
+ pit_sup <- raw_pit + mean(sel_sup)
+ } else {
+ pit_sup <- raw_pit + exp(posterior:::log_sum_exp(weights[sel_sup, j]))
+ }
+ raw_pit <- runif(1, raw_pit, pit_sup)
+ }
+
+ # --- GPD tail refinement ---
+ if (!gpd_ok || posterior:::should_return_NA(draws)) {
+ return(raw_pit)
+ }
+
+ # sort draws and carry weights along
+ ord <- sort.int(draws, index.return = TRUE)
+ sorted <- ord$x
+ log_wt_sorted <- if (!is.null(weights)) weights[ord$ix, j] else NULL
+
+ # tail proportion: sum of (normalized) weights in the tail
+ if (!is.null(log_wt_sorted)) {
+ tail_proportion <- exp(posterior:::log_sum_exp(log_wt_sorted[tail_ids]))
+ } else {
+ tail_proportion <- ndraws_tail / ndraws
+ }
+
+ # --- right tail ---
+ right_replaced <- FALSE
+ right_tail <- sorted[tail_ids]
+ if (!posterior:::is_constant(right_tail)) {
+ right_cutoff <- sorted[min(tail_ids) - 1]
+ if (right_cutoff == right_tail[1]) {
+ right_cutoff <- right_cutoff - .Machine$double.eps
+ }
+ right_tail_wt <- if (!is.null(log_wt_sorted)) {
+ exp(log_wt_sorted[tail_ids])
+ } else {
+ NULL
+ }
+ right_fit <- gpdfit(right_tail - right_cutoff, sort_x = FALSE,
+ weights = right_tail_wt)
+ if (is.finite(right_fit$k) && !is.na(right_fit$sigma)) {
+ if (y[j] > right_cutoff) {
+ gpd_cdf <- pgeneralized_pareto(
+ y[j], mu = right_cutoff, sigma = right_fit$sigma, k = right_fit$k
+ )
+ raw_pit <- 1 - tail_proportion * (1 - gpd_cdf)
+ right_replaced <- TRUE
+ }
+ }
+ }
+
+ # --- left tail (negate trick, same as ps_tail) ---
+ if (!right_replaced) {
+ left_draws <- -draws
+ left_ord <- sort.int(left_draws, index.return = TRUE)
+ left_sorted <- left_ord$x
+ log_wt_left_sorted <- if (!is.null(weights)) weights[left_ord$ix, j] else NULL
+
+ left_tail <- left_sorted[tail_ids]
+ if (!posterior:::is_constant(left_tail)) {
+ left_cutoff <- left_sorted[min(tail_ids) - 1]
+ if (left_cutoff == left_tail[1]) {
+ left_cutoff <- left_cutoff - .Machine$double.eps
+ }
+ left_tail_wt <- if (!is.null(log_wt_left_sorted)) {
+ exp(log_wt_left_sorted[tail_ids])
+ } else {
+ NULL
+ }
+ left_fit <- gpdfit(left_tail - left_cutoff, sort_x = FALSE,
+ weights = left_tail_wt)
+ if (is.finite(left_fit$k) && !is.na(left_fit$sigma)) {
+ if (-y[j] > left_cutoff) {
+ gpd_cdf <- pgeneralized_pareto(
+ -y[j], mu = left_cutoff, sigma = left_fit$sigma, k = left_fit$k
+ )
+ if (!is.null(log_wt_left_sorted)) {
+ left_tail_proportion <- exp(posterior:::log_sum_exp(log_wt_left_sorted[tail_ids]))
+ } else {
+ left_tail_proportion <- tail_proportion
+ }
+ raw_pit <- left_tail_proportion * (1 - gpd_cdf)
+ }
+ }
+ }
+ }
+
+ raw_pit
+ }, FUN.VALUE = 1.0)
+
+ min_tail_prob <- 1/ndraws/1e4
+ pit_values <- pmin(pmax(pit_values, min_tail_prob), 1-min_tail_prob)
+ setNames(pit_values, posterior::variables(x))
+}
+
+#' @rdname pareto_pit
+#' @export
+pareto_pit.rvar <- function(x, y, weights = NULL, log = FALSE,
+ ndraws_tail = NULL, ...) {
+ y <- pareto_pit_validate_y(y, x)
+ out <- array(
+ data = pareto_pit(
+ x = posterior::as_draws_matrix(c(x)),
+ y = c(y),
+ weights = weights,
+ log = log,
+ ndraws_tail = ndraws_tail,
+ ...
+ ),
+ dim = dim(x),
+ dimnames = dimnames(x)
+ )
+ out
+}
+
+# internal ----------------------------------------------------------------
+
+pareto_pit_validate_y <- function(y, x = NULL) {
+ if (!is.numeric(y)) {
+ posterior:::stop_no_call("`y` must be numeric.")
+ }
+ if (anyNA(y)) {
+ posterior:::stop_no_call("NAs not allowed in `y`.")
+ }
+ if (posterior::is_rvar(x)) {
+ if (length(x) != length(y) || any(dim(y) != dim(x))) {
+ posterior:::stop_no_call("`dim(y)` must match `dim(x)`.")
+ }
+ } else if (posterior::is_draws(x)) {
+ if (!is.vector(y, mode = "numeric") || length(y) != posterior::nvariables(x)) {
+ posterior:::stop_no_call("`y` must be a vector of length `nvariables(x)`.")
+ }
+ }
+ y
+}
+
+normalize_log_weights <- function(log_weights) {
+ apply(log_weights, 2, function(col) col - posterior:::log_sum_exp(col))
+}
+
+
+#' Quantile function for the generalized Pareto distribution
+#'
+#' Computes the quantile function for a generalized Pareto distribution
+#' with location `mu`, scale `sigma`, and shape `k`.
+#'
+#' @param p Numeric vector of probabilities.
+#' @param mu Location parameter.
+#' @param sigma Scale parameter (must be positive).
+#' @param k Shape parameter.
+#' @param lower.tail Logical; if `TRUE` (default), probabilities are `P[X <= x]`.
+#' @param log.p Logical; if `TRUE`, probabilities are given as `log(p)`.
+#' @return A numeric vector of quantiles.
+#' @keywords internal
+#' @export
+#' @examples
+#' qgeneralized_pareto(p = c(0.1, 0.5, 0.9), mu = 0, sigma = 1, k = 0.2)
+qgeneralized_pareto <- function(p, mu = 0, sigma = 1, k = 0, lower.tail = TRUE, log.p = FALSE) {
+ stopifnot(length(mu) == 1 && length(sigma) == 1 && length(k) == 1)
+ if (is.na(sigma) || sigma <= 0) {
+ return(rep(NaN, length(p)))
+ }
+ if (log.p) {
+ p <- exp(p)
+ }
+ if (!lower.tail) {
+ p <- 1 - p
+ }
+ if (k == 0) {
+ q <- mu - sigma * log1p(-p)
+ } else {
+ q <- mu + sigma * expm1(-k * log1p(-p)) / k
+ }
+ q
+}
+
+#' Distribution function for the generalized Pareto distribution
+#'
+#' Computes the cumulative distribution function (CDF) for a generalized
+#' Pareto distribution with location `mu`, scale `sigma`, and shape `k`.
+#'
+#' @param q Numeric vector of quantiles.
+#' @param mu Location parameter.
+#' @param sigma Scale parameter (must be positive).
+#' @param k Shape parameter.
+#' @param lower.tail Logical; if `TRUE` (default), probabilities are `P[X <= x]`.
+#' @param log.p Logical; if `TRUE`, probabilities are returned as `log(p)`.
+#' @return A numeric vector of probabilities.
+#' @keywords internal
+#' @export
+#' @examples
+#' pgeneralized_pareto(q = c(1, 2, 5), mu = 0, sigma = 1, k = 0.2)
+pgeneralized_pareto <- function(q, mu = 0, sigma = 1, k = 0, lower.tail = TRUE, log.p = FALSE) {
+ stopifnot(length(mu) == 1 && length(sigma) == 1 && length(k) == 1)
+ if (is.na(sigma) || sigma <= 0) {
+ return(rep(NaN, length(q)))
+ }
+ z <- (q - mu) / sigma
+ if (abs(k) < 1e-15) {
+ # for very small values of indistinguishable in floating point accuracy from the case k=0
+ p <- -expm1(-z)
+ } else {
+ # pmax handles values outside the support
+ p <- -expm1(log1p(pmax(k * z, -1)) / -k)
+ }
+ # force to [0, 1] for values outside the support
+ p <- pmin(pmax(p, 0), 1)
+ if (!lower.tail) {
+ p <- 1 - p
+ }
+ if (log.p) {
+ p <- log(p)
+ }
+ p
+}
+
+#' Estimate parameters of the Generalized Pareto distribution
+#'
+#' Given a sample \eqn{x}, Estimate the parameters \eqn{k} and
+#' \eqn{\sigma} of the generalized Pareto distribution (GPD), assuming
+#' the location parameter is 0. By default the fit uses a prior for
+#' \eqn{k} (this is in addition to the prior described by Zhang and
+#' Stephens, 2009), which will stabilize estimates for very small
+#' sample sizes (and low effective sample sizes in the case of MCMC
+#' samples). The weakly informative prior is a Gaussian prior centered
+#' at 0.5 (see details in Vehtari et al., 2024). This is used
+#' internally but is exported for use by other packages.
+#' @family helper-functions
+#' @param x A numeric vector. The sample from which to estimate the
+#' parameters.
+#' @param wip Logical indicating whether to adjust \eqn{k} based on a
+#' weakly informative Gaussian prior centered on 0.5. Defaults to
+#' `TRUE`.
+#' @param min_grid_pts The minimum number of grid points used in the
+#' fitting algorithm. The actual number used is `min_grid_pts +
+#' floor(sqrt(length(x)))`.
+#' @param sort_x If `TRUE` (the default), the first step in the
+#' fitting algorithm is to sort the elements of `x`. If `x` is
+#' already sorted in ascending order then `sort_x` can be set to
+#' `FALSE` to skip the initial sorting step.
+#' @param weights An optional numeric vector of positive weights the same
+#' length as `x`. If `NULL` (the default), all observations are
+#' weighted equally and the result is identical to the unweighted fit.
+#' Weights are normalized internally to sum to `length(x)`.
+#' @return A named list with components `k` and `sigma`.
+#'
+#' @details Here the parameter \eqn{k} is the negative of \eqn{k} in Zhang &
+#' Stephens (2009).
+#'
+#'
+#' @references
+#' Zhang, J., and Stephens, M. A. (2009). A new and efficient estimation method
+#' for the generalized Pareto distribution. *Technometrics* **51**, 316-325.
+#'
+#' @keywords internal
+#' @export
+gpdfit <- function(x, wip = TRUE, min_grid_pts = 30, sort_x = TRUE,
+ weights = NULL) {
+ # see section 4 of Zhang and Stephens (2009)
+ if (sort_x) {
+ if (!is.null(weights)) {
+ ord <- sort.int(x, index.return = TRUE)
+ x <- ord$x
+ weights <- weights[ord$ix]
+ } else {
+ x <- sort.int(x)
+ }
+ }
+ N <- length(x)
+
+ # normalize weights to sum to N so the log-likelihood scale is preserved
+ if (!is.null(weights)) {
+ weights <- weights / sum(weights) * N
+ }
+
+ prior <- 3
+ M <- min_grid_pts + floor(sqrt(N))
+ jj <- seq_len(M)
+ xstar <- x[floor(N / 4 + 0.5)] # first quartile of sample
+ if (xstar > x[1]) {
+ # first quantile is bigger than the minimum
+ theta <- 1 / x[N] + (1 - sqrt(M / (jj - 0.5))) / prior / xstar
+
+ # log1p(-theta %o% x) is M x N matrix
+ log1p_mat <- log1p(-theta %o% x)
+
+ if (!is.null(weights)) {
+ # weighted mean across observations for each theta value
+ k <- drop(log1p_mat %*% weights) / N
+ } else {
+ k <- rowMeans(log1p_mat)
+ }
+
+ l_theta <- N * (log(-theta / k) - k - 1) # profile log-lik
+ w_theta <- exp(l_theta - posterior:::log_sum_exp(l_theta)) # normalize
+ theta_hat <- sum(theta * w_theta)
+
+ if (!is.null(weights)) {
+ k_hat <- sum(weights * log1p(-theta_hat * x)) / N
+ } else {
+ k_hat <- mean.default(log1p(-theta_hat * x))
+ }
+ sigma_hat <- -k_hat / theta_hat
+
+ # adjust k_hat based on weakly informative prior, Gaussian centered on 0.5.
+ # this stabilizes estimates for very small Monte Carlo sample sizes and low ess
+ # (see Vehtari et al., 2024 for details)
+ if (wip) {
+ k_hat <- (k_hat * N + 0.5 * 10) / (N + 10)
+ }
+ if (is.na(k_hat)) {
+ k_hat <- Inf
+ sigma_hat <- NaN
+ }
+ } else {
+ # first quantile is not bigger than the minimum, which indicates
+ # that the distribution is far from a generalized Pareto
+ # distribution
+ k_hat <- NA
+ sigma_hat <- NA
+ }
+
+ list(k = k_hat, sigma = sigma_hat)
+}
+
+# helper function for formatting p-value
+fmt_p <- function(x) {
+ dplyr::if_else(x < 0.0005, "0.000", as.character(round(signif(x, 2) + 1e-10, 3)))
+}
\ No newline at end of file
diff --git a/R/ppc-distributions.R b/R/ppc-distributions.R
index d86a8b4c..c82fefa8 100644
--- a/R/ppc-distributions.R
+++ b/R/ppc-distributions.R
@@ -648,6 +648,36 @@ 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.
#' @rdname PPC-distributions
#'
ppc_pit_ecdf <- function(y,
@@ -657,59 +687,317 @@ 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))
+
+ pit <- 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.
+ test_fn <- switch(test,
+ POT = .pot_test,
+ PIET = .piet_test,
+ PRIT = .prit_test
+ )
+
+ pit_sorted <- sort(pit)
+ std_cauchy_values <- .compute_cauchy(test_fn(pit_sorted))
+ p_value_CCT <- .cauchy_combination_test(
+ test_fn(pit),
+ truncate = test != "PIET"
+ )
+ pointwise_contrib <- .compute_shapley_values(std_cauchy_values)
+
+ # 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).
+ 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..8a16458f 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,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 not normally 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,
@@ -409,63 +436,333 @@ 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)))
+
+ pit <- 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.
+ test_fn <- switch(test,
+ POT = .pot_test,
+ PIET = .piet_test,
+ PRIT = .prit_test
+ )
+
+ pit_sorted <- sort(pit)
+ std_cauchy_values <- .compute_cauchy(test_fn(pit_sorted))
+ p_value_CCT <- .cauchy_combination_test(
+ test_fn(pit),
+ truncate = test != "PIET"
+ )
+ pointwise_contrib <- .compute_shapley_values(std_cauchy_values)
+
+ # 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).
+ 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 +1133,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/man/gpdfit.Rd b/man/gpdfit.Rd
new file mode 100644
index 00000000..e0567025
--- /dev/null
+++ b/man/gpdfit.Rd
@@ -0,0 +1,53 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/helpers-ppc.R
+\name{gpdfit}
+\alias{gpdfit}
+\title{Estimate parameters of the Generalized Pareto distribution}
+\usage{
+gpdfit(x, wip = TRUE, min_grid_pts = 30, sort_x = TRUE, weights = NULL)
+}
+\arguments{
+\item{x}{A numeric vector. The sample from which to estimate the
+parameters.}
+
+\item{wip}{Logical indicating whether to adjust \eqn{k} based on a
+weakly informative Gaussian prior centered on 0.5. Defaults to
+\code{TRUE}.}
+
+\item{min_grid_pts}{The minimum number of grid points used in the
+fitting algorithm. The actual number used is \code{min_grid_pts + floor(sqrt(length(x)))}.}
+
+\item{sort_x}{If \code{TRUE} (the default), the first step in the
+fitting algorithm is to sort the elements of \code{x}. If \code{x} is
+already sorted in ascending order then \code{sort_x} can be set to
+\code{FALSE} to skip the initial sorting step.}
+
+\item{weights}{An optional numeric vector of positive weights the same
+length as \code{x}. If \code{NULL} (the default), all observations are
+weighted equally and the result is identical to the unweighted fit.
+Weights are normalized internally to sum to \code{length(x)}.}
+}
+\value{
+A named list with components \code{k} and \code{sigma}.
+}
+\description{
+Given a sample \eqn{x}, Estimate the parameters \eqn{k} and
+\eqn{\sigma} of the generalized Pareto distribution (GPD), assuming
+the location parameter is 0. By default the fit uses a prior for
+\eqn{k} (this is in addition to the prior described by Zhang and
+Stephens, 2009), which will stabilize estimates for very small
+sample sizes (and low effective sample sizes in the case of MCMC
+samples). The weakly informative prior is a Gaussian prior centered
+at 0.5 (see details in Vehtari et al., 2024). This is used
+internally but is exported for use by other packages.
+}
+\details{
+Here the parameter \eqn{k} is the negative of \eqn{k} in Zhang &
+Stephens (2009).
+}
+\references{
+Zhang, J., and Stephens, M. A. (2009). A new and efficient estimation method
+for the generalized Pareto distribution. \emph{Technometrics} \strong{51}, 316-325.
+}
+\concept{helper-functions}
+\keyword{internal}
diff --git a/man/pareto_pit.Rd b/man/pareto_pit.Rd
new file mode 100644
index 00000000..b99bfc68
--- /dev/null
+++ b/man/pareto_pit.Rd
@@ -0,0 +1,81 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/helpers-ppc.R
+\name{pareto_pit}
+\alias{pareto_pit}
+\alias{pareto_pit.default}
+\alias{pareto_pit.draws_matrix}
+\alias{pareto_pit.rvar}
+\title{Pareto-smoothed probability integral transform}
+\usage{
+pareto_pit(x, y, ...)
+
+\method{pareto_pit}{default}(x, y, weights = NULL, log = FALSE, ndraws_tail = NULL, ...)
+
+\method{pareto_pit}{draws_matrix}(x, y, weights = NULL, log = FALSE, ndraws_tail = NULL, ...)
+
+\method{pareto_pit}{rvar}(x, y, weights = NULL, log = FALSE, ndraws_tail = NULL, ...)
+}
+\arguments{
+\item{x}{(draws) A \code{posterior::draws_matrix} object or one coercible to a
+\code{posterior::draws_matrix} object, or an \code{posterior::rvar} object.}
+
+\item{y}{(observations) A 1D vector, or an array of dim(x), if x is \code{posterior::rvar}.
+Each element of \code{y} corresponds to a variable in \code{x}.}
+
+\item{...}{Arguments passed to individual methods (if applicable).}
+
+\item{weights}{A matrix of weights for each draw and variable. \code{weights}
+should have one column per variable in \code{x}, and \code{posterior::ndraws(x)} rows.}
+
+\item{log}{(logical) Are the weights passed already on the log scale? The
+default is \code{FALSE}, that is, expecting \code{weights} to be on the standard
+(non-log) scale.}
+
+\item{ndraws_tail}{(integer) Number of tail draws to use for GPD
+fitting. If \code{NULL} (the default), computed using \code{posterior::ps_tail_length()}.}
+}
+\value{
+A numeric vector of length \code{length(y)} containing the PIT values, or
+an array of shape \code{dim(y)}, if \code{x} is an \code{rvar}.
+}
+\description{
+Compute PIT values using the empirical CDF, then refine values in
+the tails by fitting a generalized Pareto distribution (GPD) to
+the tail draws. This gives smoother, more accurate PIT values in
+the tails where the ECDF is coarse, and avoids PIT values of 0 and 1.
+Due to use of generalized Pareto distribution CDF in tails, the
+PIT values are not anymore rank based and continuous uniformity
+test is appropriate.
+}
+\details{
+The function first computes raw PIT values identically to
+\code{posterior::pit()} (including support for weighted draws). It then fits a
+GPD to both tails of the draws (using the same approach as
+\code{posterior::pareto_smooth()}) and replaces PIT values for observations falling in
+the tail regions:
+
+For a right-tail observation \eqn{y_i > c_R} (where \eqn{c_R} is
+the right-tail cutoff):
+
+\deqn{PIT(y_i) = 1 - p_{tail}(1 - F_{GPD}(y_i; c_R, \sigma_R, k_R))}
+
+For a left-tail observation \eqn{y_i < c_L}:
+
+\deqn{PIT(y_i) = p_{tail}(1 - F_{GPD}(-y_i; -c_L, \sigma_L, k_L))}
+
+where \eqn{p_{tail}} is the proportion of (weighted) mass in the
+tail.
+
+When (log-)weights in \code{weights} are provided, they are used for
+the raw PIT computation (as in \code{posterior::pit()}) and for GPD fit.
+}
+\examples{
+x <- posterior::example_draws()
+y <- rnorm(posterior::nvariables(x), 5, 5)
+pareto_pit(x, y)
+
+}
+\seealso{
+\code{posterior::pit()} for the unsmoothed version, \code{posterior::pareto_smooth()} for
+Pareto smoothing of draws.
+}
diff --git a/man/pgeneralized_pareto.Rd b/man/pgeneralized_pareto.Rd
new file mode 100644
index 00000000..f723656f
--- /dev/null
+++ b/man/pgeneralized_pareto.Rd
@@ -0,0 +1,39 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/helpers-ppc.R
+\name{pgeneralized_pareto}
+\alias{pgeneralized_pareto}
+\title{Distribution function for the generalized Pareto distribution}
+\usage{
+pgeneralized_pareto(
+ q,
+ mu = 0,
+ sigma = 1,
+ k = 0,
+ lower.tail = TRUE,
+ log.p = FALSE
+)
+}
+\arguments{
+\item{q}{Numeric vector of quantiles.}
+
+\item{mu}{Location parameter.}
+
+\item{sigma}{Scale parameter (must be positive).}
+
+\item{k}{Shape parameter.}
+
+\item{lower.tail}{Logical; if \code{TRUE} (default), probabilities are \code{P[X <= x]}.}
+
+\item{log.p}{Logical; if \code{TRUE}, probabilities are returned as \code{log(p)}.}
+}
+\value{
+A numeric vector of probabilities.
+}
+\description{
+Computes the cumulative distribution function (CDF) for a generalized
+Pareto distribution with location \code{mu}, scale \code{sigma}, and shape \code{k}.
+}
+\examples{
+pgeneralized_pareto(q = c(1, 2, 5), mu = 0, sigma = 1, k = 0.2)
+}
+\keyword{internal}
diff --git a/man/qgeneralized_pareto.Rd b/man/qgeneralized_pareto.Rd
new file mode 100644
index 00000000..c589071f
--- /dev/null
+++ b/man/qgeneralized_pareto.Rd
@@ -0,0 +1,39 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/helpers-ppc.R
+\name{qgeneralized_pareto}
+\alias{qgeneralized_pareto}
+\title{Quantile function for the generalized Pareto distribution}
+\usage{
+qgeneralized_pareto(
+ p,
+ mu = 0,
+ sigma = 1,
+ k = 0,
+ lower.tail = TRUE,
+ log.p = FALSE
+)
+}
+\arguments{
+\item{p}{Numeric vector of probabilities.}
+
+\item{mu}{Location parameter.}
+
+\item{sigma}{Scale parameter (must be positive).}
+
+\item{k}{Shape parameter.}
+
+\item{lower.tail}{Logical; if \code{TRUE} (default), probabilities are \code{P[X <= x]}.}
+
+\item{log.p}{Logical; if \code{TRUE}, probabilities are given as \code{log(p)}.}
+}
+\value{
+A numeric vector of quantiles.
+}
+\description{
+Computes the quantile function for a generalized Pareto distribution
+with location \code{mu}, scale \code{sigma}, and shape \code{k}.
+}
+\examples{
+qgeneralized_pareto(p = c(0.1, 0.5, 0.9), mu = 0, sigma = 1, k = 0.2)
+}
+\keyword{internal}
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 @@
+
+
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 @@
+
+
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 @@
+
+
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 @@
+
+
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 @@
+
+
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 @@
+
+
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 @@
+
+
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 @@
+
+
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 @@
+
+
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 @@
+
+
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 @@
+
+
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 @@
+
+
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 @@
+
+
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 @@
+
+
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 @@
+
+
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 @@
+
+
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 @@
+
+
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 @@
+
+
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 @@
+
+
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 @@
+
+
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 @@
+
+
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 @@
+
+
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 @@
+
+
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 @@
+
+
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 @@
+
+
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..67bc194a 100644
--- a/tests/testthat/test-helpers-ppc.R
+++ b/tests/testthat/test-helpers-ppc.R
@@ -124,3 +124,179 @@ test_that("ecdf_intervals returns right dimensions and values", {
expect_equal(min(lims$lower), 0)
expect_equal(max(lims$lower), 100)
})
+
+# Tests for the dependence-aware uniformity tests ------------------------------
+test_that("piet_test computes correct p-values", {
+ x <- c(0.1, 0.25, 0.5, 0.75, 0.9)
+ # Equivalent to 2 * min(1-x, x)
+ expected <- c(0.2, 0.5, 1.0, 0.5, 0.2)
+
+ expect_equal(.piet_test(x), expected, tolerance = 1e-7)
+})
+
+test_that("piet_test handles boundary values of 0 and 1", {
+ x <- c(0, 1)
+ expected <- c(0, 0)
+
+ expect_equal(.piet_test(x), expected)
+})
+
+test_that("piet_test handles extreme values stably", {
+ # Testing values very close to 0 and 1
+ x <- c(1e-17, 1 - 1e-17)
+ expected <- c(2e-17, 2e-17)
+
+ # Tolerance needs to be adjusted for very small numbers
+ expect_equal(.piet_test(x), expected, tolerance = 1e-16)
+})
+
+test_that("piet_test handles NA, NaN, and empty inputs correctly", {
+ # NA and NaN propagation
+ x_na <- c(0.5, NA, NaN)
+ res_na <- .piet_test(x_na)
+
+ expect_equal(res_na[1], 1.0)
+ expect_true(is.na(res_na[2]))
+ expect_true(is.nan(res_na[3]))
+})
+
+test_that("pot_test calculates correct p-values", {
+ # Manually calculating expected values for x = c(0.2, 0.8), n = 2.
+ # Note: If X ~ Beta(1, b) then X ~ Kumaraswamy(1, b) with CDF 1 - (1 - x)^b
+ # and if X ~ Beta(a, 1) then X ~ Kumaraswamy(a, 1) with CDF to x^a
+ # Beta(1, 2) CDF at 0.2 is 1 - (1 - 0.2)^2 = 0.36
+ # Beta(2, 1) CDF at 0.8 is 0.8^2 = 0.64
+ # p-values: 2 * min(0.36, 1-0.36) = 0.72; 2 * min(0.64, 1-0.64) = 0.72
+ x <- c(0.8, 0.2)
+ expected <- c(0.72, 0.72)
+
+ expect_equal(.pot_test(x), expected)
+})
+
+test_that("pot_test handles boundary values correctly", {
+ x <- c(0, 1)
+ expected <- c(0, 0)
+
+ expect_equal(.pot_test(x), expected)
+})
+
+test_that("pot_test bounds p-values between 0 and 1 for extreme out-of-bounds inputs", {
+ # pbeta handles values outside [0, 1] by returning 0
+ x <- c(-0.5, 1.5)
+ expected <- c(0, 0)
+
+ expect_equal(.pot_test(x), expected)
+})
+
+test_that("pot_test handles NAs", {
+ x_na <- c(0.5, NA) # resulting in a = 2, b = 1
+ # Beta(2, 1) at 0.5 = 0.25 -> 2 * min(0.25, 0.75) = 0.5
+ expected <- c(0.5, NA)
+
+ expect_equal(.pot_test(x_na), expected)
+})
+
+test_that("prit_test computes correct p-values", {
+ # Let n = 2, x = c(0.5, 0.5)
+ # scaled_ecdf = 2 * c(1, 1) = c(2, 2)
+ # probs1 = pbinom(1, 2, 0.5) = 0.75
+ # probs2 = pbinom(2, 2, 0.5) = 1.00
+ # p_val = 2 * min(1 - 0.75, 1.00) = 2 * 0.25 = 0.5
+
+ x <- c(0.5, 0.5)
+ expect_equal(.prit_test(x), c(0.5, 0.5))
+})
+
+# Test for computation of Shapley values -------------------------------------
+
+test_that("compute_shapley_values handles empty vector and single element", {
+ result_empty <- .compute_shapley_values(numeric(0))
+ result_single <- .compute_shapley_values(5)
+
+ expect_equal(result_empty, numeric(0))
+ expect_equal(length(result_empty), 0)
+ expect_equal(result_single, 0)
+ expect_equal(length(result_single), 1)
+})
+
+test_that("compute_shapley_values for simple case", {
+ x <- c(1, 2)
+ result <- .compute_shapley_values(x)
+
+ # Manual calculation for n=2:
+ # harmonic_number = 1 + 1/2 = 1.5
+ # For i=1: mean_others = 2/1 = 2
+ # shapley[1] = (1/2)*1 + ((1.5-1)/2)*(1-2) = 0.5 - 0.25 = 0.25
+ # For i=2: mean_others = 1/1 = 1
+ # shapley[2] = (1/2)*2 + ((1.5-1)/2)*(2-1) = 1 + 0.25 = 1.25
+
+ expected <- c(0.25, 1.25)
+ expect_equal(result, expected, tolerance = 1e-10)
+ expect_equal(length(result), 2)
+})
+
+test_that("compute_shapley_values handles mixed input values", {
+ x <- c(-0.2, 0, 2, 3.1, 4.2)
+ result <- .compute_shapley_values(x)
+ expect_equal(length(result), 5)
+ expect_true(all(is.finite(result)))
+})
+
+# Test for (truncated) Cauchy combination test ----------------------------------
+test_that("cauchy_combination_test handles truncate = FALSE", {
+ # avg = mean(-qcauchy(x))
+ # result = 1 - pcauchy(avg)
+
+ x <- c(0.1, 0.2, 0.3)
+ result <- .cauchy_combination_test(x, truncate = FALSE)
+ expected <- 1 - pcauchy(mean(-qcauchy(x)))
+
+ expect_equal(result, expected, tolerance = 1e-10)
+ expect_true(is.finite(result))
+ expect_true(result >= 0 && result <= 1)
+})
+
+test_that("cauchy_combination_test handles truncate = TRUE", {
+ # avg = mean(-qcauchy(x))
+ # result = 1 - pcauchy(avg)
+
+ x <- c(0.1, 0.2, 0.3, 0.4, 0.7, 0.8)
+ result <- .cauchy_combination_test(x, truncate = TRUE)
+ expected <- 1 - pcauchy(mean(-qcauchy(c(0.1, 0.2, 0.3, 0.4))))
+
+ expect_equal(result, expected, tolerance = 1e-10)
+ expect_true(is.finite(result))
+ expect_true(result >= 0 && result <= 1)
+})
+
+test_that("cauchy_combination_test handles boundary values", {
+ # x = 0: -qdf(0) = Inf and cdf(Inf) = 1 -> 1 - 1 = 0
+ # x = 1: -qdf(1) = -Inf and cdf(-Inf) = 0 -> 1 - 0 = 1
+
+ expect_equal(.cauchy_combination_test(0, truncate = FALSE), 0)
+ expect_equal(.cauchy_combination_test(1, truncate = FALSE), 1)
+ expect_true(is.nan(.cauchy_combination_test(c(0, 1), truncate = FALSE)))
+ # TODO: if 1 included in vector, CCT will always evaluate to 0
+ # as the mean evaluates to Inf and 1 - cdf(Inf) = 1 - 1 = 0
+ expect_equal(.cauchy_combination_test(c(0, 0.3, 0.4, 1), truncate = TRUE), 0)
+})
+
+# Test for compute_cauchy -----------------------------------------------------
+
+test_that("compute_cauchy computes correct transformations", {
+ # For x = 0.5: tan((0.5 - 0.5) * pi) = tan(0) = 0
+ expect_equal(.compute_cauchy(0.5), 0, tolerance = 1e-10)
+
+ # For x = 0.25: tan((0.5 - 0.25) * pi) = tan(0.25 * pi) = 1
+ expect_equal(.compute_cauchy(0.25), 1, tolerance = 1e-10)
+
+ # For x = 0.75: tan((0.5 - 0.75) * pi) = tan(-0.25 * pi) = -1
+ expect_equal(.compute_cauchy(0.75), -1, tolerance = 1e-10)
+})
+
+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..edda6cfb 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
+
+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 <- 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)
+ }
+})
+
+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")
+
+ # | 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..78bed29b 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
+
+ppc_loo_pit_ecdf_patched <- ppc_loo_pit_ecdf # copy
+
+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 <- 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)
+ }
})
+
+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")
+
+ # | 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"
+ )
+})
\ No newline at end of file