From 763f8fdd17ee4e947a8ec1621781ac6fbb3ffa53 Mon Sep 17 00:00:00 2001 From: Malcolm Barrett Date: Fri, 21 Aug 2026 18:59:03 +0200 Subject: [PATCH 1/7] add function-valued interventions to check_edp() values --- R/check-edp.R | 381 ++++++++++- R/plot-edp.R | 39 +- man/autoplot.edp_result.Rd | 4 +- man/check_edp.Rd | 29 +- tests/testthat/_snaps/check-edp.md | 108 +++ .../edp-ecdf-by-intervention-value.svg | 46 +- .../edp-histogram-by-intervention-value.svg | 6 +- .../positivity-check-panel.svg | 6 +- tests/testthat/test-check-edp.R | 633 +++++++++++++++++- 9 files changed, 1168 insertions(+), 84 deletions(-) diff --git a/R/check-edp.R b/R/check-edp.R index c0d1a7d..e49f3de 100644 --- a/R/check-edp.R +++ b/R/check-edp.R @@ -82,10 +82,17 @@ edp_result <- new_class( #' estimator variant, selected with tidyselect. Defaults to `.covariates`. #' @param .treatment_covariates The treatment-model covariate columns for the #' estimator variant, selected with tidyselect. Defaults to `.covariates`. -#' @param values The intervention values \eqn{a^*}. `NULL` (the default) uses -#' every level of a factor exposure and every observed value of a binary or -#' character exposure, and the deciles of the observed exposure for continuous -#' exposures. +#' @param values The interventions. `NULL` (the default) uses every level of a +#' factor exposure and every observed value of a binary or character exposure, +#' and the deciles of the observed exposure for continuous exposures. A vector +#' holds one intervention value \eqn{a^*} per element. A list holds the same +#' single values and, for a continuous exposure, functions of the observed +#' exposure such as `list(3, 5, lower = function(.x) .x - 5)`. A function is +#' called once on the whole exposure vector and returns one intervened +#' exposure \eqn{d(a_i)} per observation, or a single value used for every +#' observation. Each function needs a name, which labels the intervention; an +#' unnamed value is labeled by the value itself, and every label in a list +#' must be distinct. #' @param variant One of `"data"` (report a single `edp` per point, the default) #' or `"estimator"` (report `edp_outcome`, `edp_treatment`, and #' `ideal_weight`). @@ -106,14 +113,16 @@ edp_result <- new_class( #' nothing of the column. #' #' @return An `edp_result` object, an S7 subclass of [positivity_diagnostic]. -#' Its `@results` tibble has one row per observation and intervention value. -#' The data variant carries columns `.id` (the observation row index), `value` -#' (the intervention value \eqn{a^*}), and `edp`. The estimator variant -#' replaces `edp` with `edp_outcome`, `edp_treatment`, and `ideal_weight`. It -#' also carries the `@variant` and `@bandwidths` properties. +#' Its `@results` tibble has one row per observation and intervention. The data +#' variant carries columns `.id` (the observation row index), `intervention` +#' (the intervention label), `value` (the intervened exposure, the constant +#' \eqn{a^*} for a single value and \eqn{d(a_i)} for a function), and `edp`. +#' The estimator variant replaces `edp` with `edp_outcome`, `edp_treatment`, +#' and `ideal_weight`. It also carries the `@variant` and `@bandwidths` +#' properties. #' #' [generics::glance()] returns a one-row tibble with `n` (the sample size), -#' `variant`, `n_values` (the number of intervention values), and the observed +#' `variant`, `n_values` (the number of interventions), and the observed #' range of every measure the variant reports. That is `edp_min` and `edp_max` #' for the data variant, and `edp_outcome_min`, `edp_outcome_max`, #' `edp_treatment_min`, `edp_treatment_max`, `ideal_weight_min`, and @@ -222,9 +231,7 @@ check_edp <- function( fn = "check_edp" ) - if (!is.null(values)) { - validate_values(values, exposure_type) - } + interventions <- resolve_interventions(values, exposure_vec, exposure_type) compute_names <- if (variant == "estimator") { unique(c(outcome_names, treatment_names)) @@ -244,14 +251,13 @@ check_edp <- function( } else { NA_real_ } - values_grid <- resolve_values(values, exposure_vec, exposure_type) exposure_kernels <- lapply( - values_grid, - function(a) { - exposure_kernel( + interventions, + function(intervention) { + exposure_kernel_target( exposure_vec, - a, + intervention$target, exposure_type, bw_exposure_resolved, categorical_similarity @@ -259,8 +265,13 @@ check_edp <- function( } ) - ids <- rep(seq_len(n), times = length(values_grid)) - value_column <- rep(values_grid, each = n) + labels <- vapply(interventions, function(x) x$label, character(1)) + ids <- rep(seq_len(n), times = length(interventions)) + intervention_column <- rep(labels, each = n) + value_column <- unlist( + lapply(interventions, function(x) rep_len(x$target, n)), + use.names = FALSE + ) if (variant == "data") { cov_bw <- resolve_covariate_bandwidths( @@ -275,10 +286,15 @@ check_edp <- function( categorical_similarity ) edp <- unlist( - lapply(exposure_kernels, function(e) as.double(similarity %*% e)), + lapply(exposure_kernels, function(e) edp_kernel_sum(similarity, e)), use.names = FALSE ) - results <- tibble::tibble(.id = ids, value = value_column, edp = edp) + results <- tibble::tibble( + .id = ids, + intervention = intervention_column, + value = value_column, + edp = edp + ) bandwidths <- list(exposure = bw_exposure_resolved, covariates = cov_bw) } else { outcome_bw <- resolve_covariate_bandwidths( @@ -304,15 +320,19 @@ check_edp <- function( categorical_similarity ) edp_outcome <- unlist( - lapply(exposure_kernels, function(e) as.double(outcome_similarity %*% e)), + lapply( + exposure_kernels, + function(e) edp_kernel_sum(outcome_similarity, e) + ), use.names = FALSE ) edp_treatment <- rep( rowSums(treatment_similarity), - times = length(values_grid) + times = length(interventions) ) results <- tibble::tibble( .id = ids, + intervention = intervention_column, value = value_column, edp_outcome = edp_outcome, edp_treatment = edp_treatment, @@ -333,7 +353,10 @@ check_edp <- function( params = list( variant = variant, kernel = kernel, - values = values_grid, + values = stats::setNames( + lapply(interventions, function(x) x$spec), + labels + ), categorical_similarity = categorical_similarity, bw_exposure = bw_exposure_resolved, bw_covariates = bw_covariates, @@ -375,6 +398,252 @@ resolve_covariate_names <- function( names(selection) } +#' Resolve every spelling of `values` into labeled intervention records +#' +#' Each record is a `list(label, spec, target)`: `label` names the intervention +#' in the results and on a plot, `spec` is the element as the user wrote it, and +#' `target` is the intervened exposure, either one value for the whole sample or +#' one value per observation. +#' +#' A vector, and the default grid, give one record per value labeled by the +#' value itself. A list may hold single values under the same rule and, for a +#' continuous exposure, functions of the observed exposure. Every label is +#' derived before any function runs, so a list that cannot be labeled fails +#' without evaluating user code. +#' +#' @param values The user-supplied values, or `NULL` for the default grid. +#' @param exposure_vec The observed exposure vector. +#' @param exposure_type The resolved exposure type. +#' @param arg_name The argument name used in error messages. +#' @param call The calling environment, used to build the error's call. +#' +#' @return A list of intervention records. +#' @keywords internal +#' @noRd +resolve_interventions <- function( + values, + exposure_vec, + exposure_type, + arg_name = "values", + call = rlang::caller_env() +) { + if (!is.null(values) && length(values) == 0) { + abort( + "{.arg {arg_name}} must contain at least one value.", + error_class = "positively_empty_error", + call = call + ) + } + if (!is.list(values)) { + if (!is.null(values)) { + validate_values(values, exposure_type, arg_name = arg_name, call = call) + } + grid <- resolve_values(values, exposure_vec, exposure_type) + return(lapply(seq_along(grid), function(i) { + value <- grid[[i]] + list(label = as.character(value), spec = value, target = value) + })) + } + + supplied <- names(values) %||% rep("", length(values)) + supplied[is.na(supplied)] <- "" + labels <- intervention_labels(values, supplied, arg_name, call) + validate_distinct_labels(labels, supplied, arg_name, call) + lapply(seq_along(values), function(i) { + resolve_one_intervention( + values[[i]], + labels[[i]], + exposure_vec, + exposure_type, + arg_name, + call + ) + }) +} + +#' Derive one label per element of a list of interventions +#' +#' A name is the label wherever the user wrote one. An unnamed single value is +#' labeled by the value itself, and an unnamed function has nothing to be +#' called by, so it is rejected here rather than reaching a plot without an +#' axis label. +#' +#' @param values A list of interventions. +#' @param supplied The names the user wrote, one per element and `""` where the +#' element has none. +#' @param arg_name The argument name used in error messages. +#' @param call The calling environment, used to build the error's call. +#' +#' @return A character vector of labels, one per element. +#' @keywords internal +#' @noRd +intervention_labels <- function(values, supplied, arg_name, call) { + vapply( + seq_along(values), + function(i) { + spec <- values[[i]] + if (!is.function(spec) && !(is.atomic(spec) && length(spec) == 1)) { + spec_class <- class(spec)[1] + abort( + c( + "{.arg {arg_name}} must contain single values or functions of the exposure.", + x = "Element {i} is a {.cls {spec_class}} of length {length(spec)}." + ), + error_class = "positively_type_error", + call = call + ) + } + if (nzchar(supplied[[i]])) { + return(supplied[[i]]) + } + if (is.function(spec)) { + abort( + c( + "Every function in {.arg {arg_name}} must be named.", + i = "The name labels the intervention in the results and on a plot.", + x = "Element {i} is an unnamed function." + ), + error_class = "positively_unnamed_error", + call = call + ) + } + as.character(spec) + }, + character(1) + ) +} + +#' Validate that no two interventions share a label +#' +#' @param labels The derived labels. +#' @param supplied The names the user wrote, one per element and `""` where the +#' element has none. +#' @param arg_name The argument name used in error messages. +#' @param call The calling environment, used to build the error's call. +#' +#' @return `labels`, invisibly, when every label is distinct. +#' @keywords internal +#' @noRd +validate_distinct_labels <- function(labels, supplied, arg_name, call) { + repeated <- unique(labels[duplicated(labels)]) + if (length(repeated) == 0) { + return(invisible(labels)) + } + bullets <- c( + "{.arg {arg_name}} must label every intervention distinctly.", + x = "Duplicated label{?s}: {.val {repeated}}." + ) + if (!all(nzchar(supplied[labels %in% repeated]))) { + bullets <- c( + bullets, + i = "An element without a name is labeled by its own value." + ) + } + abort(bullets, error_class = "positively_duplicate_error", call = call) +} + +#' Resolve one element of a list of interventions into a record +#' +#' @param spec The element as the user supplied it. +#' @param label The intervention's label. +#' @param exposure_vec The observed exposure vector. +#' @param exposure_type The resolved exposure type. +#' @param arg_name The argument name used in error messages. +#' @param call The calling environment, used to build the error's call. +#' +#' @return An intervention record. +#' @keywords internal +#' @noRd +resolve_one_intervention <- function( + spec, + label, + exposure_vec, + exposure_type, + arg_name, + call +) { + if (!is.function(spec)) { + validate_values(spec, exposure_type, arg_name = arg_name, call = call) + return(list(label = label, spec = spec, target = spec)) + } + if (exposure_type != "continuous") { + abort( + c( + "A function in {.arg {arg_name}} needs a continuous exposure.", + x = "The exposure is {exposure_type}.", + i = "A function maps each observed exposure to its intervened value." + ), + error_class = "positively_type_error", + call = call + ) + } + target <- spec(exposure_vec) + list( + label = label, + spec = spec, + target = validate_intervention_target( + target, + label, + length(exposure_vec), + arg_name, + call + ) + ) +} + +#' Validate what an intervention function returned +#' +#' The result stands in for the exposure, so it has to meet what the exposure +#' column itself must meet: numeric, complete, and finite. A single value is +#' recycled over the sample, the one output length that does not already say +#' which observation each target belongs to. +#' +#' @param target The value the function returned. +#' @param label The intervention's label. +#' @param n The number of observations. +#' @param arg_name The argument name used in error messages. +#' @param call The calling environment, used to build the error's call. +#' +#' @return A double vector of length `n`. +#' @keywords internal +#' @noRd +validate_intervention_target <- function(target, label, n, arg_name, call) { + if (!is.numeric(target)) { + target_class <- class(target)[1] + abort( + "The {.arg {arg_name}} function {.val {label}} must return a numeric vector, not a {.cls {target_class}}.", + error_class = "positively_type_error", + call = call + ) + } + if (!length(target) %in% c(1L, n)) { + abort( + c( + "The {.arg {arg_name}} function {.val {label}} must return one value per observation or a single value.", + x = "It returned {length(target)} value{?s} for {n} observation{?s}." + ), + error_class = "positively_size_error", + call = call + ) + } + target <- rep_len(as.double(target), n) + if (anyNA(target)) { + abort( + "The {.arg {arg_name}} function {.val {label}} must not return missing values.", + error_class = "positively_missing_error", + call = call + ) + } + if (!all(is.finite(target))) { + abort( + "The {.arg {arg_name}} function {.val {label}} must not return non-finite values.", + error_class = "positively_range_error", + call = call + ) + } + target +} + #' Resolve the intervention-value grid #' #' @param values The user-supplied values, or `NULL` for the default grid. @@ -627,20 +896,76 @@ exposure_kernel <- function( } } +#' The exposure-dimension kernel for one intervention target +#' +#' A single target is one point every observation is moved to, so the kernel is +#' one value per observation. A target per observation makes the kernel an +#' `n` by `n` matrix instead, rows indexing the intervened-on point and columns +#' the observed sample. +#' +#' @param exposure_vec The observed exposure vector. +#' @param target A single intervention value, or one per observation. +#' @param exposure_type The resolved exposure type. +#' @param bw The continuous-exposure half-distance (ignored when discrete). +#' @param categorical_similarity The non-matching kernel value for discrete +#' exposures. +#' +#' @return A numeric vector of length `n`, or an `n` by `n` numeric matrix. +#' @keywords internal +#' @noRd +exposure_kernel_target <- function( + exposure_vec, + target, + exposure_type, + bw, + categorical_similarity +) { + if (length(target) == 1) { + return(exposure_kernel( + exposure_vec, + target, + exposure_type, + bw, + categorical_similarity + )) + } + continuous_kernel(outer(target, as.double(exposure_vec), "-"), bw) +} + +#' Sum the product kernel over the sample for every observation +#' +#' @param similarity The `n` by `n` covariate similarity matrix. +#' @param kernel The exposure kernel: one value per observation for a single +#' target, or an `n` by `n` matrix for one target per observation. +#' +#' @return A numeric vector of length `n`. +#' @keywords internal +#' @noRd +edp_kernel_sum <- function(similarity, kernel) { + if (is.matrix(kernel)) { + return(rowSums(similarity * kernel)) + } + as.double(similarity %*% kernel) +} + #' The reparameterized Gaussian kernel over a difference vector #' #' Returns `0.5 ^ ((delta / bw) ^ 2)`, which equals `1` at `delta == 0` and #' `0.5` at `abs(delta) == bw`. A half-distance of `0` counts exact matches. #' -#' @param delta A numeric vector of differences. +#' @param delta A numeric vector or matrix of differences. #' @param bw The half-distance. #' -#' @return A numeric vector in `[0, 1]`. +#' @return A numeric object in `[0, 1]`, shaped like `delta`. #' @keywords internal #' @noRd continuous_kernel <- function(delta, bw) { if (bw == 0) { - return(as.double(delta == 0)) + matches <- delta == 0 + # Assigning the storage mode keeps a matrix a matrix, which as.double() + # would flatten. + storage.mode(matches) <- "double" + return(matches) } 0.5^((delta / bw)^2) } diff --git a/R/plot-edp.R b/R/plot-edp.R index f07e3fc..7cf4091 100644 --- a/R/plot-edp.R +++ b/R/plot-edp.R @@ -3,8 +3,8 @@ #' Plot an effective-data-points diagnostic #' #' Draws one view of a [check_edp()] result. `type = "histogram"` shows the -#' distribution of EDP faceted by intervention value, and `type = "ecdf"` shows -#' its empirical cumulative distribution colored by intervention value. For the +#' distribution of EDP faceted by intervention, and `type = "ecdf"` shows its +#' empirical cumulative distribution colored by intervention. For the #' estimator variant, `type = "scatter"` plots `edp_outcome` against #' `edp_treatment` colored by `ideal_weight`; it aborts for the data variant. #' @@ -44,6 +44,25 @@ method(autoplot, edp_result) <- function( } } +#' The results tibble keyed on the intervention label +#' +#' Every view that separates the interventions keys on the label rather than on +#' `value`, which a function intervention varies from observation to +#' observation. The levels follow the order the interventions were given in. +#' +#' @param results An `edp_result` results tibble. +#' +#' @return `results` with `intervention` as a factor. +#' @keywords internal +#' @noRd +edp_plot_data <- function(results) { + results$intervention <- factor( + results$intervention, + levels = unique(results$intervention) + ) + results +} + #' The histogram view of an EDP diagnostic #' #' @param object An `edp_result`. @@ -53,18 +72,17 @@ method(autoplot, edp_result) <- function( #' @noRd autoplot_edp_histogram <- function(object) { measure <- edp_primary_column(object@results) - plot_data <- object@results - plot_data$value <- factor(plot_data$value) + plot_data <- edp_plot_data(object@results) ggplot2::ggplot(plot_data, ggplot2::aes(x = .data[[measure]])) + ggplot2::geom_histogram(bins = 30, fill = "grey70", color = "white") + ggplot2::facet_wrap( - ggplot2::vars(.data$value), + ggplot2::vars(.data$intervention), labeller = ggplot2::label_both ) + ggplot2::labs( x = "Effective data points", y = "Observations", - title = "Effective data points by intervention value" + title = "Effective data points by intervention" ) } @@ -77,18 +95,17 @@ autoplot_edp_histogram <- function(object) { #' @noRd autoplot_edp_ecdf <- function(object) { measure <- edp_primary_column(object@results) - plot_data <- object@results - plot_data$value <- factor(plot_data$value) + plot_data <- edp_plot_data(object@results) ggplot2::ggplot( plot_data, - ggplot2::aes(x = .data[[measure]], color = .data$value) + ggplot2::aes(x = .data[[measure]], color = .data$intervention) ) + ggplot2::stat_ecdf() + ggplot2::labs( x = "Effective data points", y = "Cumulative proportion", - color = "Intervention value", - title = "Effective data points by intervention value" + color = "Intervention", + title = "Effective data points by intervention" ) } diff --git a/man/autoplot.edp_result.Rd b/man/autoplot.edp_result.Rd index f8adb4b..5516133 100644 --- a/man/autoplot.edp_result.Rd +++ b/man/autoplot.edp_result.Rd @@ -16,8 +16,8 @@ A \link[ggplot2:ggplot]{ggplot2::ggplot} object. } \description{ Draws one view of a \code{\link[=check_edp]{check_edp()}} result. \code{type = "histogram"} shows the -distribution of EDP faceted by intervention value, and \code{type = "ecdf"} shows -its empirical cumulative distribution colored by intervention value. For the +distribution of EDP faceted by intervention, and \code{type = "ecdf"} shows its +empirical cumulative distribution colored by intervention. For the estimator variant, \code{type = "scatter"} plots \code{edp_outcome} against \code{edp_treatment} colored by \code{ideal_weight}; it aborts for the data variant. } diff --git a/man/check_edp.Rd b/man/check_edp.Rd index 94d6654..8b69aa1 100644 --- a/man/check_edp.Rd +++ b/man/check_edp.Rd @@ -32,10 +32,17 @@ estimator variant, selected with tidyselect. Defaults to \code{.covariates}.} \item{.treatment_covariates}{The treatment-model covariate columns for the estimator variant, selected with tidyselect. Defaults to \code{.covariates}.} -\item{values}{The intervention values \eqn{a^*}. \code{NULL} (the default) uses -every level of a factor exposure and every observed value of a binary or -character exposure, and the deciles of the observed exposure for continuous -exposures.} +\item{values}{The interventions. \code{NULL} (the default) uses every level of a +factor exposure and every observed value of a binary or character exposure, +and the deciles of the observed exposure for continuous exposures. A vector +holds one intervention value \eqn{a^*} per element. A list holds the same +single values and, for a continuous exposure, functions of the observed +exposure such as \code{list(3, 5, lower = function(.x) .x - 5)}. A function is +called once on the whole exposure vector and returns one intervened +exposure \eqn{d(a_i)} per observation, or a single value used for every +observation. Each function needs a name, which labels the intervention; an +unnamed value is labeled by the value itself, and every label in a list +must be distinct.} \item{variant}{One of \code{"data"} (report a single \code{edp} per point, the default) or \code{"estimator"} (report \code{edp_outcome}, \code{edp_treatment}, and @@ -63,14 +70,16 @@ nothing of the column.} } \value{ An \code{edp_result} object, an S7 subclass of \link{positivity_diagnostic}. -Its \verb{@results} tibble has one row per observation and intervention value. -The data variant carries columns \code{.id} (the observation row index), \code{value} -(the intervention value \eqn{a^*}), and \code{edp}. The estimator variant -replaces \code{edp} with \code{edp_outcome}, \code{edp_treatment}, and \code{ideal_weight}. It -also carries the \verb{@variant} and \verb{@bandwidths} properties. +Its \verb{@results} tibble has one row per observation and intervention. The data +variant carries columns \code{.id} (the observation row index), \code{intervention} +(the intervention label), \code{value} (the intervened exposure, the constant +\eqn{a^*} for a single value and \eqn{d(a_i)} for a function), and \code{edp}. +The estimator variant replaces \code{edp} with \code{edp_outcome}, \code{edp_treatment}, +and \code{ideal_weight}. It also carries the \verb{@variant} and \verb{@bandwidths} +properties. \code{\link[generics:glance]{generics::glance()}} returns a one-row tibble with \code{n} (the sample size), -\code{variant}, \code{n_values} (the number of intervention values), and the observed +\code{variant}, \code{n_values} (the number of interventions), and the observed range of every measure the variant reports. That is \code{edp_min} and \code{edp_max} for the data variant, and \code{edp_outcome_min}, \code{edp_outcome_max}, \code{edp_treatment_min}, \code{edp_treatment_max}, \code{ideal_weight_min}, and diff --git a/tests/testthat/_snaps/check-edp.md b/tests/testthat/_snaps/check-edp.md index 0008329..081ae9c 100644 --- a/tests/testthat/_snaps/check-edp.md +++ b/tests/testthat/_snaps/check-edp.md @@ -184,6 +184,101 @@ Error in `check_edp()`: ! `.treatment_covariates` must select at least one column. +# the function-intervention error messages are stable + + Code + check_edp(data, exposure, x1, values = list(shift), exposure_type = "continuous") + Condition + Error in `check_edp()`: + ! Every function in `values` must be named. + i The name labels the intervention in the results and on a plot. + x Element 1 is an unnamed function. + +--- + + Code + check_edp(data, exposure, x1, values = stats::setNames(list(shift, shift), c( + "down", "down")), exposure_type = "continuous") + Condition + Error in `check_edp()`: + ! `values` must label every intervention distinctly. + x Duplicated label: "down". + +--- + + Code + check_edp(data, exposure, x1, values = stats::setNames(list(5, 1), c("", "5")), + exposure_type = "continuous") + Condition + Error in `check_edp()`: + ! `values` must label every intervention distinctly. + x Duplicated label: "5". + i An element without a name is labeled by its own value. + +--- + + Code + check_edp(data, exposure, x1, values = list(bad = function(.x) rep("a", length( + .x))), exposure_type = "continuous") + Condition + Error in `check_edp()`: + ! The `values` function "bad" must return a numeric vector, not a . + +--- + + Code + check_edp(data, exposure, x1, values = list(bad = function(.x) .x[1:2]), + exposure_type = "continuous") + Condition + Error in `check_edp()`: + ! The `values` function "bad" must return one value per observation or a single value. + x It returned 2 values for 60 observations. + +--- + + Code + check_edp(data, exposure, x1, values = list(bad = function(.x) replace(.x, 1, + NA_real_)), exposure_type = "continuous") + Condition + Error in `check_edp()`: + ! The `values` function "bad" must not return missing values. + +--- + + Code + check_edp(data, exposure, x1, values = list(bad = function(.x) replace(.x, 1, + Inf)), exposure_type = "continuous") + Condition + Error in `check_edp()`: + ! The `values` function "bad" must not return non-finite values. + +--- + + Code + check_edp(data, exposure, x1, values = list(bad = NULL), exposure_type = "continuous") + Condition + Error in `check_edp()`: + ! `values` must contain single values or functions of the exposure. + x Element 1 is a of length 0. + +--- + + Code + check_edp(data, exposure, x1, values = list(), exposure_type = "continuous") + Condition + Error in `check_edp()`: + ! `values` must contain at least one value. + +--- + + Code + check_edp(binary, exposure, c(x1, x2), values = list(down = shift)) + Condition + Error in `check_edp()`: + ! A function in `values` needs a continuous exposure. + x The exposure is binary. + i A function maps each observed exposure to its intervened value. + # the unused-bandwidth warnings are stable Code @@ -253,6 +348,19 @@ edp_treatment range: 8.995 to 98.272 ideal_weight range: 2.127 to 2.795 +# the print method is stable for a function intervention + + Code + print(res) + Output + + -- Effective data points ------------------------------------------------------- + Exposure: "exposure" (continuous) + Observations: 150 + Variant: data + Intervention values: 2 + edp range: 0.039 to 38.715 + # the EDP label and headline are stable Code diff --git a/tests/testthat/_snaps/check-edp/edp-ecdf-by-intervention-value.svg b/tests/testthat/_snaps/check-edp/edp-ecdf-by-intervention-value.svg index de592d1..77ae575 100644 --- a/tests/testthat/_snaps/check-edp/edp-ecdf-by-intervention-value.svg +++ b/tests/testthat/_snaps/check-edp/edp-ecdf-by-intervention-value.svg @@ -21,15 +21,15 @@ - - + + - - - - - + + + + + 0.00 @@ -42,22 +42,22 @@ - - - -10 -20 -30 -Effective data points + + + +10 +20 +30 +Effective data points Cumulative proportion - -Intervention value - - - - -0 -1 -Effective data points by intervention value + +Intervention + + + + +0 +1 +Effective data points by intervention diff --git a/tests/testthat/_snaps/check-edp/edp-histogram-by-intervention-value.svg b/tests/testthat/_snaps/check-edp/edp-histogram-by-intervention-value.svg index 3769268..edd7ce4 100644 --- a/tests/testthat/_snaps/check-edp/edp-histogram-by-intervention-value.svg +++ b/tests/testthat/_snaps/check-edp/edp-histogram-by-intervention-value.svg @@ -109,7 +109,7 @@ -value: 0 +intervention: 0 @@ -120,7 +120,7 @@ -value: 1 +intervention: 1 @@ -149,6 +149,6 @@ Effective data points Observations -Effective data points by intervention value +Effective data points by intervention diff --git a/tests/testthat/_snaps/check-positivity/positivity-check-panel.svg b/tests/testthat/_snaps/check-positivity/positivity-check-panel.svg index fa1f954..9e9b41a 100644 --- a/tests/testthat/_snaps/check-positivity/positivity-check-panel.svg +++ b/tests/testthat/_snaps/check-positivity/positivity-check-panel.svg @@ -174,7 +174,7 @@ -value: 0 +intervention: 0 @@ -185,7 +185,7 @@ -value: 1 +intervention: 1 @@ -210,7 +210,7 @@ Effective data points Observations -Effective data points by intervention value +Effective data points by intervention diff --git a/tests/testthat/test-check-edp.R b/tests/testthat/test-check-edp.R index a6be7c8..795c62f 100644 --- a/tests/testthat/test-check-edp.R +++ b/tests/testthat/test-check-edp.R @@ -9,8 +9,11 @@ # similarity kernel; only continuous exposures use the half-distance kernel. # Numeric covariates are continuous dimensions, factor and character covariates # are categorical dimensions. The fixed data-variant results columns are -# .id / value / edp; the estimator variant replaces edp with edp_outcome / -# edp_treatment / ideal_weight. EDP is meaningful relative across observations +# .id / intervention / value / edp; the estimator variant replaces edp with +# edp_outcome / edp_treatment / ideal_weight. An intervention is either a static +# value, where every row's `value` is that constant, or a function of the +# observed exposure, where each row's `value` is its own intervened exposure +# d(a_i). EDP is meaningful relative across observations # for a fixed covariate set, never against a universal threshold, so magnitude # claims are anchored to hand computations or to stratum-versus-stratum ratios. @@ -392,7 +395,7 @@ test_that("the data variant carries the fixed results columns and variant proper ) expect_s3_class(res@results, "tbl_df") - expect_setequal(names(res@results), c(".id", "value", "edp")) + expect_setequal(names(res@results), c(".id", "intervention", "value", "edp")) expect_type(res@results$edp, "double") expect_identical(nrow(res@results), 2L * 150L) expect_identical(res@variant, "data") @@ -453,7 +456,14 @@ test_that("the estimator variant swaps in the outcome and treatment columns", { expect_setequal( names(res@results), - c(".id", "value", "edp_outcome", "edp_treatment", "ideal_weight") + c( + ".id", + "intervention", + "value", + "edp_outcome", + "edp_treatment", + "ideal_weight" + ) ) expect_false("edp" %in% names(res@results)) expect_identical(res@variant, "estimator") @@ -936,6 +946,464 @@ test_that("ideal_weight is infinite where outcome support is zero", { expect_true(all(is.infinite(res@results$ideal_weight[subgroup]))) }) +# ---- Function-valued interventions ---------------------------------------- + +# Three observations five units apart on one constant categorical covariate, so +# the covariate similarity is 1 for every pair and EDP is the exposure kernel +# alone. At a half-distance of 5 every kernel weight is a negative power of two. +edp_shift_data <- function() { + tibble::tibble(exposure = c(0, 5, 10), s = c("a", "a", "a")) +} + +test_that("a shift intervention matches the hand-computed EDP", { + # d(a) = a - 5 moves the three targets to -5, 0, and 5. Against the observed + # exposures 0, 5, and 10 the weights 0.5 ^ ((delta / 5) ^ 2) sum to + # 0.5 + 0.5 ^ 4 + 0.5 ^ 9, then 1 + 0.5 + 0.5 ^ 4, then 0.5 + 1 + 0.5. + res <- check_edp( + edp_shift_data(), + exposure, + s, + bw_exposure = 5, + values = list(down = function(.x) .x - 5), + exposure_type = "continuous" + ) + + ordered <- res@results[order(res@results$.id), , drop = FALSE] + expect_equal(ordered$value, c(-5, 0, 5), tolerance = 1e-8) + expect_equal( + ordered$edp, + c(0.564453125, 1.5625, 2), + tolerance = 1e-8 + ) +}) + +test_that("a shift intervention matches the hand-computed two-dimension product", { + # covariate x1 with h = 1, exposure with h = 0.5, obs (x1, a) = + # (0, 0), (0, 1), (1, 0), shifted by d(a) = a + 0.5. Rows 1 and 3 are targeted + # at 0.5, the static value of the two-dimension test above, and reproduce its + # 1.25 and 1.0. Row 2 is targeted at 1.5, giving + # 0.5 ^ 9 + 0.5 + 0.5 * 0.5 ^ 9 = 0.5029296875. + data <- tibble::tibble(exposure = c(0, 1, 0), x1 = c(0, 0, 1)) + res <- check_edp( + data, + exposure, + x1, + bw_exposure = 0.5, + bw_covariates = 1, + values = list(up = function(.x) .x + 0.5), + exposure_type = "continuous" + ) + + ordered <- res@results[order(res@results$.id), , drop = FALSE] + expect_equal(ordered$value, c(0.5, 1.5, 0.5), tolerance = 1e-8) + expect_equal( + ordered$edp, + c(1.25, 0.5029296875, 1), + tolerance = 1e-8 + ) +}) + +test_that("both variants carry .id, intervention, value, then the measures", { + data <- sim_edp_gaussian(40) + n <- nrow(data) + down <- function(.x) .x - 1 + + res <- check_edp( + data, + exposure, + x1, + values = list(down = down, 0), + exposure_type = "continuous" + ) + expect_identical(names(res@results), c(".id", "intervention", "value", "edp")) + expect_type(res@results$intervention, "character") + expect_identical(res@results$intervention, rep(c("down", "0"), each = n)) + expect_identical(res@results$.id, rep(seq_len(n), times = 2)) + # A function intervention gives every observation its own target; a static one + # sends them all to the same place. + expect_identical(res@results$value[seq_len(n)], data$exposure - 1) + expect_identical(res@results$value[n + seq_len(n)], rep(0, n)) + + estimator <- check_edp( + data, + exposure, + x1, + values = list(down = down, 0), + variant = "estimator", + exposure_type = "continuous" + ) + expect_identical( + names(estimator@results), + c( + ".id", + "intervention", + "value", + "edp_outcome", + "edp_treatment", + "ideal_weight" + ) + ) + expect_type(estimator@results$intervention, "character") + expect_identical( + estimator@results$intervention, + rep(c("down", "0"), each = n) + ) +}) + +test_that("the static path keeps its exact numbers under the new schema", { + # Every kernel factor here is a negative power of two, so the sums are exact + # in double precision and the released numbers pin bit for bit. + res <- check_edp( + edp_hand_data(), + exposure, + x1, + bw_exposure = 1, + bw_covariates = 1, + values = c(0, 1), + exposure_type = "continuous" + ) + + expect_identical( + res@results, + tibble::tibble( + .id = rep(1:3, times = 2), + intervention = rep(c("0", "1"), each = 3), + value = rep(c(0, 1), each = 3), + edp = rep(c(1.5625, 2), each = 3) + ) + ) +}) + +test_that("a list of scalars matches the equivalent vector", { + data <- sim_edp_gaussian(60) + from_list <- check_edp( + data, + exposure, + x1, + values = list(1, 5), + exposure_type = "continuous" + ) + from_vector <- check_edp( + data, + exposure, + x1, + values = c(1, 5), + exposure_type = "continuous" + ) + + expect_identical(from_list@results, from_vector@results) + expect_identical(unique(from_list@results$intervention), c("1", "5")) + # Both spellings of the same grid have to leave the same record behind, or + # glance() and print() would report on the spelling rather than the grid. + expect_identical( + from_vector@params$values, + stats::setNames(list(1, 5), c("1", "5")) + ) + expect_identical(from_list@params$values, from_vector@params$values) +}) + +test_that("a named scalar takes its name as the intervention label", { + data <- sim_edp_gaussian(60) + n <- nrow(data) + res <- check_edp( + data, + exposure, + x1, + values = list(low = 1), + exposure_type = "continuous" + ) + + # A name labels a scalar as readily as it labels a function, and the label is + # all that changes: the target is still the number given. + expect_identical(res@results$intervention, rep("low", n)) + expect_identical(res@results$value, rep(1, n)) + expect_identical(names(res@params$values), "low") +}) + +test_that("a function returning a constant matches the equivalent static value", { + data <- sim_edp_gaussian(60) + from_function <- check_edp( + data, + exposure, + x1, + values = list(static = function(.x) rep(3, length(.x))), + exposure_type = "continuous" + ) + from_static <- check_edp( + data, + exposure, + x1, + values = 3, + exposure_type = "continuous" + ) + + expect_identical(from_function@results$.id, from_static@results$.id) + expect_equal(from_function@results$value, from_static@results$value) + expect_equal( + from_function@results$edp, + from_static@results$edp, + tolerance = 1e-12 + ) + expect_identical( + from_function@results$intervention, + rep("static", nrow(data)) + ) +}) + +test_that("a function returning one value is recycled over the observations", { + # Length-1 output is the one output length that does not already say which + # observation each target belongs to, so it has to be recycled rather than + # rejected for being shorter than the sample. + data <- sim_edp_gaussian(60) + n <- nrow(data) + from_function <- check_edp( + data, + exposure, + x1, + values = list(flat = function(.x) 3), + exposure_type = "continuous" + ) + from_static <- check_edp( + data, + exposure, + x1, + values = 3, + exposure_type = "continuous" + ) + + expect_identical(from_function@results$intervention, rep("flat", n)) + expect_identical(from_function@results$value, rep(3, n)) + expect_identical(from_function@results$value, from_static@results$value) + expect_equal( + from_function@results$edp, + from_static@results$edp, + tolerance = 1e-12 + ) +}) + +test_that("an intervention function is called once on the whole exposure", { + data <- sim_edp_gaussian(40) + calls <- 0L + received <- NULL + shift <- function(.x) { + calls <<- calls + 1L + received <<- .x + .x - 5 + } + + res <- check_edp( + data, + exposure, + x1, + values = list(down = shift), + exposure_type = "continuous" + ) + + # Checking the output and computing with it must share one invocation: + # applying the function a second time would double any side effect the user + # wrote and any cost they paid for. + expect_identical(calls, 1L) + # The argument is the observed exposure in full, not one observation at a time + # and not a subset. + expect_identical(received, data$exposure) + expect_identical(nrow(res@results), nrow(data)) +}) + +test_that("a function intervention leaves edp_treatment alone", { + # Exact matching (bw_exposure = 0) on four evenly spaced exposures sharing one + # constant covariate: d(a) = a + 1 places rows 1 to 3 on an observed exposure + # and row 4 one unit past the largest, where the outcome model has no support. + data <- tibble::tibble(exposure = c(0, 1, 2, 3), s = rep("a", 4)) + res <- check_edp( + data, + exposure, + s, + variant = "estimator", + bw_exposure = 0, + values = list(up = function(.x) .x + 1, 2), + exposure_type = "continuous" + ) + + expect_identical(res@results$intervention, rep(c("up", "2"), each = 4)) + # The treatment model carries no exposure dimension, so no intervention moves + # it. + expect_length(unique(res@results$edp_treatment), 1) + expect_equal(unique(res@results$edp_treatment), 4) + expect_equal(res@results$edp_outcome, c(1, 1, 1, 0, 1, 1, 1, 1)) + expect_equal(res@results$ideal_weight, c(4, 4, 4, Inf, 4, 4, 4, 4)) +}) + +test_that("glance() counts a list of interventions and keeps the normalized list", { + data <- sim_edp_gaussian(60) + down <- function(.x) .x - 1 + res <- check_edp( + data, + exposure, + x1, + values = list(down = down, 0), + exposure_type = "continuous" + ) + + expect_identical(generics::glance(res)$n_values, 2L) + expect_type(res@params$values, "list") + expect_identical(names(res@params$values), c("down", "0")) + expect_identical(res@params$values$down, down) + expect_identical(res@params$values[["0"]], 0) +}) + +test_that("check_edp() rejects a function intervention it cannot label or apply", { + local_quiet() + data <- sim_edp_gaussian(60) + shift <- function(.x) .x - 1 + + # A label becomes axis text, so an unnamed function leaves nothing to draw. + expect_error( + check_edp( + data, + exposure, + x1, + values = list(shift), + exposure_type = "continuous" + ), + class = "positively_unnamed_error" + ) + expect_error( + check_edp( + data, + exposure, + x1, + values = stats::setNames(list(shift, shift), c("down", "down")), + exposure_type = "continuous" + ), + class = "positively_duplicate_error" + ) + expect_error( + check_edp( + data, + exposure, + x1, + values = stats::setNames(list(1, 2), c("a", "a")), + exposure_type = "continuous" + ), + class = "positively_duplicate_error" + ) + + # A function maps observed exposure to intervened exposure, which only a + # continuous exposure has. + binary <- dgp_good_positivity(n = 60, seed = 1) + expect_error( + check_edp(binary, exposure, c(x1, x2), values = list(down = shift)), + class = "positively_type_error" + ) + categorical <- sim_edp_categorical(60) + expect_error( + check_edp( + categorical, + exposure, + z2, + values = list(down = shift), + exposure_type = "categorical" + ), + class = "positively_type_error" + ) +}) + +test_that("check_edp() rejects a function whose output cannot be an exposure", { + data <- sim_edp_gaussian(60) + + expect_error( + check_edp( + data, + exposure, + x1, + values = list(bad = function(.x) rep("a", length(.x))), + exposure_type = "continuous" + ), + class = "positively_type_error" + ) + expect_error( + check_edp( + data, + exposure, + x1, + values = list(bad = function(.x) .x[1:2]), + exposure_type = "continuous" + ), + class = "positively_size_error" + ) + expect_error( + check_edp( + data, + exposure, + x1, + values = list(bad = function(.x) replace(.x, 1, NA_real_)), + exposure_type = "continuous" + ), + class = "positively_missing_error" + ) + expect_error( + check_edp( + data, + exposure, + x1, + values = list(bad = function(.x) replace(.x, 1, Inf)), + exposure_type = "continuous" + ), + class = "positively_range_error" + ) +}) + +test_that("check_edp() rejects a list element that is neither a scalar nor a function", { + data <- sim_edp_gaussian(60) + + expect_error( + check_edp( + data, + exposure, + x1, + values = list(bad = NULL), + exposure_type = "continuous" + ), + class = "positively_type_error" + ) + expect_error( + check_edp( + data, + exposure, + x1, + values = list(bad = c(1, 2)), + exposure_type = "continuous" + ), + class = "positively_type_error" + ) + expect_error( + check_edp( + data, + exposure, + x1, + values = list(), + exposure_type = "continuous" + ), + class = "positively_empty_error" + ) +}) + +test_that("an error raised inside a user function propagates unchanged", { + data <- sim_edp_gaussian(60) + # The user's own failure is the one worth reading, so nothing may relabel it. + err <- expect_error( + check_edp( + data, + exposure, + x1, + values = list(boom = function(.x) stop("kaboom")), + exposure_type = "continuous" + ), + "kaboom" + ) + expect_false(inherits(err, "positively_error")) +}) + # ---- Autoplot contract ---------------------------------------------------- test_that("autoplot() returns a ggplot for each data-variant type", { @@ -1094,6 +1562,45 @@ test_that("the scatter view builds without a colour label when every weight is i expect_false(is.null(plot$labels$subtitle)) }) +test_that("the histogram view facets on the intervention label", { + # Under a function intervention `value` differs for nearly every observation, + # so a facet keyed on it would draw one panel per row. + data <- sim_edp_gaussian(60) + res <- check_edp( + data, + exposure, + x1, + values = list(down = function(.x) .x - 1, 0), + exposure_type = "continuous" + ) + + layout <- ggplot2::ggplot_build( + ggplot2::autoplot(res, type = "histogram") + )$layout$layout + + expect_true("intervention" %in% names(layout)) + expect_false("value" %in% names(layout)) + expect_identical(nrow(layout), 2L) + expect_setequal(as.character(layout$intervention), c("down", "0")) +}) + +test_that("the ecdf view colors by the intervention label", { + data <- sim_edp_gaussian(60) + res <- check_edp( + data, + exposure, + x1, + values = list(down = function(.x) .x - 1, 0), + exposure_type = "continuous" + ) + + built <- ggplot2::ggplot_build(ggplot2::autoplot(res, type = "ecdf")) + + # One curve per intervention, not one per distinct intervened exposure. + expect_length(unique(built$data[[1]]$group), 2) + expect_length(unique(built$data[[1]]$colour), 2) +}) + test_that("EDP autoplot views render as expected", { local_quiet() announce_doppelganger( @@ -1451,6 +1958,111 @@ test_that("the value and estimator-selection error messages are stable", { )) }) +test_that("the function-intervention error messages are stable", { + local_quiet() + data <- sim_edp_gaussian(60) + shift <- function(.x) .x - 1 + + expect_snapshot_abort( + check_edp( + data, + exposure, + x1, + values = list(shift), + exposure_type = "continuous" + ), + class = "positively_unnamed_error" + ) + expect_snapshot_abort( + check_edp( + data, + exposure, + x1, + values = stats::setNames(list(shift, shift), c("down", "down")), + exposure_type = "continuous" + ), + class = "positively_duplicate_error" + ) + # An unnamed scalar is labelled by its own value, so it can collide with a + # name the user wrote out even though the two elements look nothing alike. + expect_snapshot_abort( + check_edp( + data, + exposure, + x1, + values = stats::setNames(list(5, 1), c("", "5")), + exposure_type = "continuous" + ), + class = "positively_duplicate_error" + ) + expect_snapshot_abort( + check_edp( + data, + exposure, + x1, + values = list(bad = function(.x) rep("a", length(.x))), + exposure_type = "continuous" + ), + class = "positively_type_error" + ) + expect_snapshot_abort( + check_edp( + data, + exposure, + x1, + values = list(bad = function(.x) .x[1:2]), + exposure_type = "continuous" + ), + class = "positively_size_error" + ) + expect_snapshot_abort( + check_edp( + data, + exposure, + x1, + values = list(bad = function(.x) replace(.x, 1, NA_real_)), + exposure_type = "continuous" + ), + class = "positively_missing_error" + ) + expect_snapshot_abort( + check_edp( + data, + exposure, + x1, + values = list(bad = function(.x) replace(.x, 1, Inf)), + exposure_type = "continuous" + ), + class = "positively_range_error" + ) + expect_snapshot_abort( + check_edp( + data, + exposure, + x1, + values = list(bad = NULL), + exposure_type = "continuous" + ), + class = "positively_type_error" + ) + expect_snapshot_abort( + check_edp( + data, + exposure, + x1, + values = list(), + exposure_type = "continuous" + ), + class = "positively_empty_error" + ) + + binary <- dgp_good_positivity(n = 60, seed = 1) + expect_snapshot_abort( + check_edp(binary, exposure, c(x1, x2), values = list(down = shift)), + class = "positively_type_error" + ) +}) + test_that("the unused-bandwidth warnings are stable", { local_quiet() withr::local_options(warn = 0) @@ -1513,6 +2125,19 @@ test_that("the print method is stable", { expect_snapshot(print(estimator)) }) +test_that("the print method is stable for a function intervention", { + local_quiet() + + res <- check_edp( + sim_edp_gaussian(150), + exposure, + x1, + values = list(down = function(.x) .x - 1, 0), + exposure_type = "continuous" + ) + expect_snapshot(print(res)) +}) + # ---- Display methods ------------------------------------------------------- # Seven intervention values under the data variant, which is the variant whose From 4f5ea643613515f351e047e3bff7a66900458134 Mon Sep 17 00:00:00 2001 From: Malcolm Barrett Date: Fri, 21 Aug 2026 19:32:19 +0200 Subject: [PATCH 2/7] add a pull_plot_data() accessor returning what each EDP view draws --- NAMESPACE | 1 + R/plot-edp.R | 37 +-- R/pull-plot-data.R | 128 ++++++++++ _pkgdown.yml | 4 +- man/autoplot.edp_result.Rd | 1 + man/pull_plot_data.Rd | 57 +++++ tests/testthat/_snaps/pull-plot-data.md | 18 ++ tests/testthat/test-pull-plot-data.R | 318 ++++++++++++++++++++++++ 8 files changed, 531 insertions(+), 33 deletions(-) create mode 100644 R/pull-plot-data.R create mode 100644 man/pull_plot_data.Rd create mode 100644 tests/testthat/_snaps/pull-plot-data.md create mode 100644 tests/testthat/test-pull-plot-data.R diff --git a/NAMESPACE b/NAMESPACE index 694ca2a..a24611a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -21,6 +21,7 @@ export(hdr_density_normal) export(new_hdr_density) export(positivity_check) export(positivity_diagnostic) +export(pull_plot_data) export(sniff_violations) export(tidy) import(S7) diff --git a/R/plot-edp.R b/R/plot-edp.R index 7cf4091..c17bb4c 100644 --- a/R/plot-edp.R +++ b/R/plot-edp.R @@ -7,6 +7,7 @@ #' empirical cumulative distribution colored by intervention. For the #' estimator variant, `type = "scatter"` plots `edp_outcome` against #' `edp_treatment` colored by `ideal_weight`; it aborts for the data variant. +#' [pull_plot_data()] returns the tibble any of these views draws. #' #' @param object An `edp_result` from [check_edp()]. #' @param type One of `"histogram"`, `"ecdf"`, or `"scatter"`. @@ -44,25 +45,6 @@ method(autoplot, edp_result) <- function( } } -#' The results tibble keyed on the intervention label -#' -#' Every view that separates the interventions keys on the label rather than on -#' `value`, which a function intervention varies from observation to -#' observation. The levels follow the order the interventions were given in. -#' -#' @param results An `edp_result` results tibble. -#' -#' @return `results` with `intervention` as a factor. -#' @keywords internal -#' @noRd -edp_plot_data <- function(results) { - results$intervention <- factor( - results$intervention, - levels = unique(results$intervention) - ) - results -} - #' The histogram view of an EDP diagnostic #' #' @param object An `edp_result`. @@ -72,7 +54,7 @@ edp_plot_data <- function(results) { #' @noRd autoplot_edp_histogram <- function(object) { measure <- edp_primary_column(object@results) - plot_data <- edp_plot_data(object@results) + plot_data <- pull_plot_data(object, type = "histogram") ggplot2::ggplot(plot_data, ggplot2::aes(x = .data[[measure]])) + ggplot2::geom_histogram(bins = 30, fill = "grey70", color = "white") + ggplot2::facet_wrap( @@ -95,7 +77,7 @@ autoplot_edp_histogram <- function(object) { #' @noRd autoplot_edp_ecdf <- function(object) { measure <- edp_primary_column(object@results) - plot_data <- edp_plot_data(object@results) + plot_data <- pull_plot_data(object, type = "ecdf") ggplot2::ggplot( plot_data, ggplot2::aes(x = .data[[measure]], color = .data$intervention) @@ -118,17 +100,8 @@ autoplot_edp_ecdf <- function(object) { #' @keywords internal #' @noRd autoplot_edp_scatter <- function(object, call = rlang::caller_env()) { - if (object@variant != "estimator") { - abort( - c( - "The scatter view needs the estimator variant.", - i = "Rerun {.fn check_edp} with {.code variant = \"estimator\"}." - ), - error_class = "positively_variant_error", - call = call - ) - } - results <- object@results + edp_require_estimator(object, call = call) + results <- pull_plot_data(object, type = "scatter") finite_rows <- results[is.finite(results$ideal_weight), , drop = FALSE] infinite_rows <- results[is.infinite(results$ideal_weight), , drop = FALSE] diff --git a/R/pull-plot-data.R b/R/pull-plot-data.R new file mode 100644 index 0000000..4096c74 --- /dev/null +++ b/R/pull-plot-data.R @@ -0,0 +1,128 @@ +# S7 methods defined here are wired up at load time by S7::methods_register(). +# The file sorts after every class it dispatches on. + +#' Pull the data behind a plot +#' +#' `pull_plot_data()` returns the tibble one view of a diagnostic result draws. +#' A reader can rebuild the figure with [ggplot2::ggplot()] directly, take the +#' numbers somewhere else, or read what a view is showing without rendering it. +#' +#' @details +#' A view returns exactly one tibble, whatever the figure does with it. Where a +#' view splits its rows across layers, the whole frame comes back and the split +#' stays derivable from it. +#' +#' `type` names the view. Each method mirrors its class's [autoplot()] menu, +#' including the view that menu defaults to, so `pull_plot_data(x)` returns the +#' frame `autoplot(x)` draws. Every view is built from the frame this returns, +#' so the figure and the frame cannot disagree, and a view a result cannot draw +#' is refused here exactly as `autoplot()` refuses it. +#' +#' For a [check_edp()] result: +#' +#' * `"histogram"` and `"ecdf"` return the results at their own grain, one row +#' per observation and intervention, with `intervention` as a factor whose +#' levels follow the order the interventions were given in. The label is what +#' separates the interventions rather than `value`, which a function +#' intervention varies from observation to observation. +#' * `"scatter"` returns the whole estimator results tibble. The view draws the +#' rows of finite `ideal_weight` in one layer and the infinite rows in +#' another, and both are read off `ideal_weight`. It needs the estimator +#' variant, and asking for it after a data-variant run is an error. +#' +#' @param x A [positivity_diagnostic]. +#' @param ... Passed to methods. Every method takes `type`, the name of the view +#' whose data to return, matched against that class's [autoplot()] menu. +#' +#' @return A [tibble][tibble::tibble] holding one view's data. +#' +#' @examples +#' set.seed(1) +#' n <- 100 +#' x1 <- rnorm(n) +#' dose <- rnorm(n, mean = x1) +#' df <- data.frame(dose = dose, x1 = x1) +#' result <- check_edp(df, dose, x1, values = c(0, 1), exposure_type = "continuous") +#' +#' pull_plot_data(result) +#' pull_plot_data(result, type = "ecdf") +#' @export +pull_plot_data <- new_generic("pull_plot_data", "x", function(x, ...) { + S7::S7_dispatch() +}) + +# ---- Effective data points ------------------------------------------------- + +#' The results tibble keyed on the intervention label +#' +#' Every view that separates the interventions keys on the label rather than on +#' `value`, which a function intervention varies from observation to +#' observation. The levels come from the normalized `values` rather than from +#' the rows, so a facet strip and an axis take the order the run was given +#' instead of an order the results happen to carry. +#' +#' Two interventions can carry the same label, either from repeated values or +#' from values that round to one label, so the levels are deduplicated. Such +#' interventions share a level and a view draws their rows together; every row +#' of the results is kept either way. +#' +#' @param x An `edp_result`. +#' +#' @return The results with `intervention` as a factor whose levels are the +#' distinct labels in grid order. +#' @keywords internal +#' @noRd +edp_observation_frame <- function(x) { + results <- x@results + results$intervention <- factor( + results$intervention, + levels = unique(names(x@params$values)) + ) + results +} + +#' Require the estimator variant for the scatter view +#' +#' The scatter view reads `edp_outcome`, `edp_treatment`, and `ideal_weight`, +#' which only the estimator variant computes. Both the plot and the accessor +#' gate here, so a reader is told the same thing whichever they asked for and is +#' shown the call they made. +#' +#' @param x An `edp_result`. +#' @param call The call to blame, or the calling environment. +#' +#' @return `NULL`, invisibly, when the run is the estimator variant. +#' @keywords internal +#' @noRd +edp_require_estimator <- function(x, call = rlang::caller_env()) { + if (x@variant != "estimator") { + abort( + c( + "The scatter view needs the estimator variant.", + i = "Rerun {.fn check_edp} with {.code variant = \"estimator\"}." + ), + error_class = "positively_variant_error", + call = call + ) + } + invisible(NULL) +} + +method(pull_plot_data, edp_result) <- function( + x, + type = c("histogram", "ecdf", "scatter"), + ... +) { + # The reader's own call rather than the dispatched method, which names a + # function nobody typed. + type <- resolve_arg_match( + rlang::arg_match(type), + call = rlang::caller_call() + ) + if (type == "scatter") { + edp_require_estimator(x, call = rlang::caller_call()) + x@results + } else { + edp_observation_frame(x) + } +} diff --git a/_pkgdown.yml b/_pkgdown.yml index 59bea78..973c2c5 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -37,8 +37,10 @@ reference: - check_density_ratios - title: Plots desc: > - Autoplot methods that visualize each diagnostic result. + Autoplot methods that visualize each diagnostic result, and the accessor + that returns the data behind a view. contents: + - pull_plot_data - autoplot.positivity_check - autoplot.edp_result - autoplot.port_result diff --git a/man/autoplot.edp_result.Rd b/man/autoplot.edp_result.Rd index 5516133..1d0eee1 100644 --- a/man/autoplot.edp_result.Rd +++ b/man/autoplot.edp_result.Rd @@ -20,6 +20,7 @@ distribution of EDP faceted by intervention, and \code{type = "ecdf"} shows its empirical cumulative distribution colored by intervention. For the estimator variant, \code{type = "scatter"} plots \code{edp_outcome} against \code{edp_treatment} colored by \code{ideal_weight}; it aborts for the data variant. +\code{\link[=pull_plot_data]{pull_plot_data()}} returns the tibble any of these views draws. } \examples{ set.seed(1) diff --git a/man/pull_plot_data.Rd b/man/pull_plot_data.Rd new file mode 100644 index 0000000..c2f8819 --- /dev/null +++ b/man/pull_plot_data.Rd @@ -0,0 +1,57 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pull-plot-data.R +\name{pull_plot_data} +\alias{pull_plot_data} +\title{Pull the data behind a plot} +\usage{ +pull_plot_data(x, ...) +} +\arguments{ +\item{x}{A \link{positivity_diagnostic}.} + +\item{...}{Passed to methods. Every method takes \code{type}, the name of the view +whose data to return, matched against that class's \code{\link[=autoplot]{autoplot()}} menu.} +} +\value{ +A \link[tibble:tibble]{tibble} holding one view's data. +} +\description{ +\code{pull_plot_data()} returns the tibble one view of a diagnostic result draws. +A reader can rebuild the figure with \code{\link[ggplot2:ggplot]{ggplot2::ggplot()}} directly, take the +numbers somewhere else, or read what a view is showing without rendering it. +} +\details{ +A view returns exactly one tibble, whatever the figure does with it. Where a +view splits its rows across layers, the whole frame comes back and the split +stays derivable from it. + +\code{type} names the view. Each method mirrors its class's \code{\link[=autoplot]{autoplot()}} menu, +including the view that menu defaults to, so \code{pull_plot_data(x)} returns the +frame \code{autoplot(x)} draws. Every view is built from the frame this returns, +so the figure and the frame cannot disagree, and a view a result cannot draw +is refused here exactly as \code{autoplot()} refuses it. + +For a \code{\link[=check_edp]{check_edp()}} result: +\itemize{ +\item \code{"histogram"} and \code{"ecdf"} return the results at their own grain, one row +per observation and intervention, with \code{intervention} as a factor whose +levels follow the order the interventions were given in. The label is what +separates the interventions rather than \code{value}, which a function +intervention varies from observation to observation. +\item \code{"scatter"} returns the whole estimator results tibble. The view draws the +rows of finite \code{ideal_weight} in one layer and the infinite rows in +another, and both are read off \code{ideal_weight}. It needs the estimator +variant, and asking for it after a data-variant run is an error. +} +} +\examples{ +set.seed(1) +n <- 100 +x1 <- rnorm(n) +dose <- rnorm(n, mean = x1) +df <- data.frame(dose = dose, x1 = x1) +result <- check_edp(df, dose, x1, values = c(0, 1), exposure_type = "continuous") + +pull_plot_data(result) +pull_plot_data(result, type = "ecdf") +} diff --git a/tests/testthat/_snaps/pull-plot-data.md b/tests/testthat/_snaps/pull-plot-data.md new file mode 100644 index 0000000..bb57815 --- /dev/null +++ b/tests/testthat/_snaps/pull-plot-data.md @@ -0,0 +1,18 @@ +# pull_plot_data() and autoplot() refuse the scatter view alike + + Code + pull_plot_data(plain, type = "scatter") + Condition + Error in `pull_plot_data()`: + ! The scatter view needs the estimator variant. + i Rerun `check_edp()` with `variant = "estimator"`. + +--- + + Code + ggplot2::autoplot(plain, type = "scatter") + Condition + Error in `autoplot.positively::edp_result`: + ! The scatter view needs the estimator variant. + i Rerun `check_edp()` with `variant = "estimator"`. + diff --git a/tests/testthat/test-pull-plot-data.R b/tests/testthat/test-pull-plot-data.R new file mode 100644 index 0000000..dc4d9b5 --- /dev/null +++ b/tests/testthat/test-pull-plot-data.R @@ -0,0 +1,318 @@ +# `pull_plot_data()` hands back exactly the frame one autoplot view draws, so +# nearly every block here is a drift test: the accessor's frame against what +# ggplot2 reports the corresponding layers drawing. A frame that stopped +# matching its view would still look reasonable on its own, which is what these +# comparisons are for. + +# ---- Scenario generators -------------------------------------------------- + +# Gaussian exposure with an independent Gaussian covariate. Both intervention +# values sit inside the support, so every EDP is positive and finite. +sim_pull_gaussian <- function(n, seed = 1) { + withr::local_seed(seed) + tibble::tibble(exposure = stats::rnorm(n), x1 = stats::rnorm(n)) +} + +# Three strata under exact categorical matching: two are treated at different +# rates and get distinct finite ideal weights, and the third is never treated, +# so its outcome support at a* = 1 is zero and its ideal weight is infinite. The +# scatter view draws the finite and the infinite rows in separate layers. +sim_pull_mixed_support <- function() { + data.frame( + exposure = c( + rep(c(0L, 1L), c(4L, 4L)), + rep(c(0L, 1L), c(4L, 2L)), + rep(0L, 5L) + ), + s = factor(rep(c("a", "b", "c"), c(8L, 6L, 5L))) + ) +} + +# ---- View menus, read off the methods ------------------------------------- + +# The menus are read from the methods rather than restated, so a block cannot +# pass against a view the plot no longer offers. The S7 methods on `autoplot()` +# register as S3 methods on ggplot2's generic; `pull_plot_data()` is an S7 +# generic and its methods are reached through S7. +autoplot_types <- function(object) { + method <- utils::getS3method("autoplot", class(object)[[1]]) + eval(formals(method)$type) +} + +pull_plot_data_types <- function(class) { + eval(formals(S7::method(pull_plot_data, class))$type) +} + +deparsed_call <- function(condition) { + paste(deparse(conditionCall(condition)), collapse = " ") +} + +# ---- The generic ---------------------------------------------------------- + +test_that("pull_plot_data() is an exported generic with an edp_result method", { + expect_true("pull_plot_data" %in% getNamespaceExports("positively")) + expect_s3_class(pull_plot_data, "S7_generic") + expect_identical(names(formals(pull_plot_data))[[1]], "x") + + data <- sim_pull_gaussian(60) + res <- check_edp( + data, + exposure, + x1, + values = c(0, 1), + exposure_type = "continuous" + ) + + expect_true(is.function(S7::method(pull_plot_data, edp_result))) + expect_s3_class(pull_plot_data(res), "tbl_df") +}) + +test_that("pull_plot_data() defaults to the frame the default view draws", { + data <- sim_pull_gaussian(60) + res <- check_edp( + data, + exposure, + x1, + values = c(0, 1), + exposure_type = "continuous" + ) + + # Whichever view heads the menu is the default of both, so the two defaults + # are compared rather than named. + default_type <- autoplot_types(res)[[1]] + expect_identical( + pull_plot_data(res), + pull_plot_data(res, type = default_type) + ) + + # The default view draws a single frame at the plot level, so what the plot + # holds is what it draws. + expect_identical(pull_plot_data(res), ggplot2::autoplot(res)$data) +}) + +test_that("the type menu mirrors autoplot's and refuses an unknown view", { + data <- sim_pull_gaussian(60) + res <- check_edp( + data, + exposure, + x1, + values = c(0, 1), + exposure_type = "continuous" + ) + + # The accessor exists to hand back what a view draws, so it offers exactly the + # views autoplot() offers, in the same order. + expect_identical(pull_plot_data_types(edp_result), autoplot_types(res)) + + # A view name is chosen from a fixed menu like any other argument, and asking + # for the frame rather than the figure does not make the refusal a lesser one. + # The reader is shown the call they made rather than the method dispatch + # reached. + menu_error <- rlang::catch_cnd(pull_plot_data(res, type = "bogus")) + expect_s3_class(menu_error, "positively_args_error") + expect_match(deparsed_call(menu_error), "pull_plot_data", fixed = TRUE) +}) + +# ---- Frames against what the views draw ----------------------------------- + +test_that("the histogram frame is what the histogram view draws", { + data <- sim_pull_gaussian(60) + res <- check_edp( + data, + exposure, + x1, + values = c(0, 1), + exposure_type = "continuous" + ) + + frame <- pull_plot_data(res, type = "histogram") + + # The observation-grain results, keyed on the intervention label in grid + # order. A function intervention varies `value` from row to row, so the label + # is the only key that separates the interventions. Pinning the whole frame + # rather than a column at a time leaves no room for a column to be dropped, + # reordered, or quietly recomputed. + expected <- res@results + expected$intervention <- factor( + expected$intervention, + levels = names(res@params$values) + ) + expect_identical(frame, expected) + + plot <- ggplot2::autoplot(res, type = "histogram") + expect_identical(plot$data, frame) + + built <- ggplot2::ggplot_build(plot) + layout <- built$layout$layout + layer <- built$data[[1]] + # One panel per intervention, in level order, and every row of the frame lands + # in the bins of its own panel. + expect_identical( + as.character(layout$intervention), + levels(frame$intervention) + ) + panel_counts <- as.numeric(tapply(layer$count, layer$PANEL, sum)) + expect_identical(panel_counts, as.numeric(table(frame$intervention))) +}) + +test_that("the ecdf frame is what the ecdf view draws", { + data <- sim_pull_gaussian(60) + res <- check_edp( + data, + exposure, + x1, + values = c(0, 1), + exposure_type = "continuous" + ) + + frame <- pull_plot_data(res, type = "ecdf") + + expect_s3_class(frame$intervention, "factor") + expect_identical(levels(frame$intervention), names(res@params$values)) + expect_identical(nrow(frame), nrow(res@results)) + # Anchored to the results rather than only to the curve drawn below, so a + # frame that agreed with its own view but not with the diagnostic would + # still be caught. + expect_identical(frame$edp, res@results$edp) + + plot <- ggplot2::autoplot(res, type = "ecdf") + expect_identical(plot$data, frame) + + built <- ggplot2::ggplot_build(plot) + layer <- built$data[[1]] + # One curve per intervention, each stepping at exactly that intervention's EDP + # values. `stat_ecdf()` pads a curve's ends with infinities. + expect_length(unique(layer$group), nlevels(frame$intervention)) + drawn <- lapply(split(layer$x, layer$group), function(x) { + sort(x[is.finite(x)]) + }) + expected <- lapply( + split(frame$edp, frame$intervention), + function(edp) sort(unique(edp)) + ) + expect_equal(unname(drawn), unname(expected)) +}) + +test_that("the scatter frame is the full estimator results and covers both layers", { + local_quiet() + data <- sim_pull_mixed_support() + res <- check_edp( + data, + exposure, + s, + variant = "estimator", + values = 1, + categorical_similarity = 0 + ) + + frame <- pull_plot_data(res, type = "scatter") + + # The view splits on `ideal_weight` itself, so the whole results tibble goes + # back and the split the plot draws stays derivable from it. + expect_identical(frame, res@results) + n_finite <- sum(is.finite(frame$ideal_weight)) + n_infinite <- sum(is.infinite(frame$ideal_weight)) + expect_gt(n_finite, 0L) + expect_gt(n_infinite, 0L) + + built <- ggplot2::ggplot_build(ggplot2::autoplot(res, type = "scatter")) + layer_rows <- vapply(built$data, nrow, integer(1)) + expect_length(built$data, 2) + expect_setequal(layer_rows, c(n_finite, n_infinite)) + expect_identical(sum(layer_rows), nrow(frame)) + + # The two layers together draw the frame's rows and nothing else. + drawn <- do.call( + rbind, + lapply(built$data, function(layer) layer[c("x", "y")]) + ) + expect_equal(sort(drawn$x), sort(frame$edp_treatment)) + expect_equal(sort(drawn$y), sort(frame$edp_outcome)) +}) + +test_that("a function intervention keys the frame by its label", { + data <- sim_pull_gaussian(60) + res <- check_edp( + data, + exposure, + x1, + values = list(down = function(.x) .x - 1, 0), + exposure_type = "continuous" + ) + + frame <- pull_plot_data(res, type = "histogram") + + expect_identical(levels(frame$intervention), c("down", "0")) + expect_identical(as.integer(table(frame$intervention)), rep(nrow(data), 2L)) + # The shift gives nearly every observation its own intervened exposure, so + # `value` cannot separate the interventions and the label must. + shifted <- frame$value[frame$intervention == "down"] + expect_gt(length(unique(shifted)), 1L) + expect_identical(unique(frame$value[frame$intervention == "0"]), 0) +}) + +test_that("interventions sharing a label share a level", { + data <- sim_pull_gaussian(60) + # The vector path labels each intervention by its own value, so a repeated + # value gives two interventions one label. + res <- check_edp( + data, + exposure, + x1, + values = c(1, 1), + exposure_type = "continuous" + ) + + frame <- pull_plot_data(res) + + # The shared label is one level, and neither intervention's rows are dropped + # or turned into missing levels along the way. + expect_identical(levels(frame$intervention), "1") + expect_identical(nrow(frame), 2L * nrow(data)) + expect_false(anyNA(frame$intervention)) + + for (type in c("histogram", "ecdf")) { + expect_identical(pull_plot_data(res, type = type), frame) + } + + # Both observation views draw the merged rows together, one panel and one + # curve rather than two of either. + histogram <- ggplot2::ggplot_build(ggplot2::autoplot(res, type = "histogram")) + expect_identical(nrow(histogram$layout$layout), 1L) + expect_identical(sum(histogram$data[[1]]$count), as.numeric(nrow(frame))) + + ecdf <- ggplot2::ggplot_build(ggplot2::autoplot(res, type = "ecdf")) + expect_length(unique(ecdf$data[[1]]$group), 1L) +}) + +# ---- Gate parity ---------------------------------------------------------- + +test_that("pull_plot_data() and autoplot() refuse the scatter view alike", { + local_quiet() + data <- sim_pull_gaussian(60) + plain <- check_edp( + data, + exposure, + x1, + values = 0, + exposure_type = "continuous" + ) + + pull_error <- rlang::catch_cnd(pull_plot_data(plain, type = "scatter")) + plot_error <- rlang::catch_cnd(ggplot2::autoplot(plain, type = "scatter")) + + # One gate behind two entry points: the reader is told the same thing either + # way and is shown the call they made. + expect_identical(conditionMessage(pull_error), conditionMessage(plot_error)) + expect_match(deparsed_call(pull_error), "pull_plot_data", fixed = TRUE) + expect_match(deparsed_call(plot_error), "autoplot", fixed = TRUE) + + expect_snapshot_abort( + pull_plot_data(plain, type = "scatter"), + class = "positively_variant_error" + ) + expect_snapshot_abort( + ggplot2::autoplot(plain, type = "scatter"), + class = "positively_variant_error" + ) +}) From 0cb73a9e6c2c9d0813799ff7d624520e517b344d Mon Sep 17 00:00:00 2001 From: Malcolm Barrett Date: Fri, 21 Aug 2026 20:39:38 +0200 Subject: [PATCH 3/7] make the EDP boxplot the default view and add a density view --- DESCRIPTION | 1 + R/plot-edp.R | 171 +++++++++++- R/pull-plot-data.R | 104 +++++++- inst/WORDLIST | 1 + man/autoplot.edp_result.Rd | 18 +- man/pull_plot_data.Rd | 18 +- tests/testthat/_snaps/check-edp.md | 9 + .../check-edp/edp-boxplot-by-intervention.svg | 87 ++++++ .../check-edp/edp-density-by-intervention.svg | 55 ++++ .../edp-estimator-boxplot-by-measure.svg | 168 ++++++++++++ .../edp-histogram-shift-intervention.svg | 158 +++++++++++ .../positivity-check-panel.svg | 220 +++++---------- tests/testthat/test-check-edp.R | 251 ++++++++++++++++++ tests/testthat/test-pull-plot-data.R | 207 ++++++++++++++- 14 files changed, 1290 insertions(+), 178 deletions(-) create mode 100644 tests/testthat/_snaps/check-edp/edp-boxplot-by-intervention.svg create mode 100644 tests/testthat/_snaps/check-edp/edp-density-by-intervention.svg create mode 100644 tests/testthat/_snaps/check-edp/edp-estimator-boxplot-by-measure.svg create mode 100644 tests/testthat/_snaps/check-edp/edp-histogram-shift-intervention.svg diff --git a/DESCRIPTION b/DESCRIPTION index 70ba72a..d3e4d5a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -34,6 +34,7 @@ Imports: vctrs Suggests: dplyr (>= 1.1.0), + ggridges, halfmoon, knitr, lmtp, diff --git a/R/plot-edp.R b/R/plot-edp.R index c17bb4c..0e5540b 100644 --- a/R/plot-edp.R +++ b/R/plot-edp.R @@ -2,15 +2,22 @@ #' Plot an effective-data-points diagnostic #' -#' Draws one view of a [check_edp()] result. `type = "histogram"` shows the -#' distribution of EDP faceted by intervention, and `type = "ecdf"` shows its -#' empirical cumulative distribution colored by intervention. For the -#' estimator variant, `type = "scatter"` plots `edp_outcome` against +#' Draws one view of a [check_edp()] result. `type = "boxplot"`, the default, +#' boxes EDP by intervention. Its whiskers reach the fifth and the ninety-fifth +#' percentile, following Ring and Schomaker (2026) rather than Tukey's fences, +#' and the values outside them are drawn as points. For the estimator variant +#' it draws the outcome-model and the treatment-model measure side by side in a +#' facet per intervention. `type = "histogram"` shows the distribution of EDP +#' faceted by intervention, `type = "ecdf"` shows its empirical cumulative +#' distribution colored by intervention, and `type = "density"` draws one +#' density ridgeline per intervention, which needs the \pkg{ggridges} package. +#' For the estimator variant, `type = "scatter"` plots `edp_outcome` against #' `edp_treatment` colored by `ideal_weight`; it aborts for the data variant. #' [pull_plot_data()] returns the tibble any of these views draws. #' #' @param object An `edp_result` from [check_edp()]. -#' @param type One of `"histogram"`, `"ecdf"`, or `"scatter"`. +#' @param type One of `"boxplot"`, `"histogram"`, `"ecdf"`, `"density"`, or +#' `"scatter"`. #' @param ... Not used. #' #' @return A [ggplot2::ggplot] object. @@ -23,6 +30,7 @@ #' df <- data.frame(dose = dose, x1 = x1) #' result <- check_edp(df, dose, x1, values = c(0, 1), exposure_type = "continuous") #' +#' autoplot(result) #' autoplot(result, type = "histogram") #' autoplot(result, type = "ecdf") #' @@ -32,17 +40,97 @@ NULL method(autoplot, edp_result) <- function( object, - type = c("histogram", "ecdf", "scatter"), + type = c("boxplot", "histogram", "ecdf", "density", "scatter"), ... ) { type <- resolve_arg_match(rlang::arg_match(type)) - if (type == "histogram") { - autoplot_edp_histogram(object) - } else if (type == "ecdf") { - autoplot_edp_ecdf(object) - } else { - autoplot_edp_scatter(object, call = rlang::current_call()) + switch( + type, + boxplot = autoplot_edp_boxplot(object), + histogram = autoplot_edp_histogram(object), + ecdf = autoplot_edp_ecdf(object), + density = autoplot_edp_density(object, call = rlang::current_call()), + scatter = autoplot_edp_scatter(object, call = rlang::current_call()) + ) +} + +#' The boxplot view of an EDP diagnostic +#' +#' The boxes are drawn from summaries computed ahead of the plot, because the +#' whiskers reach the fifth and the ninety-fifth percentile rather than Tukey's +#' fences, and the values outside them are a point layer of their own. That +#' layer draws the `outliers` list-column unchopped, so it holds no rows at all +#' when every value lies within the whiskers. +#' +#' @param object An `edp_result`. +#' +#' @return A [ggplot2::ggplot] object. +#' @keywords internal +#' @noRd +autoplot_edp_boxplot <- function(object) { + estimator <- object@variant == "estimator" + plot_data <- pull_plot_data(object, type = "boxplot") + axis <- if (estimator) "measure" else "intervention" + plot <- ggplot2::ggplot( + plot_data, + ggplot2::aes( + x = .data[[axis]], + ymin = .data$ymin, + lower = .data$lower, + middle = .data$middle, + upper = .data$upper, + ymax = .data$ymax + ) + ) + + ggplot2::geom_boxplot(stat = "identity") + + ggplot2::geom_point( + data = edp_boxplot_outliers(plot_data), + mapping = ggplot2::aes(x = .data[[axis]], y = .data$outlier), + inherit.aes = FALSE + ) + if (!estimator) { + return( + plot + + ggplot2::labs( + x = "Intervention", + y = "Effective data points", + title = "Effective data points by intervention" + ) + ) } + plot + + ggplot2::facet_wrap( + ggplot2::vars(.data$intervention), + labeller = ggplot2::label_both + ) + + ggplot2::scale_x_discrete( + labels = c( + edp_outcome = "Outcome model", + edp_treatment = "Treatment model" + ) + ) + + ggplot2::labs( + x = "Model", + y = "Effective data points", + title = "Effective data points by model and intervention" + ) +} + +#' The values outside the whiskers, one row each +#' +#' @param plot_data The box statistics from [pull_plot_data()]. +#' +#' @return The keys of `plot_data` repeated once per outlier they carry, with +#' the outliers themselves in `outlier`. +#' @keywords internal +#' @noRd +edp_boxplot_outliers <- function(plot_data) { + statistics <- c("n", "ymin", "lower", "middle", "upper", "ymax", "outliers") + keys <- plot_data[setdiff(names(plot_data), statistics)] + outliers <- vctrs::vec_rep_each(keys, lengths(plot_data$outliers)) + # unlist() on a list of empty vectors returns NULL rather than an empty one. + outliers$outlier <- as.double(unlist(plot_data$outliers, use.names = FALSE)) + outliers } #' The histogram view of an EDP diagnostic @@ -91,6 +179,65 @@ autoplot_edp_ecdf <- function(object) { ) } +#' The ridgeline density view of an EDP diagnostic +#' +#' @param object An `edp_result`. +#' @param call The calling environment, used to build the error's call. +#' +#' @return A [ggplot2::ggplot] object. +#' @keywords internal +#' @noRd +autoplot_edp_density <- function(object, call = rlang::caller_env()) { + if (!rlang::is_installed("ggridges")) { + abort( + c( + "The density view requires the {.pkg ggridges} package.", + i = "Install {.pkg ggridges}, or draw another view: {.code autoplot(x, type = \"histogram\")}." + ), + error_class = "positively_missing_package_error", + call = call + ) + } + measure <- edp_primary_column(object@results) + plot_data <- pull_plot_data(object, type = "density") + ggplot2::ggplot( + plot_data, + ggplot2::aes(x = .data[[measure]], y = .data$intervention) + ) + + ggridges::geom_density_ridges( + bandwidth = edp_joint_bandwidth( + plot_data[[measure]], + plot_data$intervention + ) + ) + + ggplot2::labs( + x = "Effective data points", + y = "Intervention", + title = "Effective data points by intervention" + ) +} + +#' The bandwidth the ridgelines share +#' +#' `ggridges::geom_density_ridges()` gives every ridgeline in a panel one +#' bandwidth, the mean of the per-group `stats::bw.nrd0()` rules, and reports +#' the value it settled on. Reproducing the rule here and passing the result in +#' keeps that report out of every transcript that draws the view without +#' changing the figure. Each intervention contributes one row per observation +#' and `check_edp()` requires at least two observations, so no group is ever too +#' small for `stats::bw.nrd0()`. +#' +#' @param values The measure the ridgelines smooth. +#' @param intervention A factor grouping `values`, one ridgeline per level. +#' +#' @return The joint bandwidth, a single positive number. +#' @keywords internal +#' @noRd +edp_joint_bandwidth <- function(values, intervention) { + groups <- split(values, intervention) + mean(vapply(groups, stats::bw.nrd0, numeric(1))) +} + #' The estimator scatter view of an EDP diagnostic #' #' @param object An `edp_result`. diff --git a/R/pull-plot-data.R b/R/pull-plot-data.R index 4096c74..ec7c108 100644 --- a/R/pull-plot-data.R +++ b/R/pull-plot-data.R @@ -20,11 +20,19 @@ #' #' For a [check_edp()] result: #' -#' * `"histogram"` and `"ecdf"` return the results at their own grain, one row -#' per observation and intervention, with `intervention` as a factor whose -#' levels follow the order the interventions were given in. The label is what -#' separates the interventions rather than `value`, which a function -#' intervention varies from observation to observation. +#' * `"boxplot"` returns the box statistics themselves, one row per +#' intervention, or one row per measure and intervention for the estimator +#' variant. `ymin` and `ymax` are the fifth and the ninety-fifth percentile, +#' `lower`, `middle`, and `upper` are the quartiles, `n` is the number of +#' observations summarized, and `outliers` is a list-column of the values +#' lying outside the whiskers, which the view draws as points. The estimator variant covers +#' `edp_outcome` and `edp_treatment`; `ideal_weight` shares no scale with +#' them and is left to the scatter view. +#' * `"histogram"`, `"ecdf"`, and `"density"` return the results at their own +#' grain, one row per observation and intervention, with `intervention` as a +#' factor whose levels follow the order the interventions were given in. The +#' label is what separates the interventions rather than `value`, which a +#' function intervention varies from observation to observation. #' * `"scatter"` returns the whole estimator results tibble. The view draws the #' rows of finite `ideal_weight` in one layer and the infinite rows in #' another, and both are read off `ideal_weight`. It needs the estimator @@ -81,6 +89,86 @@ edp_observation_frame <- function(x) { results } +#' Box statistics at the fifth and ninety-fifth percentiles +#' +#' Ring and Schomaker (2026) end the whiskers at the fifth and the ninety-fifth +#' percentile rather than at Tukey's fences, so the summaries are computed here +#' and the view draws them as they stand. Type 7 is the rule +#' `stats::quantile()` and [ggplot2::geom_boxplot()] both take by default, so +#' the quartiles are the ones an ordinary boxplot of the same values would +#' show. +#' +#' @param values The measure to summarize. +#' @param intervention A factor grouping `values`, whose levels fix the rows and +#' their order. +#' +#' @return A tibble of one row per level of `intervention`, holding `n`, the +#' five box statistics, and the outliers as a list-column. +#' @keywords internal +#' @noRd +edp_box_stats <- function(values, intervention) { + groups <- unname(split(values, intervention)) + quantiles <- lapply(groups, function(group) { + stats::quantile( + group, + c(0.05, 0.25, 0.5, 0.75, 0.95), + type = 7, + names = FALSE + ) + }) + statistic <- function(position) { + vapply(quantiles, `[[`, numeric(1), position) + } + tibble::tibble( + n = lengths(groups), + ymin = statistic(1L), + lower = statistic(2L), + middle = statistic(3L), + upper = statistic(4L), + ymax = statistic(5L), + outliers = Map( + function(group, quantile) { + group[group < quantile[[1L]] | group > quantile[[5L]]] + }, + groups, + quantiles + ) + ) +} + +#' The box statistics one row per box +#' +#' @param x An `edp_result`. +#' +#' @return A tibble keyed on the intervention, and on the measure as well for +#' the estimator variant. +#' @keywords internal +#' @noRd +edp_boxplot_frame <- function(x) { + results <- edp_observation_frame(x) + interventions <- results$intervention + keys <- factor(levels(interventions), levels = levels(interventions)) + if (x@variant != "estimator") { + return(vctrs::vec_cbind( + tibble::tibble(intervention = keys), + edp_box_stats(results$edp, interventions) + )) + } + # `ideal_weight` is a ratio rather than a count of effective data points, so + # it shares no scale with the two EDP measures and no axis with them either. + measures <- c("edp_outcome", "edp_treatment") + boxes <- lapply(measures, function(measure) { + vctrs::vec_cbind( + tibble::tibble( + intervention = keys, + measure = factor(measure, levels = measures) + ), + edp_box_stats(results[[measure]], interventions) + ) + }) + vctrs::vec_rbind(!!!boxes) +} + #' Require the estimator variant for the scatter view #' #' The scatter view reads `edp_outcome`, `edp_treatment`, and `ideal_weight`, @@ -110,7 +198,7 @@ edp_require_estimator <- function(x, call = rlang::caller_env()) { method(pull_plot_data, edp_result) <- function( x, - type = c("histogram", "ecdf", "scatter"), + type = c("boxplot", "histogram", "ecdf", "density", "scatter"), ... ) { # The reader's own call rather than the dispatched method, which names a @@ -119,7 +207,9 @@ method(pull_plot_data, edp_result) <- function( rlang::arg_match(type), call = rlang::caller_call() ) - if (type == "scatter") { + if (type == "boxplot") { + edp_boxplot_frame(x) + } else if (type == "scatter") { edp_require_estimator(x, call = rlang::caller_call()) x@results } else { diff --git a/inst/WORDLIST b/inst/WORDLIST index b1570e0..e19d6cf 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -56,6 +56,7 @@ reframes renormalizes reweight reweights +ridgeline sPoRT stochasticity tibble diff --git a/man/autoplot.edp_result.Rd b/man/autoplot.edp_result.Rd index 1d0eee1..27fc339 100644 --- a/man/autoplot.edp_result.Rd +++ b/man/autoplot.edp_result.Rd @@ -7,7 +7,8 @@ \arguments{ \item{object}{An \code{edp_result} from \code{\link[=check_edp]{check_edp()}}.} -\item{type}{One of \code{"histogram"}, \code{"ecdf"}, or \code{"scatter"}.} +\item{type}{One of \code{"boxplot"}, \code{"histogram"}, \code{"ecdf"}, \code{"density"}, or +\code{"scatter"}.} \item{...}{Not used.} } @@ -15,10 +16,16 @@ A \link[ggplot2:ggplot]{ggplot2::ggplot} object. } \description{ -Draws one view of a \code{\link[=check_edp]{check_edp()}} result. \code{type = "histogram"} shows the -distribution of EDP faceted by intervention, and \code{type = "ecdf"} shows its -empirical cumulative distribution colored by intervention. For the -estimator variant, \code{type = "scatter"} plots \code{edp_outcome} against +Draws one view of a \code{\link[=check_edp]{check_edp()}} result. \code{type = "boxplot"}, the default, +boxes EDP by intervention. Its whiskers reach the fifth and the ninety-fifth +percentile, following Ring and Schomaker (2026) rather than Tukey's fences, +and the values outside them are drawn as points. For the estimator variant +it draws the outcome-model and the treatment-model measure side by side in a +facet per intervention. \code{type = "histogram"} shows the distribution of EDP +faceted by intervention, \code{type = "ecdf"} shows its empirical cumulative +distribution colored by intervention, and \code{type = "density"} draws one +density ridgeline per intervention, which needs the \pkg{ggridges} package. +For the estimator variant, \code{type = "scatter"} plots \code{edp_outcome} against \code{edp_treatment} colored by \code{ideal_weight}; it aborts for the data variant. \code{\link[=pull_plot_data]{pull_plot_data()}} returns the tibble any of these views draws. } @@ -30,6 +37,7 @@ dose <- rnorm(n, mean = x1) df <- data.frame(dose = dose, x1 = x1) result <- check_edp(df, dose, x1, values = c(0, 1), exposure_type = "continuous") +autoplot(result) autoplot(result, type = "histogram") autoplot(result, type = "ecdf") diff --git a/man/pull_plot_data.Rd b/man/pull_plot_data.Rd index c2f8819..6571f4e 100644 --- a/man/pull_plot_data.Rd +++ b/man/pull_plot_data.Rd @@ -33,11 +33,19 @@ is refused here exactly as \code{autoplot()} refuses it. For a \code{\link[=check_edp]{check_edp()}} result: \itemize{ -\item \code{"histogram"} and \code{"ecdf"} return the results at their own grain, one row -per observation and intervention, with \code{intervention} as a factor whose -levels follow the order the interventions were given in. The label is what -separates the interventions rather than \code{value}, which a function -intervention varies from observation to observation. +\item \code{"boxplot"} returns the box statistics themselves, one row per +intervention, or one row per measure and intervention for the estimator +variant. \code{ymin} and \code{ymax} are the fifth and the ninety-fifth percentile, +\code{lower}, \code{middle}, and \code{upper} are the quartiles, \code{n} is the number of +observations summarized, and \code{outliers} is a list-column of the values +lying outside the whiskers, which the view draws as points. The estimator variant covers +\code{edp_outcome} and \code{edp_treatment}; \code{ideal_weight} shares no scale with +them and is left to the scatter view. +\item \code{"histogram"}, \code{"ecdf"}, and \code{"density"} return the results at their own +grain, one row per observation and intervention, with \code{intervention} as a +factor whose levels follow the order the interventions were given in. The +label is what separates the interventions rather than \code{value}, which a +function intervention varies from observation to observation. \item \code{"scatter"} returns the whole estimator results tibble. The view draws the rows of finite \code{ideal_weight} in one layer and the infinite rows in another, and both are read off \code{ideal_weight}. It needs the estimator diff --git a/tests/testthat/_snaps/check-edp.md b/tests/testthat/_snaps/check-edp.md index 081ae9c..cd83718 100644 --- a/tests/testthat/_snaps/check-edp.md +++ b/tests/testthat/_snaps/check-edp.md @@ -7,6 +7,15 @@ ! The scatter view needs the estimator variant. i Rerun `check_edp()` with `variant = "estimator"`. +# the density view says what it needs when ggridges is absent + + Code + ggplot2::autoplot(res, type = "density") + Condition + Error in `autoplot.positively::edp_result`: + ! The density view requires the ggridges package. + i Install ggridges, or draw another view: `autoplot(x, type = "histogram")`. + # the argument validation messages are stable Code diff --git a/tests/testthat/_snaps/check-edp/edp-boxplot-by-intervention.svg b/tests/testthat/_snaps/check-edp/edp-boxplot-by-intervention.svg new file mode 100644 index 0000000..f19a07b --- /dev/null +++ b/tests/testthat/_snaps/check-edp/edp-boxplot-by-intervention.svg @@ -0,0 +1,87 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +10 +20 +30 + + + + + +0 +1 +Intervention +Effective data points +Effective data points by intervention + + diff --git a/tests/testthat/_snaps/check-edp/edp-density-by-intervention.svg b/tests/testthat/_snaps/check-edp/edp-density-by-intervention.svg new file mode 100644 index 0000000..cb35e77 --- /dev/null +++ b/tests/testthat/_snaps/check-edp/edp-density-by-intervention.svg @@ -0,0 +1,55 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + +0 +1 + + + + + + + +0 +10 +20 +30 +40 +Effective data points +Intervention +Effective data points by intervention + + diff --git a/tests/testthat/_snaps/check-edp/edp-estimator-boxplot-by-measure.svg b/tests/testthat/_snaps/check-edp/edp-estimator-boxplot-by-measure.svg new file mode 100644 index 0000000..154155c --- /dev/null +++ b/tests/testthat/_snaps/check-edp/edp-estimator-boxplot-by-measure.svg @@ -0,0 +1,168 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +intervention: 0 + + + + + + + + + + +intervention: 1 + + + + +Outcome model +Treatment model + + +Outcome model +Treatment model +0 +25 +50 +75 +100 + + + + + +Model +Effective data points +Effective data points by model and intervention + + diff --git a/tests/testthat/_snaps/check-edp/edp-histogram-shift-intervention.svg b/tests/testthat/_snaps/check-edp/edp-histogram-shift-intervention.svg new file mode 100644 index 0000000..b5812d1 --- /dev/null +++ b/tests/testthat/_snaps/check-edp/edp-histogram-shift-intervention.svg @@ -0,0 +1,158 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +intervention: down + + + + + + + + + + +intervention: 0 + + + + + + + +0 +10 +20 +30 +40 + + + + + +0 +10 +20 +30 +40 +0 +10 +20 +30 + + + + +Effective data points +Observations +Effective data points by intervention + + diff --git a/tests/testthat/_snaps/check-positivity/positivity-check-panel.svg b/tests/testthat/_snaps/check-positivity/positivity-check-panel.svg index 9e9b41a..4bd3b37 100644 --- a/tests/testthat/_snaps/check-positivity/positivity-check-panel.svg +++ b/tests/testthat/_snaps/check-positivity/positivity-check-panel.svg @@ -18,7 +18,7 @@ - + @@ -60,156 +60,82 @@ - - + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -intervention: 0 - - - - - - - - - - -intervention: 1 - - - - - -0 -20 -40 - - - -0 -20 -40 -0 -5 -10 -15 - - - - -Effective data points -Observations +0 +20 +40 + + + + + +0 +1 +Intervention +Effective data points Effective data points by intervention diff --git a/tests/testthat/test-check-edp.R b/tests/testthat/test-check-edp.R index 795c62f..1a40be0 100644 --- a/tests/testthat/test-check-edp.R +++ b/tests/testthat/test-check-edp.R @@ -1416,8 +1416,12 @@ test_that("autoplot() returns a ggplot for each data-variant type", { exposure_type = "continuous" ) + expect_s3_class(ggplot2::autoplot(res, type = "boxplot"), "ggplot") expect_s3_class(ggplot2::autoplot(res, type = "histogram"), "ggplot") expect_s3_class(ggplot2::autoplot(res, type = "ecdf"), "ggplot") + + skip_if_not_installed("ggridges") + expect_s3_class(ggplot2::autoplot(res, type = "density"), "ggplot") }) test_that("autoplot() rejects an unknown type as a classed error", { @@ -1451,6 +1455,132 @@ test_that("plot() draws the view and returns the result invisibly", { expect_identical(plot(res), res) }) +test_that("the default view boxes EDP by intervention", { + data <- sim_edp_gaussian(60) + res <- check_edp( + data, + exposure, + x1, + values = list(down = function(.x) .x - 1, 0), + exposure_type = "continuous" + ) + + plot <- ggplot2::autoplot(res) + frame <- pull_plot_data(res, type = "boxplot") + expect_identical(plot$data, frame) + + # The whiskers reach the 5th and the 95th percentile rather than Tukey's + # fences, so the box is drawn from summaries computed ahead of it and the + # points beyond the whiskers are a layer of their own. + expect_length(plot$layers, 2L) + expect_s3_class(plot$layers[[1]]$geom, "GeomBoxplot") + expect_s3_class(plot$layers[[1]]$stat, "StatIdentity") + expect_s3_class(plot$layers[[2]]$geom, "GeomPoint") + + built <- ggplot2::ggplot_build(plot) + box <- built$data[[1]] + expect_identical(nrow(box), nrow(frame)) + for (stat in c("ymin", "lower", "middle", "upper", "ymax")) { + expect_equal(box[[stat]], frame[[stat]]) + } + + # One panel, with the interventions along x in the order they were given. + expect_identical(nrow(built$layout$layout), 1L) + expect_identical(as.numeric(box$x), c(1, 2)) + expect_identical( + built$layout$panel_params[[1]]$x$get_labels(), + c("down", "0") + ) + + points <- built$data[[2]] + expect_gt(nrow(points), 0L) + expect_identical(nrow(points), length(unlist(frame$outliers))) + expect_equal( + unname(lapply(split(points$y, as.numeric(points$x)), sort)), + unname(lapply(frame$outliers, sort)) + ) +}) + +test_that("the boxplot draws an empty point layer when nothing lies outside", { + local_quiet() + # Two strata under exact categorical matching: one of five observations with + # three treated, one of three with one treated. At the treated target every + # EDP is the number treated in its own stratum, so the eight values tie + # five-deep at 3 and three-deep at 1. The type-7 5th and 95th percentiles + # fall inside those ties, which leaves nothing outside the whiskers. + data <- data.frame( + exposure = c(0L, 0L, 1L, 1L, 1L, 0L, 0L, 1L), + s = factor(rep(c("a", "b"), c(5L, 3L))) + ) + res <- check_edp( + data, + exposure, + s, + values = 1, + categorical_similarity = 0 + ) + frame <- pull_plot_data(res, type = "boxplot") + expect_identical(frame$outliers[[1]], double(0)) + + plot <- ggplot2::autoplot(res) + expect_no_condition(ggplot2::ggplot_build(plot)) + + built <- ggplot2::ggplot_build(plot) + expect_length(built$data, 2L) + expect_identical(nrow(built$data[[2]]), 0L) +}) + +test_that("the estimator boxplot draws both EDP measures by intervention", { + data <- sim_edp_gaussian(60) + res <- check_edp( + data, + exposure, + x1, + variant = "estimator", + values = list(down = function(.x) .x - 1, 0), + exposure_type = "continuous" + ) + + plot <- ggplot2::autoplot(res) + frame <- pull_plot_data(res, type = "boxplot") + expect_identical(plot$data, frame) + expect_s3_class(plot$layers[[1]]$geom, "GeomBoxplot") + expect_s3_class(plot$layers[[1]]$stat, "StatIdentity") + + built <- ggplot2::ggplot_build(plot) + # One facet per intervention, in the order the interventions were given, with + # the two measures along x. The scale's data values pin the drawn order; the + # display text the axis shows is the figure's business. + expect_identical( + as.character(built$layout$layout$intervention), + c("down", "0") + ) + measures <- built$layout$panel_params[[1]]$x$limits + expect_identical(measures, c("edp_outcome", "edp_treatment")) + + # The ideal weight shares no scale with the two EDP measures, so this view + # does not draw it. + expect_false("ideal_weight" %in% measures) + expect_false("ideal_weight" %in% names(plot$data)) + expect_false(any(vapply( + built$data, + function(layer) "ideal_weight" %in% names(layer), + logical(1) + ))) + + box <- built$data[[1]] + expect_identical(nrow(box), nrow(frame)) + drawn <- paste( + levels(frame$intervention)[as.integer(box$PANEL)], + levels(frame$measure)[as.numeric(box$x)] + ) + position <- match(paste(frame$intervention, frame$measure), drawn) + expect_false(anyNA(position)) + for (stat in c("ymin", "lower", "middle", "upper", "ymax")) { + expect_equal(box[[stat]][position], frame[[stat]]) + } +}) + test_that("the estimator scatter view is a ggplot and aborts for the data variant", { data <- sim_edp_gaussian(150) estimator <- check_edp( @@ -1601,6 +1731,60 @@ test_that("the ecdf view colors by the intervention label", { expect_length(unique(built$data[[1]]$colour), 2) }) +test_that("the density view draws one ridgeline per intervention", { + skip_if_not_installed("ggridges") + data <- sim_edp_gaussian(60) + res <- check_edp( + data, + exposure, + x1, + values = list(down = function(.x) .x - 1, 0), + exposure_type = "continuous" + ) + + plot <- ggplot2::autoplot(res, type = "density") + expect_identical(plot$data, pull_plot_data(res, type = "density")) + expect_s3_class(plot$layers[[1]]$geom, "GeomDensityRidges") + + # ggridges reports the joint bandwidth it settles on unless it is handed one, + # so a view that leaves it to pick talks in every transcript that draws it. + expect_no_message(ggplot2::ggplot_build(plot)) + + built <- ggplot2::ggplot_build(plot) + # The interventions run up y in the order they were given; EDP runs along x, + # the measure the histogram view draws. + expect_identical( + built$layout$panel_params[[1]]$y$get_labels(), + c("down", "0") + ) + x_range <- built$layout$panel_params[[1]]$x.range + expect_lte(x_range[[1]], min(res@results$edp)) + expect_gte(x_range[[2]], max(res@results$edp)) +}) + +test_that("the density view says what it needs when ggridges is absent", { + local_quiet() + data <- sim_edp_gaussian(60) + # The result is built before the mock, because resolving the exposure type + # asks whether a package is installed. + res <- check_edp( + data, + exposure, + x1, + values = c(0, 1), + exposure_type = "continuous" + ) + + # Mocked rather than skipped, so the reader without ggridges sees a message + # this suite has read, rather than one that fires only on machines it never + # runs on. + local_mocked_bindings(is_installed = function(...) FALSE, .package = "rlang") + expect_snapshot_abort( + ggplot2::autoplot(res, type = "density"), + class = "positively_missing_package_error" + ) +}) + test_that("EDP autoplot views render as expected", { local_quiet() announce_doppelganger( @@ -1661,6 +1845,73 @@ test_that("EDP autoplot views render as expected", { ) }) +test_that("EDP boxplot views render as expected", { + local_quiet() + announce_doppelganger( + "EDP boxplot by intervention", + "EDP estimator boxplot by measure" + ) + data <- sim_edp_gaussian(150) + res <- check_edp( + data, + exposure, + x1, + values = c(0, 1), + exposure_type = "continuous" + ) + # The figure has to show the outlier layer it pins. + expect_gt(length(unlist(pull_plot_data(res, type = "boxplot")$outliers)), 0L) + expect_doppelganger("EDP boxplot by intervention", ggplot2::autoplot(res)) + + estimator <- check_edp( + data, + exposure, + x1, + variant = "estimator", + values = c(0, 1), + exposure_type = "continuous" + ) + expect_doppelganger( + "EDP estimator boxplot by measure", + ggplot2::autoplot(estimator) + ) +}) + +test_that("the EDP density view renders as expected", { + local_quiet() + announce_doppelganger("EDP density by intervention") + skip_if_not_installed("ggridges") + data <- sim_edp_gaussian(150) + res <- check_edp( + data, + exposure, + x1, + values = c(0, 1), + exposure_type = "continuous" + ) + expect_doppelganger( + "EDP density by intervention", + ggplot2::autoplot(res, type = "density") + ) +}) + +test_that("a function intervention renders in the histogram view", { + local_quiet() + announce_doppelganger("EDP histogram shift intervention") + data <- sim_edp_gaussian(150) + res <- check_edp( + data, + exposure, + x1, + values = list(down = function(.x) .x - 1, 0), + exposure_type = "continuous" + ) + expect_doppelganger( + "EDP histogram shift intervention", + ggplot2::autoplot(res, type = "histogram") + ) +}) + # ---- Value and estimator-selection validation ----------------------------- test_that("check_edp() validates a user-supplied intervention-value grid", { diff --git a/tests/testthat/test-pull-plot-data.R b/tests/testthat/test-pull-plot-data.R index dc4d9b5..537b8a6 100644 --- a/tests/testthat/test-pull-plot-data.R +++ b/tests/testthat/test-pull-plot-data.R @@ -102,6 +102,10 @@ test_that("the type menu mirrors autoplot's and refuses an unknown view", { # The accessor exists to hand back what a view draws, so it offers exactly the # views autoplot() offers, in the same order. + expect_identical( + autoplot_types(res), + c("boxplot", "histogram", "ecdf", "density", "scatter") + ) expect_identical(pull_plot_data_types(edp_result), autoplot_types(res)) # A view name is chosen from a fixed menu like any other argument, and asking @@ -115,6 +119,205 @@ test_that("the type menu mirrors autoplot's and refuses an unknown view", { # ---- Frames against what the views draw ----------------------------------- +test_that("the boxplot frame is the hand-computed five-number summary", { + # Five observations whose exposure and covariate both run 0 to 4, with both + # half-distances at 1. Every kernel factor is then + # 0.5 ^ ((a_j - a*) ^ 2 + (x_j - x_i) ^ 2), a negative power of two, so each + # EDP is a sum of powers of two and is exact in double precision. + data <- tibble::tibble(exposure = c(0, 1, 2, 3, 4), x1 = c(0, 1, 2, 3, 4)) + res <- check_edp( + data, + exposure, + x1, + bw_exposure = 1, + bw_covariates = 1, + values = 0, + exposure_type = "continuous" + ) + edp <- c( + 2^0 + 2^-2 + 2^-8 + 2^-18 + 2^-32, + 2^0 + 2^-5 + 2^-13 + 2^-25, + 2^-2 + 2^-3 + 2^-10 + 2^-20, + 2^-4 + 2^-8 + 2^-17, + 2^-8 + 2^-9 + 2^-15 + ) + expect_identical(res@results$edp, edp) + + frame <- pull_plot_data(res, type = "boxplot") + + expect_identical( + names(frame), + c( + "intervention", + "n", + "ymin", + "lower", + "middle", + "upper", + "ymax", + "outliers" + ) + ) + expect_identical(nrow(frame), 1L) + expect_s3_class(frame$intervention, "factor") + expect_identical(levels(frame$intervention), "0") + expect_type(frame$n, "integer") + expect_identical(frame$n, 5L) + for (stat in c("ymin", "lower", "middle", "upper", "ymax")) { + expect_type(frame[[stat]], "double") + } + expect_type(frame$outliers, "list") + expect_type(frame$outliers[[1]], "double") + + # Type 7 reads a quantile at h = (n - 1)p + 1, which for these five + # observations is 1.2, 2, 3, 4, and 4.8: the quartiles and the median are + # order statistics outright, and each whisker interpolates a fifth of the way + # in from the end. The interpolated ends are compared rather than pinned bit + # for bit, because a fifth is not a binary fraction. + sorted <- sort(edp) + expect_equal( + frame$ymin, + sorted[[1]] + 0.2 * (sorted[[2]] - sorted[[1]]), + tolerance = 1e-12 + ) + expect_identical(frame$lower, sorted[[2]]) + expect_identical(frame$middle, sorted[[3]]) + expect_identical(frame$upper, sorted[[4]]) + expect_equal( + frame$ymax, + sorted[[4]] + 0.8 * (sorted[[5]] - sorted[[4]]), + tolerance = 1e-12 + ) + + # The smallest and the largest observation are the two lying strictly outside + # the whiskers. + expect_identical(sort(frame$outliers[[1]]), sorted[c(1L, 5L)]) +}) + +test_that("the boxplot frame is one row per intervention in grid order", { + data <- sim_pull_gaussian(60) + res <- check_edp( + data, + exposure, + x1, + values = list(down = function(.x) .x - 1, 0), + exposure_type = "continuous" + ) + + frame <- pull_plot_data(res, type = "boxplot") + + # The levels follow the order the interventions were given rather than the + # order their labels sort in, and so do the rows. + expect_identical( + levels(frame$intervention), + unique(names(res@params$values)) + ) + expect_identical(levels(frame$intervention), c("down", "0")) + expect_identical(as.character(frame$intervention), c("down", "0")) + expect_identical(frame$n, c(60L, 60L)) + + plot <- ggplot2::autoplot(res, type = "boxplot") + expect_identical(plot$data, frame) + + built <- ggplot2::ggplot_build(plot) + box <- built$data[[1]] + # One box per row of the frame, at the heights the frame gives rather than at + # the ones a boxplot would compute for itself. + expect_identical(nrow(box), nrow(frame)) + expect_identical(as.numeric(box$x), c(1, 2)) + for (stat in c("ymin", "lower", "middle", "upper", "ymax")) { + expect_equal(box[[stat]], frame[[stat]]) + } + + # The point layer draws the outliers list-column unchopped: every value, with + # its own box, and nothing besides. + points <- built$data[[2]] + expect_gt(nrow(points), 0L) + expect_identical(nrow(points), length(unlist(frame$outliers))) + expect_equal( + unname(lapply(split(points$y, as.numeric(points$x)), sort)), + unname(lapply(frame$outliers, sort)) + ) +}) + +test_that("the estimator boxplot frame covers the two EDP measures", { + data <- sim_pull_gaussian(60) + res <- check_edp( + data, + exposure, + x1, + variant = "estimator", + values = list(down = function(.x) .x - 1, 0), + exposure_type = "continuous" + ) + + frame <- pull_plot_data(res, type = "boxplot") + + expect_identical( + names(frame), + c( + "intervention", + "measure", + "n", + "ymin", + "lower", + "middle", + "upper", + "ymax", + "outliers" + ) + ) + expect_s3_class(frame$measure, "factor") + expect_s3_class(frame$intervention, "factor") + # The ideal weight shares no scale with the two EDP measures, so this view + # leaves it to the scatter view and to tidy(). + expect_identical(levels(frame$measure), c("edp_outcome", "edp_treatment")) + expect_identical(levels(frame$intervention), c("down", "0")) + + # One row per measure and intervention, and every pairing is present. + expect_identical(nrow(frame), 4L) + expect_setequal( + paste(frame$intervention, frame$measure), + c( + "down edp_outcome", + "0 edp_outcome", + "down edp_treatment", + "0 edp_treatment" + ) + ) + expect_identical(frame$n, rep(60L, 4L)) + expect_type(frame$outliers, "list") + + # Each row summarizes its own measure within its own intervention, which the + # median pins independently of the quantile call. + for (row in seq_len(nrow(frame))) { + measure <- as.character(frame$measure[[row]]) + label <- as.character(frame$intervention[[row]]) + values <- res@results[[measure]][res@results$intervention == label] + expect_identical(length(values), frame$n[[row]]) + expect_equal(frame$middle[[row]], stats::median(values)) + } +}) + +test_that("the density frame is the observation frame and needs no ggridges", { + data <- sim_pull_gaussian(60) + res <- check_edp( + data, + exposure, + x1, + values = c(0, 1), + exposure_type = "continuous" + ) + + frame <- pull_plot_data(res, type = "density") + expect_identical(frame, pull_plot_data(res, type = "histogram")) + + # The frame is the numbers rather than the figure, so the package that draws + # the ridgelines is not needed to read them. + local_mocked_bindings(is_installed = function(...) FALSE, .package = "rlang") + expect_identical(pull_plot_data(res, type = "density"), frame) +}) + test_that("the histogram frame is what the histogram view draws", { data <- sim_pull_gaussian(60) res <- check_edp( @@ -263,7 +466,7 @@ test_that("interventions sharing a label share a level", { exposure_type = "continuous" ) - frame <- pull_plot_data(res) + frame <- pull_plot_data(res, type = "histogram") # The shared label is one level, and neither intervention's rows are dropped # or turned into missing levels along the way. @@ -271,7 +474,7 @@ test_that("interventions sharing a label share a level", { expect_identical(nrow(frame), 2L * nrow(data)) expect_false(anyNA(frame$intervention)) - for (type in c("histogram", "ecdf")) { + for (type in c("histogram", "ecdf", "density")) { expect_identical(pull_plot_data(res, type = type), frame) } From 0e6487c800d245fe70cf727bcd34b4720d405ec2 Mon Sep 17 00:00:00 2001 From: Malcolm Barrett Date: Fri, 21 Aug 2026 22:00:05 +0200 Subject: [PATCH 4/7] add pull_plot_data() methods for the remaining diagnostics --- R/plot-density-ratios.R | 30 +- R/plot-eta-bias.R | 83 +-- R/plot-extrapolation.R | 29 +- R/plot-hat-values.R | 12 +- R/plot-hdr.R | 25 +- R/plot-port.R | 8 +- R/pull-plot-data.R | 505 ++++++++++++- man/autoplot.density_ratios_result.Rd | 2 +- man/autoplot.eta_bias_result.Rd | 1 + man/autoplot.extrapolation_result.Rd | 1 + man/autoplot.hat_values_result.Rd | 1 + man/autoplot.hdr_result.Rd | 3 +- man/autoplot.port_result.Rd | 3 +- man/pull_plot_data.Rd | 65 +- tests/testthat/_snaps/pull-plot-data.md | 88 +++ tests/testthat/test-pull-plot-data.R | 944 +++++++++++++++++++++++- 16 files changed, 1655 insertions(+), 145 deletions(-) diff --git a/R/plot-density-ratios.R b/R/plot-density-ratios.R index 1855c11..b31fcb6 100644 --- a/R/plot-density-ratios.R +++ b/R/plot-density-ratios.R @@ -7,7 +7,7 @@ #' for a point treatment and a per-time boxplot for a time-varying treatment. #' `type = "cumulative"` shows the cumulative-product quantiles across time #' points, the sequential positivity signature, and is available for matrix -#' inputs. +#' inputs. [pull_plot_data()] returns the tibble either of these views draws. #' #' @param object A `density_ratios_result` from [check_density_ratios()]. #' @param type One of `"distribution"` or `"cumulative"`. @@ -49,7 +49,7 @@ method(autoplot, density_ratios_result) <- function( #' @keywords internal #' @noRd autoplot_density_ratios_histogram <- function(object) { - plot_data <- tibble::tibble(ratio = object@ratios[[1]]) + plot_data <- pull_plot_data(object, type = "distribution") ggplot2::ggplot(plot_data, ggplot2::aes(x = .data$ratio)) + ggplot2::geom_histogram(bins = 30, fill = "grey70", color = "white") + ggplot2::labs( @@ -67,13 +67,7 @@ autoplot_density_ratios_histogram <- function(object) { #' @keywords internal #' @noRd autoplot_density_ratios_boxplot <- function(object) { - per_time <- purrr::imap( - object@ratios, - function(ratios, time) { - tibble::tibble(time = as.integer(time), ratio = ratios) - } - ) - plot_data <- vctrs::vec_rbind(!!!per_time) + plot_data <- pull_plot_data(object, type = "distribution") ggplot2::ggplot( plot_data, ggplot2::aes(x = factor(.data$time), y = .data$ratio) @@ -97,22 +91,8 @@ autoplot_density_ratios_cumulative <- function( object, call = rlang::caller_env() ) { - results <- object@results - cumulative <- results[grepl("^cumulative_", results$statistic), ] - if (nrow(cumulative) == 0) { - abort( - c( - "{.code type = \"cumulative\"} requires a time-varying (matrix) input.", - i = "This result summarizes a point treatment with a single time point." - ), - error_class = "positively_type_error", - call = call - ) - } - is_quantile <- grepl("^cumulative_quantile_", cumulative$statistic) - is_max <- cumulative$statistic == "cumulative_max" - plot_data <- cumulative[is_quantile | is_max, ] - plot_data$series <- cumulative_series(plot_data$statistic) + density_ratios_require_cumulative(object, call = call) + plot_data <- pull_plot_data(object, type = "cumulative") ggplot2::ggplot( plot_data, ggplot2::aes(x = .data$time, y = .data$value, color = .data$series) diff --git a/R/plot-eta-bias.R b/R/plot-eta-bias.R index 1d454ed..afc00d4 100644 --- a/R/plot-eta-bias.R +++ b/R/plot-eta-bias.R @@ -8,6 +8,7 @@ #' is aimed at. `type = "sweep"` shows ETA.Bias with a #' two-Monte-Carlo-standard-error band across the truncation sweep, the #' bias-variance tradeoff of weight truncation, with one line per estimand term. +#' [pull_plot_data()] returns the tibble either of these views draws. #' #' An estimand of one term, read at one truncation level, draws a single facet. #' A discrete exposure truncates a fitted probability, so its sweep is drawn @@ -132,33 +133,15 @@ eta_bias_level_index <- function(n_levels, n_terms) { #' @keywords internal #' @noRd autoplot_eta_bias_bootstrap <- function(object) { - results <- object@results - term_levels <- unique(results$term) + boot_data <- pull_plot_data(object, type = "bootstrap") n_terms <- length(object@truth) - n_levels <- eta_bias_n_levels(results, n_terms) - sweep_levels <- eta_bias_levels(object) - level_index <- eta_bias_level_index(n_levels, n_terms) - - boot_data <- do.call( - rbind, - lapply(seq_len(nrow(results)), function(cell) { - tibble::tibble( - term = results$term[[cell]], - level = level_index[[cell]], - estimate = object@boot_estimates[[cell]] - ) - }) - ) - # Both keys are factors so the panels follow the sweep and the order the terms - # were built in. Left as they are, the level index would sort lexically and - # put a tenth level ahead of a second one. - boot_data$term <- factor(boot_data$term, levels = term_levels) - boot_data$level <- factor(boot_data$level, levels = seq_len(n_levels)) - # The truth is a property of the term alone, so a term's facets carry the one - # line the estimator in them was aimed at. + # The truth rides with every draw, but a line per draw is a line drawn + # thousands of times over, so the lines are drawn from one row per term. The + # truth is a property of the term alone, so a term's facets carry the one line + # the estimator in them was aimed at. truth_data <- tibble::tibble( - term = factor(names(object@truth), levels = term_levels), + term = factor(names(object@truth), levels = levels(boot_data$term)), truth = unname(object@truth) ) @@ -170,7 +153,7 @@ autoplot_eta_bias_bootstrap <- function(object) { linetype = "dashed", inherit.aes = FALSE ) + - eta_bias_bootstrap_facets(n_terms, sweep_levels$labels) + + eta_bias_bootstrap_facets(n_terms, nlevels(boot_data$label)) + ggplot2::labs( x = "Bootstrap estimate", y = "Bootstrap draws", @@ -190,26 +173,23 @@ autoplot_eta_bias_bootstrap <- function(object) { #' keeps a level's terms in separate panels: pooling them would put the bootstrap #' distributions of two different estimands in one histogram. #' +#' The levels are faceted on the `label` the frame carries, in the order the +#' sweep ran, so the strips read a column a reader can see rather than a lookup +#' behind the figure. +#' #' @param n_terms The number of estimand terms. -#' @param labels The strip label for each truncation level. +#' @param n_levels The number of truncation levels. #' #' @return A [ggplot2::ggplot] facet specification. #' @keywords internal #' @noRd -eta_bias_bootstrap_facets <- function(n_terms, labels) { - strip_labeller <- ggplot2::labeller( - level = ggplot2::as_labeller(function(value) labels[as.integer(value)]) - ) - if (n_terms == 1 && length(labels) == 1) { - return(ggplot2::facet_wrap( - ggplot2::vars(.data$level), - labeller = strip_labeller - )) +eta_bias_bootstrap_facets <- function(n_terms, n_levels) { + if (n_terms == 1 && n_levels == 1) { + return(ggplot2::facet_wrap(ggplot2::vars(.data$label))) } ggplot2::facet_grid( rows = ggplot2::vars(.data$term), - cols = ggplot2::vars(.data$level), - labeller = strip_labeller + cols = ggplot2::vars(.data$label) ) } @@ -222,35 +202,10 @@ eta_bias_bootstrap_facets <- function(n_terms, labels) { #' @keywords internal #' @noRd autoplot_eta_bias_sweep <- function(object, call = rlang::caller_env()) { - results <- object@results + eta_bias_require_sweep(object, call = call) n_terms <- length(object@truth) - n_levels <- eta_bias_n_levels(results, n_terms) - if (n_levels == 1) { - # An estimand of several terms fills several rows at one truncation level, - # so the row count reads as a sweep that never ran and the level count is - # what decides there is nothing to draw. - entry <- if (object@exposure_type == "continuous") { - "quantile levels" - } else { - "lower bounds" - } - abort( - c( - "The sweep view needs a truncation sweep of more than one level.", - i = "Rerun {.fn check_eta_bias} with a {.arg truncation_grid} of multiple {entry}." - ), - error_class = "positively_sweep_absent_error", - call = call - ) - } sweep_levels <- eta_bias_levels(object) - sweep_data <- tibble::tibble( - term = factor(results$term, levels = unique(results$term)), - truncation = sweep_levels$values[eta_bias_level_index(n_levels, n_terms)], - bias = results$bias, - lower = results$bias - 2 * results$mc_se, - upper = results$bias + 2 * results$mc_se - ) + sweep_data <- pull_plot_data(object, type = "sweep") # One term needs no key to tell terms apart, so it keeps the neutral grey band # and the plain line. The colour and fill labels belong to those mappings; # naming an aesthetic the plot never maps draws an unknown-label message. diff --git a/R/plot-extrapolation.R b/R/plot-extrapolation.R index b254f56..4fbf3e5 100644 --- a/R/plot-extrapolation.R +++ b/R/plot-extrapolation.R @@ -7,6 +7,7 @@ #' faceted by exposure group, the share of the opposite group that lies nearby #' each observation. `type = "hull"` shows convex-hull membership as stacked bars #' per exposure group, and is available only when the hull test ran. +#' [pull_plot_data()] returns the tibble either of these views draws. #' #' @param object An `extrapolation_result` from [check_extrapolation()]. #' @param type One of `"distribution"` or `"hull"`. @@ -50,8 +51,7 @@ method(autoplot, extrapolation_result) <- function( #' @keywords internal #' @noRd autoplot_extrapolation_distribution <- function(object) { - plot_data <- object@results - plot_data$exposure <- factor(plot_data$exposure) + plot_data <- pull_plot_data(object, type = "distribution") ggplot2::ggplot( plot_data, ggplot2::aes(x = .data$frac_nearby) @@ -83,29 +83,8 @@ autoplot_extrapolation_distribution <- function(object) { #' @keywords internal #' @noRd autoplot_extrapolation_hull <- function(object, call = rlang::caller_env()) { - if (!object@hull_run) { - # A rerun with hull = TRUE only helps when numeric covariates were present; - # without any, the missing precondition is what to report instead. - advice <- if (isTRUE(object@params$n_numeric > 0)) { - "Rerun {.fn check_extrapolation} with {.code hull = TRUE}." - } else { - "The hull test needs at least one numeric covariate." - } - abort( - c( - "The convex-hull view needs the hull test to have run.", - i = advice - ), - error_class = "positively_hull_absent_error", - call = call - ) - } - plot_data <- object@results - plot_data$exposure <- factor(plot_data$exposure) - plot_data$membership <- factor( - ifelse(plot_data$in_hull, "Inside hull", "Outside hull"), - levels = c("Inside hull", "Outside hull") - ) + extrapolation_require_hull(object, call = call) + plot_data <- pull_plot_data(object, type = "hull") ggplot2::ggplot( plot_data, ggplot2::aes(x = .data$exposure, fill = .data$membership) diff --git a/R/plot-hat-values.R b/R/plot-hat-values.R index dc50493..9adf680 100644 --- a/R/plot-hat-values.R +++ b/R/plot-hat-values.R @@ -7,6 +7,7 @@ #' observed \eqn{\hat{\phi}} and the `conf_level` null quantile marked. #' `type = "profile"` shows the proportion of high-leverage candidates at each #' exposure percentile, the leverage profile across the exposure range. +#' [pull_plot_data()] returns the tibble either of these views draws. #' #' @param object A `hat_values_result` from [check_hat_values()]. #' @param type One of `"null"` or `"profile"`. @@ -51,11 +52,11 @@ method(autoplot, hat_values_result) <- function( #' @keywords internal #' @noRd autoplot_hat_values_null <- function(object) { - null_data <- tibble::tibble(phi = object@null_dist) + null_data <- pull_plot_data(object, type = "null") conf_level <- object@params$conf_level # Fewer bins for small null_reps keeps the histogram from reading as ragged; # capped at 30 so large null distributions are not oversmoothed. - bins <- min(30L, max(1L, ceiling(sqrt(length(object@null_dist))))) + bins <- min(30L, max(1L, ceiling(sqrt(nrow(null_data))))) ggplot2::ggplot(null_data, ggplot2::aes(x = .data$phi)) + ggplot2::geom_histogram(bins = bins, fill = "grey70", color = "white") + ggplot2::geom_vline( @@ -87,12 +88,7 @@ autoplot_hat_values_null <- function(object) { #' @keywords internal #' @noRd autoplot_hat_values_profile <- function(object) { - results <- object@results - fraction <- tapply(results$high_leverage, results$prob, mean) - profile <- tibble::tibble( - prob = as.double(names(fraction)), - fraction = as.double(fraction) - ) + profile <- pull_plot_data(object, type = "profile") ggplot2::ggplot(profile, ggplot2::aes(x = .data$prob, y = .data$fraction)) + ggplot2::geom_line() + ggplot2::geom_point() + diff --git a/R/plot-hdr.R b/R/plot-hdr.R index 29b4df1..6bf85cf 100644 --- a/R/plot-hdr.R +++ b/R/plot-hdr.R @@ -4,7 +4,8 @@ #' #' Draws the non-overlap ratio \eqn{\hat{\tau}(a)} against the target exposure #' value `a`, with a reference line at zero. For a sequential result the curves -#' are colored by time point. +#' are colored by time point. [pull_plot_data()] returns the tibble this view +#' draws. #' #' @param object An `hdr_result` from [check_hdr()] or [check_hdr_seq()]. #' @param ... Not used. @@ -26,24 +27,25 @@ NULL method(autoplot, hdr_result) <- function(object, ...) { - results <- object@results - if ("time" %in% names(results)) { - autoplot_hdr_sequential(results) + plot_data <- pull_plot_data(object) + if ("time" %in% names(plot_data)) { + autoplot_hdr_sequential(plot_data) } else { - autoplot_hdr_point(results) + autoplot_hdr_point(plot_data) } } #' The point view of an HDR non-overlap diagnostic #' -#' @param results The `@results` tibble of a point `hdr_result`. +#' @param plot_data The frame [pull_plot_data()] returns for a point +#' `hdr_result`. #' #' @return A [ggplot2::ggplot] object. #' @keywords internal #' @noRd -autoplot_hdr_point <- function(results) { +autoplot_hdr_point <- function(plot_data) { ggplot2::ggplot( - results, + plot_data, ggplot2::aes(x = .data$value, y = .data$nonoverlap) ) + ggplot2::geom_hline(yintercept = 0, linetype = "dashed") + @@ -58,14 +60,13 @@ autoplot_hdr_point <- function(results) { #' The sequential view of an HDR non-overlap diagnostic #' -#' @param results The `@results` tibble of a sequential `hdr_result`. +#' @param plot_data The frame [pull_plot_data()] returns for a sequential +#' `hdr_result`. #' #' @return A [ggplot2::ggplot] object. #' @keywords internal #' @noRd -autoplot_hdr_sequential <- function(results) { - plot_data <- results - plot_data$time <- factor(plot_data$time) +autoplot_hdr_sequential <- function(plot_data) { ggplot2::ggplot( plot_data, ggplot2::aes( diff --git a/R/plot-port.R b/R/plot-port.R index 092eb68..8bcedaa 100644 --- a/R/plot-port.R +++ b/R/plot-port.R @@ -9,7 +9,8 @@ #' resolved thresholds. A censoring run is faceted by both time point and row #' type, so the exposure and censoring subgroups are separated, the censoring #' facets plot censoring prevalence, and each facet's reference lines are drawn -#' from the thresholds that family was judged against. +#' from the thresholds that family was judged against. [pull_plot_data()] +#' returns the tibble this view draws. #' #' A full result can report hundreds of subgroups, which crowds the y axis. #' Setting `low_support_only = TRUE` draws only the low-support subgroups and @@ -44,7 +45,10 @@ method(autoplot, port_result) <- function( ... ) { validate_flag(low_support_only, arg_name = "low_support_only") - plot_data <- port_plot_data(object, low_support_only = low_support_only) + plot_data <- pull_plot_data(object, low_support_only = low_support_only) + # The reference lines are per facet rather than per row, so a facet whose bars + # were all filtered out keeps the pair of lines its subgroups were judged + # against instead of losing them with the rows. reference <- port_beta_reference(object) sequential <- "time" %in% names(object@results) has_type <- "type" %in% names(object@results) diff --git a/R/pull-plot-data.R b/R/pull-plot-data.R index ec7c108..97621a4 100644 --- a/R/pull-plot-data.R +++ b/R/pull-plot-data.R @@ -38,9 +38,66 @@ #' another, and both are read off `ideal_weight`. It needs the estimator #' variant, and asking for it after a data-variant run is an error. #' +#' For a [check_density_ratios()] result: +#' +#' * `"distribution"` returns one row per ratio and time point, `time` and +#' `ratio`. That is the grain the point-treatment histogram bins and the +#' time-varying boxplot summarizes, and the raw ratios live nowhere else in +#' the result. +#' * `"cumulative"` returns the cumulative-product summaries the series view +#' draws, one row per time point and statistic, with `series` naming the +#' summary in the words the key uses. It needs a time-varying input. +#' +#' For a [check_eta_bias()] result: +#' +#' * `"bootstrap"` returns one row per bootstrap draw, keyed on the estimand +#' `term`, on the truncation `level`, and on the `label` its facet strip +#' reads. Each draw carries the `truth` its term was aimed at alongside the +#' `estimate`, which the view draws as one line per term. +#' * `"sweep"` returns one row per term and truncation level, holding the +#' `truncation` the sweep is drawn against, the `bias`, and the `lower` and +#' `upper` ends of the band two Monte Carlo standard errors either side of +#' it. It needs a sweep of more than one level. +#' +#' For a [check_hat_values()] result: +#' +#' * `"null"` returns the null replicates in `phi`. The three numbers the +#' figure marks beside them, `null_quantile`, `phi_hat`, and `conf_level`, +#' ride as repeated columns. +#' * `"profile"` returns one row per exposure percentile the candidates were +#' built at, `prob` and the `fraction` of that percentile's candidates +#' reading as high leverage. The results hold one row per candidate, so the +#' aggregation is what this view adds. +#' +#' For a [check_port()] or [check_port_seq()] result there is one view, so the +#' method takes `low_support_only` where a menu would take `type`. It returns +#' one row per bar: the reported subgroups, with a leaf that more than one tree +#' reported collapsed to a single row, the display `label`, the `width` +#' proportional to subgroup size, and the pair of thresholds the row's +#' prevalence was judged against in `beta_lower` and `beta_upper`. +#' +#' For a [check_hdr()] or [check_hdr_seq()] result there is one view, and the +#' frame is the results with `time` as a factor where a sequential run carries +#' one, because the wave separates the curves rather than measuring anything. +#' +#' For a [check_extrapolation()] result: +#' +#' * `"distribution"` returns the results with `exposure` as a factor, which is +#' what panels the histograms. +#' * `"hull"` returns the same rows with `membership`, the reading of `in_hull` +#' that the stacked bars and the key both use. It needs the hull test to have +#' run. +#' +#' A [positivity_check] holds one frame per child rather than a frame of its +#' own, so name the diagnostic to pull: `pull_plot_data(check, "port")`. The +#' frame is that child's, unchanged, and a `type` passed alongside reaches it. +#' #' @param x A [positivity_diagnostic]. -#' @param ... Passed to methods. Every method takes `type`, the name of the view -#' whose data to return, matched against that class's [autoplot()] menu. +#' @param ... Passed to methods. A class drawing more than one view takes +#' `type`, the name of the view whose data to return, matched against that +#' class's [autoplot()] menu. A [check_port()] result takes +#' `low_support_only` as [autoplot()] does, and a [positivity_check] takes the +#' name of the diagnostic to pull. #' #' @return A [tibble][tibble::tibble] holding one view's data. #' @@ -216,3 +273,447 @@ method(pull_plot_data, edp_result) <- function( edp_observation_frame(x) } } + +# ---- Density ratios -------------------------------------------------------- + +#' The raw ratios, one row per ratio and time point +#' +#' `@ratios` is the only place the ratios themselves live, so this unchops them +#' into the grain both distribution figures read: a point treatment bins one +#' time point's ratios into a histogram, a time-varying treatment summarizes +#' each time point's own into a box. The `time` column stays either way, so one +#' view keeps one shape whichever figure it draws. +#' +#' @param x A `density_ratios_result`. +#' +#' @return A tibble of `time` and `ratio`. +#' @keywords internal +#' @noRd +density_ratios_distribution_frame <- function(x) { + per_time <- purrr::imap( + x@ratios, + function(ratios, time) { + tibble::tibble(time = as.integer(time), ratio = ratios) + } + ) + vctrs::vec_rbind(!!!per_time) +} + +#' The cumulative-product series across time points +#' +#' Only the quantiles and the maximum, because a mean or a proportion shares no +#' scale with them. +#' +#' @param x A `density_ratios_result` whose results carry cumulative rows. +#' +#' @return The cumulative quantile and maximum rows with a `series` label. +#' @keywords internal +#' @noRd +density_ratios_cumulative_frame <- function(x) { + results <- x@results + cumulative <- results[grepl("^cumulative_", results$statistic), ] + is_quantile <- grepl("^cumulative_quantile_", cumulative$statistic) + is_max <- cumulative$statistic == "cumulative_max" + plot_data <- cumulative[is_quantile | is_max, ] + plot_data$series <- cumulative_series(plot_data$statistic) + plot_data +} + +#' Require cumulative rows for the cumulative view +#' +#' A point treatment has one time point and so no product across time points to +#' follow. Both the plot and the accessor gate here, so a reader is told the +#' same thing whichever they asked for and is shown the call they made. +#' +#' @param x A `density_ratios_result`. +#' @param call The call to blame, or the calling environment. +#' +#' @return `NULL`, invisibly, when the result carries cumulative rows. +#' @keywords internal +#' @noRd +density_ratios_require_cumulative <- function(x, call = rlang::caller_env()) { + if (!any(grepl("^cumulative_", x@results$statistic))) { + abort( + c( + "{.code type = \"cumulative\"} requires a time-varying (matrix) input.", + i = "This result summarizes a point treatment with a single time point." + ), + error_class = "positively_type_error", + call = call + ) + } + invisible(NULL) +} + +method(pull_plot_data, density_ratios_result) <- function( + x, + type = c("distribution", "cumulative"), + ... +) { + call <- rlang::caller_call() + type <- resolve_arg_match(rlang::arg_match(type), call = call) + if (type == "cumulative") { + density_ratios_require_cumulative(x, call = call) + density_ratios_cumulative_frame(x) + } else { + density_ratios_distribution_frame(x) + } +} + +# ---- ETA.Bias -------------------------------------------------------------- + +#' The bootstrap draws, one row each +#' +#' `@boot_estimates` is the only place the draws live, so this unchops them and +#' repeats each reading's keys across its own draws. A draw is dropped per +#' reported reading when its estimate is not finite, so the readings can rest on +#' different numbers of draws and the keys are repeated by length rather than by +#' a shared count. +#' +#' The term, the level, and the level's label are all factors, so the panels +#' follow the sweep and the order the terms were built in. Left as an integer, +#' the level index would sort lexically and put a tenth level ahead of a second +#' one. The `label` is what the facet strips read, and the truth is a property +#' of the term alone, so every row of a term carries the one value the estimator +#' in it was aimed at. +#' +#' @param x An `eta_bias_result`. +#' +#' @return A tibble of `term`, `level`, `label`, `estimate`, and `truth`. +#' @keywords internal +#' @noRd +eta_bias_bootstrap_frame <- function(x) { + results <- x@results + n_terms <- length(x@truth) + n_levels <- eta_bias_n_levels(results, n_terms) + labels <- eta_bias_levels(x)$labels + draws <- lengths(x@boot_estimates) + level <- rep(eta_bias_level_index(n_levels, n_terms), times = draws) + term_levels <- unique(results$term) + frame <- tibble::tibble( + term = factor(rep(results$term, times = draws), levels = term_levels), + level = factor(level, levels = seq_len(n_levels)), + label = factor(labels[level], levels = labels), + estimate = unlist(x@boot_estimates, use.names = FALSE) + ) + frame$truth <- unname(x@truth[as.character(frame$term)]) + frame +} + +#' The truncation sweep with the band drawn around it +#' +#' @param x An `eta_bias_result` read over more than one truncation level. +#' +#' @return A tibble of `term`, `truncation`, `bias`, `lower`, and `upper`. +#' @keywords internal +#' @noRd +eta_bias_sweep_frame <- function(x) { + results <- x@results + n_terms <- length(x@truth) + n_levels <- eta_bias_n_levels(results, n_terms) + values <- eta_bias_levels(x)$values + tibble::tibble( + term = factor(results$term, levels = unique(results$term)), + truncation = values[eta_bias_level_index(n_levels, n_terms)], + bias = results$bias, + lower = results$bias - 2 * results$mc_se, + upper = results$bias + 2 * results$mc_se + ) +} + +#' Require a sweep of more than one level for the sweep view +#' +#' An estimand of several terms fills several rows at one truncation level, so +#' the row count reads as a sweep that never ran and the level count is what +#' decides there is nothing to draw. Both the plot and the accessor gate here, +#' so a reader is told the same thing whichever they asked for and is shown the +#' call they made. +#' +#' @param x An `eta_bias_result`. +#' @param call The call to blame, or the calling environment. +#' +#' @return `NULL`, invisibly, when the result was read over several levels. +#' @keywords internal +#' @noRd +eta_bias_require_sweep <- function(x, call = rlang::caller_env()) { + if (eta_bias_n_levels(x@results, length(x@truth)) == 1) { + entry <- if (x@exposure_type == "continuous") { + "quantile levels" + } else { + "lower bounds" + } + abort( + c( + "The sweep view needs a truncation sweep of more than one level.", + i = "Rerun {.fn check_eta_bias} with a {.arg truncation_grid} of multiple {entry}." + ), + error_class = "positively_sweep_absent_error", + call = call + ) + } + invisible(NULL) +} + +method(pull_plot_data, eta_bias_result) <- function( + x, + type = c("bootstrap", "sweep"), + ... +) { + call <- rlang::caller_call() + type <- resolve_arg_match(rlang::arg_match(type), call = call) + if (type == "bootstrap") { + eta_bias_bootstrap_frame(x) + } else { + eta_bias_require_sweep(x, call = call) + eta_bias_sweep_frame(x) + } +} + +# ---- Hat values ------------------------------------------------------------ + +#' The null replicates and the readings marked beside them +#' +#' `@null_dist` is the only place the replicates live. The observed reading, the +#' null quantile it is judged against, and the confidence level that quantile +#' was taken at are single numbers the figure marks as lines, and they ride as +#' repeated columns because a view returns one tibble whatever the figure does +#' with it. +#' +#' @param x A `hat_values_result`. +#' +#' @return A tibble of `phi`, `null_quantile`, `phi_hat`, and `conf_level`. +#' @keywords internal +#' @noRd +hat_values_null_frame <- function(x) { + tibble::tibble( + phi = x@null_dist, + null_quantile = x@null_quantile, + phi_hat = x@phi_hat, + conf_level = x@params$conf_level + ) +} + +#' The share of high-leverage candidates at each exposure percentile +#' +#' The results hold one row per candidate and percentile, so the aggregation is +#' the whole of what the profile view adds. The percentiles are taken from the +#' results themselves rather than from the names the grouping carries, which +#' round-trip through text and back to a number that is close to the percentile +#' asked for rather than equal to it. +#' +#' @param x A `hat_values_result`. +#' +#' @return A tibble of `prob` and `fraction`, one row per percentile. +#' @keywords internal +#' @noRd +hat_values_profile_frame <- function(x) { + results <- x@results + # split() groups on the sorted distinct values, so the group means come back + # in the order the percentiles are built in below. + groups <- split(results$high_leverage, results$prob) + tibble::tibble( + prob = sort(unique(results$prob)), + fraction = vapply(groups, mean, numeric(1), USE.NAMES = FALSE) + ) +} + +method(pull_plot_data, hat_values_result) <- function( + x, + type = c("null", "profile"), + ... +) { + type <- resolve_arg_match(rlang::arg_match(type), call = rlang::caller_call()) + if (type == "null") { + hat_values_null_frame(x) + } else { + hat_values_profile_frame(x) + } +} + +# ---- PoRT ------------------------------------------------------------------ + +#' The reported subgroups with the thresholds each was judged against +#' +#' The bars are the reported subgroups, a leaf reported by more than one tree +#' collapsed to one row; the pair of thresholds a row's prevalence was read +#' against comes last, because that pair is what shows why a subgroup reads as +#' low support. A sequential run resolves its thresholds per time point, and a +#' censoring run judges the exposure and censoring families on different risk +#' sets, so a row's thresholds follow its family as well as its wave. +#' +#' @param x A `port_result`. +#' @param low_support_only Keep only the low-support rows when `TRUE`. +#' +#' @return The plotting frame with `beta_lower` and `beta_upper` appended. +#' @keywords internal +#' @noRd +port_subgroup_frame <- function(x, low_support_only = FALSE) { + frame <- port_plot_data(x, low_support_only = low_support_only) + lower <- port_row_beta(x, frame) + frame$beta_lower <- lower + frame$beta_upper <- 1 - lower + frame +} + +#' The lower threshold each row was judged against +#' +#' @param x A `port_result`. +#' @param frame The rows to resolve a threshold for. +#' +#' @return A double vector, one entry per row of `frame`. +#' @keywords internal +#' @noRd +port_row_beta <- function(x, frame) { + if (!("time" %in% names(frame))) { + return(rep(x@beta[[1]], nrow(frame))) + } + if (!("type" %in% names(frame))) { + return(as.double(x@beta[frame$time])) + } + as.double(ifelse( + frame$type == "censoring", + x@censoring_beta[frame$time], + x@beta[frame$time] + )) +} + +method(pull_plot_data, port_result) <- function( + x, + low_support_only = FALSE, + ... +) { + validate_flag( + low_support_only, + arg_name = "low_support_only", + call = rlang::caller_call() + ) + port_subgroup_frame(x, low_support_only = low_support_only) +} + +# ---- HDR ------------------------------------------------------------------- + +#' The non-overlap results keyed on the wave +#' +#' A sequential run draws one curve per time point, so the wave is a factor +#' here where the results carry it as the integer index. Left numeric it would +#' colour the curves off a continuous scale, which reads a wave as a quantity. +#' +#' @param x An `hdr_result`. +#' +#' @return The results, with `time` as a factor where there is one. +#' @keywords internal +#' @noRd +hdr_curve_frame <- function(x) { + results <- x@results + if ("time" %in% names(results)) { + results$time <- factor(results$time) + } + results +} + +method(pull_plot_data, hdr_result) <- function(x, ...) { + hdr_curve_frame(x) +} + +# ---- Extrapolation --------------------------------------------------------- + +#' The results keyed on the exposure group +#' +#' The exposure panels the histograms rather than measuring anything, so it is a +#' factor here where the results carry the raw level. +#' +#' @param x An `extrapolation_result`. +#' +#' @return The results, with `exposure` as a factor. +#' @keywords internal +#' @noRd +extrapolation_distribution_frame <- function(x) { + results <- x@results + results$exposure <- factor(results$exposure) + results +} + +#' The results with hull membership read as the key reads it +#' +#' @param x An `extrapolation_result` whose hull test ran. +#' +#' @return The distribution frame with a `membership` factor appended. +#' @keywords internal +#' @noRd +extrapolation_hull_frame <- function(x) { + results <- extrapolation_distribution_frame(x) + results$membership <- factor( + ifelse(results$in_hull, "Inside hull", "Outside hull"), + levels = c("Inside hull", "Outside hull") + ) + results +} + +#' Require the hull test for the convex-hull view +#' +#' Both the plot and the accessor gate here, so a reader is told the same thing +#' whichever they asked for and is shown the call they made. +#' +#' @param x An `extrapolation_result`. +#' @param call The call to blame, or the calling environment. +#' +#' @return `NULL`, invisibly, when the hull test ran. +#' @keywords internal +#' @noRd +extrapolation_require_hull <- function(x, call = rlang::caller_env()) { + if (!x@hull_run) { + # A rerun with hull = TRUE only helps when numeric covariates were present; + # without any, the missing precondition is what to report instead. + advice <- if (isTRUE(x@params$n_numeric > 0)) { + "Rerun {.fn check_extrapolation} with {.code hull = TRUE}." + } else { + "The hull test needs at least one numeric covariate." + } + abort( + c( + "The convex-hull view needs the hull test to have run.", + i = advice + ), + error_class = "positively_hull_absent_error", + call = call + ) + } + invisible(NULL) +} + +method(pull_plot_data, extrapolation_result) <- function( + x, + type = c("distribution", "hull"), + ... +) { + call <- rlang::caller_call() + type <- resolve_arg_match(rlang::arg_match(type), call = call) + if (type == "distribution") { + extrapolation_distribution_frame(x) + } else { + extrapolation_require_hull(x, call = call) + extrapolation_hull_frame(x) + } +} + +# ---- Positivity check ------------------------------------------------------ + +method(pull_plot_data, positivity_check) <- function(x, diagnostic, ...) { + call <- rlang::caller_call() + # `autoplot()` reads a missing diagnostic as the whole panel and takes `NULL` + # for it, so the two spellings arrive here alike. There is no panel of frames, + # so both are refused rather than one of them guessed at. + if (missing(diagnostic) || is.null(diagnostic)) { + abort( + c( + "{.arg diagnostic} must name a diagnostic in this container.", + i = "Available diagnostics are {.val {names(x)}}." + ), + error_class = "positively_diagnostic_error", + call = call + ) + } + # extract_named_check() is what `$` and `[[` refuse an unknown name with, so + # every way of naming a child spells the refusal the same. + pull_plot_data(extract_named_check(x, diagnostic, call = call), ...) +} diff --git a/man/autoplot.density_ratios_result.Rd b/man/autoplot.density_ratios_result.Rd index 38fc654..1a07bda 100644 --- a/man/autoplot.density_ratios_result.Rd +++ b/man/autoplot.density_ratios_result.Rd @@ -20,7 +20,7 @@ Draws one of two views of a \code{\link[=check_density_ratios]{check_density_rat for a point treatment and a per-time boxplot for a time-varying treatment. \code{type = "cumulative"} shows the cumulative-product quantiles across time points, the sequential positivity signature, and is available for matrix -inputs. +inputs. \code{\link[=pull_plot_data]{pull_plot_data()}} returns the tibble either of these views draws. } \examples{ m <- matrix(c(1, 1, 2, 1, 3, 1, 4, 1), nrow = 2) diff --git a/man/autoplot.eta_bias_result.Rd b/man/autoplot.eta_bias_result.Rd index 24f333a..576f2b8 100644 --- a/man/autoplot.eta_bias_result.Rd +++ b/man/autoplot.eta_bias_result.Rd @@ -21,6 +21,7 @@ estimand term and truncation level, each facet marking the \code{truth} its term is aimed at. \code{type = "sweep"} shows ETA.Bias with a two-Monte-Carlo-standard-error band across the truncation sweep, the bias-variance tradeoff of weight truncation, with one line per estimand term. +\code{\link[=pull_plot_data]{pull_plot_data()}} returns the tibble either of these views draws. } \details{ An estimand of one term, read at one truncation level, draws a single facet. diff --git a/man/autoplot.extrapolation_result.Rd b/man/autoplot.extrapolation_result.Rd index 2022ccd..3da0694 100644 --- a/man/autoplot.extrapolation_result.Rd +++ b/man/autoplot.extrapolation_result.Rd @@ -20,6 +20,7 @@ Draws one of two views of a \code{\link[=check_extrapolation]{check_extrapolatio faceted by exposure group, the share of the opposite group that lies nearby each observation. \code{type = "hull"} shows convex-hull membership as stacked bars per exposure group, and is available only when the hull test ran. +\code{\link[=pull_plot_data]{pull_plot_data()}} returns the tibble either of these views draws. } \examples{ set.seed(1) diff --git a/man/autoplot.hat_values_result.Rd b/man/autoplot.hat_values_result.Rd index f396d69..9606fde 100644 --- a/man/autoplot.hat_values_result.Rd +++ b/man/autoplot.hat_values_result.Rd @@ -20,6 +20,7 @@ shows the null distribution of \eqn{\hat{\phi}} as a histogram with the observed \eqn{\hat{\phi}} and the \code{conf_level} null quantile marked. \code{type = "profile"} shows the proportion of high-leverage candidates at each exposure percentile, the leverage profile across the exposure range. +\code{\link[=pull_plot_data]{pull_plot_data()}} returns the tibble either of these views draws. } \examples{ set.seed(1) diff --git a/man/autoplot.hdr_result.Rd b/man/autoplot.hdr_result.Rd index 8411145..fac655b 100644 --- a/man/autoplot.hdr_result.Rd +++ b/man/autoplot.hdr_result.Rd @@ -15,7 +15,8 @@ A \link[ggplot2:ggplot]{ggplot2::ggplot} object. \description{ Draws the non-overlap ratio \eqn{\hat{\tau}(a)} against the target exposure value \code{a}, with a reference line at zero. For a sequential result the curves -are colored by time point. +are colored by time point. \code{\link[=pull_plot_data]{pull_plot_data()}} returns the tibble this view +draws. } \examples{ set.seed(1) diff --git a/man/autoplot.port_result.Rd b/man/autoplot.port_result.Rd index ed53473..3be29f5 100644 --- a/man/autoplot.port_result.Rd +++ b/man/autoplot.port_result.Rd @@ -23,7 +23,8 @@ A sequential result is faceted by time point, and each facet carries its own resolved thresholds. A censoring run is faceted by both time point and row type, so the exposure and censoring subgroups are separated, the censoring facets plot censoring prevalence, and each facet's reference lines are drawn -from the thresholds that family was judged against. +from the thresholds that family was judged against. \code{\link[=pull_plot_data]{pull_plot_data()}} +returns the tibble this view draws. } \details{ A full result can report hundreds of subgroups, which crowds the y axis. diff --git a/man/pull_plot_data.Rd b/man/pull_plot_data.Rd index 6571f4e..2d557fb 100644 --- a/man/pull_plot_data.Rd +++ b/man/pull_plot_data.Rd @@ -9,8 +9,11 @@ pull_plot_data(x, ...) \arguments{ \item{x}{A \link{positivity_diagnostic}.} -\item{...}{Passed to methods. Every method takes \code{type}, the name of the view -whose data to return, matched against that class's \code{\link[=autoplot]{autoplot()}} menu.} +\item{...}{Passed to methods. A class drawing more than one view takes +\code{type}, the name of the view whose data to return, matched against that +class's \code{\link[=autoplot]{autoplot()}} menu. A \code{\link[=check_port]{check_port()}} result takes +\code{low_support_only} as \code{\link[=autoplot]{autoplot()}} does, and a \link{positivity_check} takes the +name of the diagnostic to pull.} } \value{ A \link[tibble:tibble]{tibble} holding one view's data. @@ -51,6 +54,64 @@ rows of finite \code{ideal_weight} in one layer and the infinite rows in another, and both are read off \code{ideal_weight}. It needs the estimator variant, and asking for it after a data-variant run is an error. } + +For a \code{\link[=check_density_ratios]{check_density_ratios()}} result: +\itemize{ +\item \code{"distribution"} returns one row per ratio and time point, \code{time} and +\code{ratio}. That is the grain the point-treatment histogram bins and the +time-varying boxplot summarizes, and the raw ratios live nowhere else in +the result. +\item \code{"cumulative"} returns the cumulative-product summaries the series view +draws, one row per time point and statistic, with \code{series} naming the +summary in the words the key uses. It needs a time-varying input. +} + +For a \code{\link[=check_eta_bias]{check_eta_bias()}} result: +\itemize{ +\item \code{"bootstrap"} returns one row per bootstrap draw, keyed on the estimand +\code{term}, on the truncation \code{level}, and on the \code{label} its facet strip +reads. Each draw carries the \code{truth} its term was aimed at alongside the +\code{estimate}, which the view draws as one line per term. +\item \code{"sweep"} returns one row per term and truncation level, holding the +\code{truncation} the sweep is drawn against, the \code{bias}, and the \code{lower} and +\code{upper} ends of the band two Monte Carlo standard errors either side of +it. It needs a sweep of more than one level. +} + +For a \code{\link[=check_hat_values]{check_hat_values()}} result: +\itemize{ +\item \code{"null"} returns the null replicates in \code{phi}. The three numbers the +figure marks beside them, \code{null_quantile}, \code{phi_hat}, and \code{conf_level}, +ride as repeated columns. +\item \code{"profile"} returns one row per exposure percentile the candidates were +built at, \code{prob} and the \code{fraction} of that percentile's candidates +reading as high leverage. The results hold one row per candidate, so the +aggregation is what this view adds. +} + +For a \code{\link[=check_port]{check_port()}} or \code{\link[=check_port_seq]{check_port_seq()}} result there is one view, so the +method takes \code{low_support_only} where a menu would take \code{type}. It returns +one row per bar: the reported subgroups, with a leaf that more than one tree +reported collapsed to a single row, the display \code{label}, the \code{width} +proportional to subgroup size, and the pair of thresholds the row's +prevalence was judged against in \code{beta_lower} and \code{beta_upper}. + +For a \code{\link[=check_hdr]{check_hdr()}} or \code{\link[=check_hdr_seq]{check_hdr_seq()}} result there is one view, and the +frame is the results with \code{time} as a factor where a sequential run carries +one, because the wave separates the curves rather than measuring anything. + +For a \code{\link[=check_extrapolation]{check_extrapolation()}} result: +\itemize{ +\item \code{"distribution"} returns the results with \code{exposure} as a factor, which is +what panels the histograms. +\item \code{"hull"} returns the same rows with \code{membership}, the reading of \code{in_hull} +that the stacked bars and the key both use. It needs the hull test to have +run. +} + +A \link{positivity_check} holds one frame per child rather than a frame of its +own, so name the diagnostic to pull: \code{pull_plot_data(check, "port")}. The +frame is that child's, unchanged, and a \code{type} passed alongside reaches it. } \examples{ set.seed(1) diff --git a/tests/testthat/_snaps/pull-plot-data.md b/tests/testthat/_snaps/pull-plot-data.md index bb57815..d6309da 100644 --- a/tests/testthat/_snaps/pull-plot-data.md +++ b/tests/testthat/_snaps/pull-plot-data.md @@ -16,3 +16,91 @@ ! The scatter view needs the estimator variant. i Rerun `check_edp()` with `variant = "estimator"`. +# pull_plot_data() and autoplot() refuse the cumulative view alike + + Code + pull_plot_data(res, type = "cumulative") + Condition + Error in `pull_plot_data()`: + ! `type = "cumulative"` requires a time-varying (matrix) input. + i This result summarizes a point treatment with a single time point. + +--- + + Code + ggplot2::autoplot(res, type = "cumulative") + Condition + Error in `autoplot.positively::density_ratios_result`: + ! `type = "cumulative"` requires a time-varying (matrix) input. + i This result summarizes a point treatment with a single time point. + +# pull_plot_data() and autoplot() refuse the sweep view alike + + Code + pull_plot_data(res, type = "sweep") + Condition + Error in `pull_plot_data()`: + ! The sweep view needs a truncation sweep of more than one level. + i Rerun `check_eta_bias()` with a `truncation_grid` of multiple lower bounds. + +--- + + Code + ggplot2::autoplot(res, type = "sweep") + Condition + Error in `autoplot.positively::eta_bias_result`: + ! The sweep view needs a truncation sweep of more than one level. + i Rerun `check_eta_bias()` with a `truncation_grid` of multiple lower bounds. + +# pull_plot_data() and autoplot() refuse a non-flag alike + + Code + pull_plot_data(res, low_support_only = "yes") + Condition + Error in `pull_plot_data()`: + ! `low_support_only` must be `TRUE` or `FALSE`. + +--- + + Code + ggplot2::autoplot(res, low_support_only = "yes") + Condition + Error in `ggplot2::autoplot()`: + ! `low_support_only` must be `TRUE` or `FALSE`. + +# pull_plot_data() and autoplot() refuse the hull view alike + + Code + pull_plot_data(res, type = "hull") + Condition + Error in `pull_plot_data()`: + ! The convex-hull view needs the hull test to have run. + i Rerun `check_extrapolation()` with `hull = TRUE`. + +--- + + Code + ggplot2::autoplot(res, type = "hull") + Condition + Error in `autoplot.positively::extrapolation_result`: + ! The convex-hull view needs the hull test to have run. + i Rerun `check_extrapolation()` with `hull = TRUE`. + +# a positivity check refuses to guess which diagnostic to pull + + Code + pull_plot_data(check) + Condition + Error in `pull_plot_data()`: + ! `diagnostic` must name a diagnostic in this container. + i Available diagnostics are "port" and "edp". + +--- + + Code + pull_plot_data(check, "bogus") + Condition + Error in `pull_plot_data()`: + ! "bogus" is not a diagnostic in this container. + i Available diagnostics are "port" and "edp". + diff --git a/tests/testthat/test-pull-plot-data.R b/tests/testthat/test-pull-plot-data.R index 537b8a6..9c49d70 100644 --- a/tests/testthat/test-pull-plot-data.R +++ b/tests/testthat/test-pull-plot-data.R @@ -34,9 +34,12 @@ sim_pull_mixed_support <- function() { # pass against a view the plot no longer offers. The S7 methods on `autoplot()` # register as S3 methods on ggplot2's generic; `pull_plot_data()` is an S7 # generic and its methods are reached through S7. +autoplot_method <- function(object) { + utils::getS3method("autoplot", class(object)[[1]]) +} + autoplot_types <- function(object) { - method <- utils::getS3method("autoplot", class(object)[[1]]) - eval(formals(method)$type) + eval(formals(autoplot_method(object))$type) } pull_plot_data_types <- function(class) { @@ -47,6 +50,36 @@ deparsed_call <- function(condition) { paste(deparse(conditionCall(condition)), collapse = " ") } +# ---- Reading a built plot ------------------------------------------------- + +# Built layers are read by what draws them rather than by the order they were +# added in: a histogram or a bar carries a count, a reference line an intercept, +# points a shape, a box its outliers, a ribbon a ymax. +layers_drawing <- function(built, aesthetic) { + Filter(function(layer) aesthetic %in% names(layer), built$data) +} + +layer_drawing <- function(built, aesthetic) { + layers_drawing(built, aesthetic)[[1]] +} + +# Every strip label a faceted plot renders, with its labeller applied. The +# labels are read off the built layout, one facet dimension at a time, so what +# a block compares against is the text the figure draws rather than the lookup +# behind it. +rendered_strip_labels <- function(plot) { + built <- ggplot2::ggplot_build(plot) + layout <- built$layout$layout + # ggplot2 spells its own layout bookkeeping columns in upper case, so dropping + # those leaves the variables the faceting was built on, whatever they are + # named. + keys <- layout[!grepl("^[A-Z_]+$", names(layout))] + labeller <- match.fun(built$layout$facet_params$labeller) + unlist(lapply(names(keys), function(key) { + as.character(unlist(labeller(unique(keys[key])))) + })) +} + # ---- The generic ---------------------------------------------------------- test_that("pull_plot_data() is an exported generic with an edp_result method", { @@ -519,3 +552,910 @@ test_that("pull_plot_data() and autoplot() refuse the scatter view alike", { class = "positively_variant_error" ) }) + +# ---- Density ratios ------------------------------------------------------- + +# The exact Gaussian shift-MTP density ratio, log r = k Z - k ^ 2 / 2. A vector +# of these is a point treatment; a matrix is a time-varying treatment, one +# column per time point. +sim_pull_ratios <- function(n, k = 0.5, seed = 1) { + withr::local_seed(seed) + exp(k * stats::rnorm(n) - k^2 / 2) +} + +test_that("the density-ratio menu mirrors autoplot's", { + res <- check_density_ratios(sim_pull_ratios(40)) + + expect_identical(autoplot_types(res), c("distribution", "cumulative")) + expect_identical( + pull_plot_data_types(density_ratios_result), + autoplot_types(res) + ) + expect_identical( + pull_plot_data(res), + pull_plot_data(res, type = "distribution") + ) +}) + +test_that("the distribution frame unchops the ratios by time", { + m <- matrix(sim_pull_ratios(60), nrow = 20, ncol = 3) + res <- check_density_ratios(m) + + frame <- pull_plot_data(res, type = "distribution") + + # One row per ratio per time point, the grain the boxplot view summarizes and + # the histogram view bins. `@ratios` is the only place the raw ratios live, so + # nothing else can hand a reader the numbers behind either figure. + expect_identical(names(frame), c("time", "ratio")) + expect_type(frame$time, "integer") + expect_type(frame$ratio, "double") + expect_identical(nrow(frame), 60L) + expect_identical(frame$time, rep(1:3, each = 20L)) + expect_identical(frame$ratio, unlist(res@ratios)) + + plot <- ggplot2::autoplot(res, type = "distribution") + expect_identical(plot$data, frame) + + built <- ggplot2::ggplot_build(plot) + box <- layer_drawing(built, "outliers") + # One box per time point, at the positions the times take, each spanning that + # time's ratios with its outliers included. + expect_identical(nrow(box), 3L) + expect_identical(as.numeric(box$x), c(1, 2, 3)) + by_time <- split(frame$ratio, frame$time) + expect_equal(box$ymin_final, unname(vapply(by_time, min, numeric(1)))) + expect_equal(box$ymax_final, unname(vapply(by_time, max, numeric(1)))) + expect_equal(box$middle, unname(vapply(by_time, stats::median, numeric(1)))) +}) + +test_that("the distribution frame is what the point histogram draws", { + res <- check_density_ratios(sim_pull_ratios(50, k = 1)) + + frame <- pull_plot_data(res, type = "distribution") + + # A point treatment is one time point, and the frame keeps the time column + # either way, so one view has one shape whichever figure it draws. + expect_identical(names(frame), c("time", "ratio")) + expect_identical(frame$time, rep(1L, 50L)) + expect_identical(frame$ratio, res@ratios[[1]]) + + plot <- ggplot2::autoplot(res, type = "distribution") + expect_identical(plot$data, frame) + + built <- ggplot2::ggplot_build(plot) + bins <- layer_drawing(built, "count") + # Every ratio is binned once and the bins span the frame, so the histogram + # accounts for the frame's rows and for nothing besides. + expect_identical(sum(bins$count), as.numeric(nrow(frame))) + expect_gte(min(frame$ratio), min(bins$xmin)) + expect_lte(max(frame$ratio), max(bins$xmax)) +}) + +test_that("the cumulative frame is the series the view draws", { + m <- matrix(sim_pull_ratios(60), nrow = 20, ncol = 3) + res <- check_density_ratios(m) + + frame <- pull_plot_data(res, type = "cumulative") + + expect_identical(names(frame), c("time", "statistic", "value", "series")) + expect_type(frame$time, "integer") + expect_type(frame$statistic, "character") + expect_type(frame$value, "double") + expect_s3_class(frame$series, "factor") + # The series read as text and are ordered by the quantile behind them, with + # the maximum last rather than where its label happens to sort. + expect_identical( + levels(frame$series), + c( + "50th percentile", + "90th percentile", + "95th percentile", + "99th percentile", + "Maximum" + ) + ) + # Only the cumulative rows, and among those only the quantiles and the + # maximum: a mean or a proportion shares no scale with them. + expect_true(all(startsWith(frame$statistic, "cumulative_"))) + expect_identical(nrow(frame), 15L) + + plot <- ggplot2::autoplot(res, type = "cumulative") + expect_identical(plot$data, frame) + + built <- ggplot2::ggplot_build(plot) + points <- layer_drawing(built, "shape") + # One point per row and one group per series, each at the time and the value + # the frame gives it. + expect_identical(nrow(points), nrow(frame)) + expect_length(unique(points$group), nlevels(frame$series)) + expect_equal(sort(points$x), sort(as.numeric(frame$time))) + expect_equal( + unname(lapply(split(points$y, points$group), sort)), + unname(lapply(split(frame$value, frame$series), sort)) + ) +}) + +test_that("pull_plot_data() and autoplot() refuse the cumulative view alike", { + res <- check_density_ratios(sim_pull_ratios(40, k = 1)) + + pull_error <- rlang::catch_cnd(pull_plot_data(res, type = "cumulative")) + plot_error <- rlang::catch_cnd(ggplot2::autoplot(res, type = "cumulative")) + + # One gate behind two entry points: the reader is told the same thing either + # way and is shown the call they made. + expect_identical(conditionMessage(pull_error), conditionMessage(plot_error)) + expect_match(deparsed_call(pull_error), "pull_plot_data", fixed = TRUE) + expect_match(deparsed_call(plot_error), "autoplot", fixed = TRUE) + + expect_snapshot_abort( + pull_plot_data(res, type = "cumulative"), + class = "positively_type_error" + ) + expect_snapshot_abort( + ggplot2::autoplot(res, type = "cumulative"), + class = "positively_type_error" + ) +}) + +# ---- ETA.Bias ------------------------------------------------------------- + +# A steep propensity model pushes fitted scores toward zero and one, the +# practical violation ETA.Bias is aimed at. The outcome is linear in the +# exposure, so the estimand holds one term at every truncation level. +sim_pull_eta <- function(n = 300, seed = 1) { + withr::local_seed(seed) + x1 <- stats::rnorm(n) + x2 <- stats::rnorm(n) + a <- stats::rbinom(n, 1L, stats::plogis(2 * (x1 + x2))) + tibble::tibble(a = a, y = a + x1 + x2 + stats::rnorm(n), x1 = x1, x2 = x2) +} + +# The same mechanism with a continuous exposure, which caps a stabilized weight +# rather than bounding a probability and so reports no lower bound at any level. +sim_pull_eta_continuous <- function(n = 300, seed = 1) { + withr::local_seed(seed) + x1 <- stats::rnorm(n) + x2 <- stats::rnorm(n) + a <- x1 + x2 + stats::rnorm(n) + tibble::tibble(a = a, y = a + x1 + x2 + stats::rnorm(n), x1 = x1, x2 = x2) +} + +# The bootstrap draw is seeded so the estimates a frame carries are the ones the +# figure drawn beside it holds. +fit_pull_eta <- function(data = sim_pull_eta(), ..., seed = 2024) { + withr::local_seed(seed) + check_eta_bias(data, a, y, c(x1, x2), ...) +} + +test_that("the eta-bias menu mirrors autoplot's", { + local_quiet() + res <- fit_pull_eta(n_boot = 10) + + expect_identical(autoplot_types(res), c("bootstrap", "sweep")) + expect_identical(pull_plot_data_types(eta_bias_result), autoplot_types(res)) + expect_identical(pull_plot_data(res), pull_plot_data(res, type = "bootstrap")) +}) + +test_that("the bootstrap frame unchops the bootstrap estimates", { + local_quiet() + res <- fit_pull_eta(truncation_grid = c(0, 0.05, 0.1), n_boot = 20) + + frame <- pull_plot_data(res, type = "bootstrap") + + expect_identical( + names(frame), + c("term", "level", "label", "estimate", "truth") + ) + expect_s3_class(frame$term, "factor") + expect_s3_class(frame$level, "factor") + expect_s3_class(frame$label, "factor") + expect_type(frame$estimate, "double") + expect_type(frame$truth, "double") + + # One row per draw per reported reading. `@boot_estimates` is the only place + # the draws live, so nothing else hands a reader the distributions the facets + # show. + expect_identical(nrow(frame), sum(lengths(res@boot_estimates))) + expect_identical(frame$estimate, unlist(res@boot_estimates)) + expect_identical(levels(frame$term), unique(res@results$term)) + # Both keys are factors so the panels follow the sweep and the order the terms + # were built in. Left as text the level index would sort lexically and put a + # tenth level ahead of a second one. + expect_identical(levels(frame$level), c("1", "2", "3")) + expect_identical( + levels(frame$label), + c("lower = 0", "lower = 0.05", "lower = 0.1") + ) + # The truth is a property of the term alone, so every row of a term carries + # the one value the estimator in it was aimed at. + expect_identical(frame$truth, unname(res@truth[as.character(frame$term)])) + + plot <- ggplot2::autoplot(res, type = "bootstrap") + expect_identical(plot$data, frame) + + built <- ggplot2::ggplot_build(plot) + + # The strips read the labels the frame carries, so a reader rebuilding the + # figure keys the facets on a column rather than on a lookup they cannot see. + # The term strip rides along in the same list, so the labels are counted + # rather than compared to the whole of it. + strips <- rendered_strip_labels(plot) + expect_true(all(levels(frame$label) %in% strips)) + expect_identical( + sum(strips %in% levels(frame$label)), + nlevels(frame$label) + ) + + bins <- layer_drawing(built, "count") + # Every draw is binned once, inside the panel of its own level. + expect_identical(sum(bins$count), as.numeric(nrow(frame))) + expect_identical( + unname(vapply(split(bins$count, bins$PANEL), sum, numeric(1))), + as.numeric(table(frame$level)) + ) + + # One truth per term, so the lines draw the frame's distinct truths and no + # target a panel's estimator was never aimed at. + lines <- layer_drawing(built, "xintercept") + expect_setequal(lines$xintercept, unique(frame$truth)) +}) + +test_that("the sweep frame is the band the sweep view draws", { + local_quiet() + res <- fit_pull_eta(truncation_grid = c(0, 0.05, 0.1), n_boot = 20) + + frame <- pull_plot_data(res, type = "sweep") + + expect_identical( + names(frame), + c("term", "truncation", "bias", "lower", "upper") + ) + expect_s3_class(frame$term, "factor") + expect_identical(levels(frame$term), unique(res@results$term)) + expect_type(frame$truncation, "double") + expect_identical(nrow(frame), nrow(res@results)) + expect_identical(frame$bias, res@results$bias) + # The band is two Monte Carlo standard errors either side of the bias, which + # is the interval the figure fills rather than one derived from it. + expect_identical(frame$lower, res@results$bias - 2 * res@results$mc_se) + expect_identical(frame$upper, res@results$bias + 2 * res@results$mc_se) + # A discrete exposure truncates a fitted probability, so its sweep is drawn + # against the lower bound it applied. + expect_identical(frame$truncation, res@results$truncation_lower) + + plot <- ggplot2::autoplot(res, type = "sweep") + expect_identical(plot$data, frame) + + built <- ggplot2::ggplot_build(plot) + points <- layer_drawing(built, "shape") + expect_equal(points$x, frame$truncation) + expect_equal(points$y, frame$bias) + + band <- layer_drawing(built, "ymax") + expect_equal(band$x, frame$truncation) + expect_equal(band$ymin, frame$lower) + expect_equal(band$ymax, frame$upper) +}) + +test_that("a continuous sweep frame is drawn against the quantile levels", { + local_quiet() + grid <- c(0, 0.05, 0.1) + res <- fit_pull_eta( + data = sim_pull_eta_continuous(), + exposure_type = "continuous", + truncation_grid = grid, + n_boot = 20 + ) + + frame <- pull_plot_data(res, type = "sweep") + + # A continuous run caps a stabilized weight rather than bounding a + # probability, so it reports no lower bound at any level and the quantile + # levels the caller swept are what the sweep can be drawn against. + expect_true(all(is.na(res@results$truncation_lower))) + expect_identical(frame$truncation, grid) + expect_identical(frame$bias, res@results$bias) + + built <- ggplot2::ggplot_build(ggplot2::autoplot(res, type = "sweep")) + points <- layer_drawing(built, "shape") + expect_equal(points$x, frame$truncation) + expect_equal(points$y, frame$bias) +}) + +test_that("pull_plot_data() and autoplot() refuse the sweep view alike", { + local_quiet() + res <- fit_pull_eta(n_boot = 10) + + pull_error <- rlang::catch_cnd(pull_plot_data(res, type = "sweep")) + plot_error <- rlang::catch_cnd(ggplot2::autoplot(res, type = "sweep")) + + expect_identical(conditionMessage(pull_error), conditionMessage(plot_error)) + expect_match(deparsed_call(pull_error), "pull_plot_data", fixed = TRUE) + expect_match(deparsed_call(plot_error), "autoplot", fixed = TRUE) + + expect_snapshot_abort( + pull_plot_data(res, type = "sweep"), + class = "positively_sweep_absent_error" + ) + expect_snapshot_abort( + ggplot2::autoplot(res, type = "sweep"), + class = "positively_sweep_absent_error" + ) +}) + +# ---- Hat values ----------------------------------------------------------- + +# A dose that tracks its one covariate, so leverage rises at the ends of the +# exposure range and the profile has structure to show. +sim_pull_hat <- function(n = 150, seed = 1) { + withr::local_seed(seed) + x1 <- stats::rnorm(n) + tibble::tibble(dose = stats::rnorm(n, mean = x1), x1 = x1) +} + +test_that("the hat-value menu mirrors autoplot's", { + local_quiet() + res <- check_hat_values(sim_pull_hat(), dose, x1, null_reps = 5) + + expect_identical(autoplot_types(res), c("null", "profile")) + expect_identical(pull_plot_data_types(hat_values_result), autoplot_types(res)) + expect_identical(pull_plot_data(res), pull_plot_data(res, type = "null")) +}) + +test_that("the null frame carries the null draws and the lines beside them", { + local_quiet() + # The null draw is seeded so the quantile the frame carries is the one the + # figure marks. + res <- withr::with_seed( + 2024, + check_hat_values(sim_pull_hat(), dose, x1, null_reps = 40) + ) + + frame <- pull_plot_data(res, type = "null") + + expect_identical( + names(frame), + c("phi", "null_quantile", "phi_hat", "conf_level") + ) + for (column in names(frame)) { + expect_type(frame[[column]], "double") + } + # `@null_dist` is the only place the null replicates live. The three scalars + # the figure marks ride as repeated columns, because a view returns one tibble + # whatever the figure does with it. + expect_identical(nrow(frame), 40L) + expect_identical(frame$phi, res@null_dist) + expect_identical(frame$null_quantile, rep(res@null_quantile, 40L)) + expect_identical(frame$phi_hat, rep(res@phi_hat, 40L)) + expect_identical(frame$conf_level, rep(res@params$conf_level, 40L)) + + plot <- ggplot2::autoplot(res, type = "null") + expect_identical(plot$data, frame) + + built <- ggplot2::ggplot_build(plot) + + bins <- layer_drawing(built, "count") + expect_identical(sum(bins$count), as.numeric(nrow(frame))) + expect_gte(min(frame$phi), min(bins$xmin)) + expect_lte(max(frame$phi), max(bins$xmax)) + + # Two reference lines, read by what draws them rather than by the order they + # were added in: the observed reading is the solid coloured one and the null + # quantile the dashed one. + lines <- layers_drawing(built, "xintercept") + expect_length(lines, 2) + observed <- Filter(function(line) all(line$colour == "firebrick"), lines) + dashed <- Filter(function(line) all(line$linetype == "dashed"), lines) + expect_length(observed, 1) + expect_length(dashed, 1) + expect_identical(observed[[1]]$xintercept, unique(frame$phi_hat)) + expect_identical(dashed[[1]]$xintercept, unique(frame$null_quantile)) +}) + +test_that("the profile frame is the fraction the profile view draws", { + local_quiet() + res <- check_hat_values(sim_pull_hat(), dose, x1, null_reps = 5) + + frame <- pull_plot_data(res, type = "profile") + + expect_identical(names(frame), c("prob", "fraction")) + expect_type(frame$prob, "double") + expect_type(frame$fraction, "double") + # One row per exposure percentile the candidates were built at, holding the + # share of that percentile's candidates that read as high leverage. The + # results carry one row per candidate, so the aggregation is the whole of what + # the view adds and it belongs in the frame. + expect_equal(frame$prob, sort(unique(res@results$prob))) + expect_equal( + frame$fraction, + as.double(tapply(res@results$high_leverage, res@results$prob, mean)) + ) + + plot <- ggplot2::autoplot(res, type = "profile") + expect_identical(plot$data, frame) + + built <- ggplot2::ggplot_build(plot) + points <- layer_drawing(built, "shape") + line <- layer_drawing(built, "flipped_aes") + expect_equal(points$x, frame$prob) + expect_equal(points$y, frame$fraction) + expect_equal(line$x, frame$prob) + expect_equal(line$y, frame$fraction) +}) + +# ---- PoRT ----------------------------------------------------------------- + +# A deterministic reading-rule anchor. One binary covariate g: 300 rows at +# g == 1 with 3 treated, a prevalence of 0.010, and 700 rows at g == 0 with 350 +# treated, a prevalence of 0.500. Nothing is drawn at random, so the reported +# subgroups are the same two on every run. +sim_pull_port <- function() { + tibble::tibble( + exposure = c(rep(1L, 3), rep(0L, 297), rep(1L, 350), rep(0L, 350)), + g = c(rep(1L, 300), rep(0L, 700)) + ) +} + +# A censoring result assembled directly, so the exposure and the censoring +# family carry thresholds that differ from each other and from wave to wave. +# What is under test is which pair of thresholds a row was judged against, and a +# run whose thresholds all resolve to one number cannot show it. +make_pull_port_seq <- function() { + results <- tibble::tibble( + time = c(1L, 1L, 2L, 2L), + type = c("exposure", "censoring", "exposure", "censoring"), + subgroup = rep("g", 4L), + description = c("g>=1", "g>=2", "g>=3", "g>=4"), + exposure_level = rep("1", 4L), + n = c(10L, 20L, 30L, 40L), + proportion = c(0.1, 0.2, 0.3, 0.4), + prevalence = c(0.01, 0.02, 0.03, 0.04), + low_support = c(TRUE, TRUE, FALSE, TRUE) + ) + port_result( + results = results, + exposure = c("a1", "a2"), + exposure_type = "binary", + n = 100L, + params = list(), + call = quote(check_port_seq()), + trees = list(), + alpha = 0.05, + beta = c(0.05, 0.1), + censoring_beta = c(0.2, 0.3), + gamma = 2 + ) +} + +test_that("the PoRT frame is the plotting frame and the thresholds behind it", { + local_quiet() + res <- check_port(sim_pull_port(), exposure, g) + + # One view, so no view to name. The accessor takes the argument the plot takes + # instead, and in the same position. + expect_identical( + names(formals(S7::method(pull_plot_data, port_result))), + c("x", "low_support_only", "...") + ) + expect_identical( + names(formals(autoplot_method(res)))[[2]], + "low_support_only" + ) + + frame <- pull_plot_data(res) + + expect_identical( + names(frame), + c( + "subgroup", + "description", + "exposure_level", + "n", + "proportion", + "prevalence", + "low_support", + "label", + "width", + "beta_lower", + "beta_upper" + ) + ) + expect_s3_class(frame$low_support, "factor") + expect_identical(levels(frame$low_support), c("FALSE", "TRUE")) + expect_type(frame$beta_lower, "double") + expect_type(frame$beta_upper, "double") + + # Duplicate rows, which arise when the same leaf is reported by more than one + # tree, are collapsed to one bar, so the frame is the reported subgroups + # rather than the results. + expect_identical(nrow(frame), 2L) + expect_equal(frame$prevalence, c(0.01, 0.5)) + expect_identical(as.character(frame$low_support), c("TRUE", "FALSE")) + + # Every row carries the pair of thresholds its prevalence was read against, + # which is what a reader needs to see why a subgroup reads as low support. + expect_identical(frame$beta_lower, rep(res@beta[[1]], nrow(frame))) + expect_identical(frame$beta_upper, rep(1 - res@beta[[1]], nrow(frame))) + + plot <- ggplot2::autoplot(res) + expect_identical(plot$data, frame) + + built <- ggplot2::ggplot_build(plot) + bars <- layer_drawing(built, "xmax") + expect_identical(nrow(bars), nrow(frame)) + drawn <- bars[order(bars$xmax), ] + expected <- frame[order(frame$prevalence), ] + expect_equal(drawn$xmax, expected$prevalence) + # The bars run horizontally, so the width proportional to subgroup size is the + # extent each bar takes in y. + expect_equal(drawn$ymax - drawn$ymin, expected$width) + + lines <- layer_drawing(built, "xintercept") + expect_setequal( + lines$xintercept, + c(unique(frame$beta_lower), unique(frame$beta_upper)) + ) +}) + +test_that("each PoRT row carries the thresholds its own facet is drawn with", { + res <- make_pull_port_seq() + + frame <- pull_plot_data(res) + reference <- port_beta_reference(res) + + expect_true(all(c("time", "type") %in% names(frame))) + expect_type(frame$beta_lower, "double") + expect_type(frame$beta_upper, "double") + # A censoring run judges the two families on different risk sets, so a row's + # thresholds follow its family as well as its wave. + expect_identical(frame$beta_lower, c(0.05, 0.2, 0.1, 0.3)) + expect_identical(frame$beta_upper, 1 - c(0.05, 0.2, 0.1, 0.3)) + + # autoplot() draws its reference lines from a frame built per facet. The two + # have to name the same pair wherever a facet has a bar in it, or a row would + # be read against one pair and drawn against another. + for (row in seq_len(nrow(frame))) { + facet <- reference[ + reference$time == frame$time[[row]] & + reference$type == frame$type[[row]], + ] + expect_setequal( + facet$xintercept, + c(frame$beta_lower[[row]], frame$beta_upper[[row]]) + ) + } +}) + +test_that("low_support_only filters the frame and the bars alike", { + local_quiet() + res <- check_port(sim_pull_port(), exposure, g) + + full <- pull_plot_data(res) + filtered <- pull_plot_data(res, low_support_only = TRUE) + + expect_identical(names(filtered), names(full)) + expect_identical(nrow(filtered), 1L) + expect_identical(as.character(filtered$low_support), "TRUE") + # The widths are read against the largest subgroup still shown rather than + # against one the filter dropped. + expect_equal(max(filtered$width), 1) + + plot <- ggplot2::autoplot(res, low_support_only = TRUE) + expect_identical(plot$data, filtered) + + built <- ggplot2::ggplot_build(plot) + expect_identical(nrow(layer_drawing(built, "xmax")), nrow(filtered)) +}) + +test_that("pull_plot_data() and autoplot() refuse a non-flag alike", { + local_quiet() + res <- check_port(sim_pull_port(), exposure, g) + + pull_error <- rlang::catch_cnd(pull_plot_data(res, low_support_only = "yes")) + plot_error <- rlang::catch_cnd( + ggplot2::autoplot(res, low_support_only = "yes") + ) + + # The argument the two entry points share is validated once behind both, so + # the reader is told the same thing either way and is shown the call they + # made. + expect_identical(conditionMessage(pull_error), conditionMessage(plot_error)) + expect_match(deparsed_call(pull_error), "pull_plot_data", fixed = TRUE) + expect_match(deparsed_call(plot_error), "autoplot", fixed = TRUE) + + expect_snapshot_abort( + pull_plot_data(res, low_support_only = "yes"), + class = "positively_type_error" + ) + expect_snapshot_abort( + ggplot2::autoplot(res, low_support_only = "yes"), + class = "positively_type_error" + ) +}) + +# ---- HDR ------------------------------------------------------------------ + +sim_pull_hdr <- function(n = 150, seed = 1) { + withr::local_seed(seed) + l <- stats::rnorm(n) + tibble::tibble(exposure = stats::rnorm(n, mean = l), l = l) +} + +test_that("the point HDR frame is the results the curve draws", { + local_quiet() + res <- check_hdr(sim_pull_hdr(), exposure, l, values = c(-1, 0, 1)) + + # One view, so there is no view to name and neither signature offers a menu. + expect_null(formals(S7::method(pull_plot_data, hdr_result))$type) + expect_null(formals(autoplot_method(res))$type) + + frame <- pull_plot_data(res) + + # A point result has nothing to key the curves on, so the frame is the results + # as they stand. + expect_identical(frame, res@results) + expect_identical(names(frame), c("value", "nonoverlap")) + + plot <- ggplot2::autoplot(res) + expect_identical(plot$data, frame) + + built <- ggplot2::ggplot_build(plot) + points <- layer_drawing(built, "shape") + expect_equal(points$x, frame$value) + expect_equal(points$y, frame$nonoverlap) +}) + +test_that("the sequential HDR frame keys the curves on the time point", { + local_quiet() + data <- dgp_longitudinal(n = 200, seed = 5) + res <- check_hdr_seq( + data, + c(a1, a2, a3), + list(l0, l1, l2), + values = c(-1, 0, 1) + ) + + frame <- pull_plot_data(res) + + # The wave separates the curves, so it is a factor here where the results + # carry it as the integer index. Left numeric it would colour the curves off a + # continuous scale, which reads a wave as a quantity. + expect_identical(names(frame), names(res@results)) + expect_s3_class(frame$time, "factor") + expect_identical(levels(frame$time), c("1", "2", "3")) + expect_identical(as.integer(as.character(frame$time)), res@results$time) + expect_identical(frame$value, res@results$value) + expect_identical(frame$nonoverlap, res@results$nonoverlap) + + plot <- ggplot2::autoplot(res) + expect_identical(plot$data, frame) + + built <- ggplot2::ggplot_build(plot) + points <- layer_drawing(built, "shape") + # One curve per wave, together drawing every row of the frame. + expect_length(unique(points$group), nlevels(frame$time)) + expect_identical(nrow(points), nrow(frame)) + expect_equal(sort(points$y), sort(frame$nonoverlap)) +}) + +# ---- Extrapolation -------------------------------------------------------- + +# Two exposure groups drawn from one Gaussian covariate cloud, so both groups +# cover each other and the hull test has a mix of inside and outside rows to +# report. +sim_pull_extrapolation <- function(n = 60, p = 3, seed = 1) { + withr::local_seed(seed) + covariates <- matrix(stats::rnorm(n * p), nrow = n, ncol = p) + colnames(covariates) <- paste0("x", seq_len(p)) + out <- tibble::as_tibble(as.data.frame(covariates)) + n0 <- n %/% 2L + out$exposure <- c(rep(0L, n0), rep(1L, n - n0)) + out[c("exposure", paste0("x", seq_len(p)))] +} + +test_that("the extrapolation menu mirrors autoplot's", { + local_quiet() + res <- check_extrapolation( + sim_pull_extrapolation(), + exposure, + tidyselect::starts_with("x"), + hull = FALSE + ) + + expect_identical(autoplot_types(res), c("distribution", "hull")) + expect_identical( + pull_plot_data_types(extrapolation_result), + autoplot_types(res) + ) + expect_identical( + pull_plot_data(res), + pull_plot_data(res, type = "distribution") + ) +}) + +test_that("the distribution frame groups the results by exposure", { + local_quiet() + res <- check_extrapolation( + sim_pull_extrapolation(), + exposure, + tidyselect::starts_with("x"), + hull = FALSE + ) + + frame <- pull_plot_data(res, type = "distribution") + + # The exposure panels the view rather than measuring anything, so it is a + # factor here where the results carry the raw level. + expect_identical(names(frame), names(res@results)) + expect_s3_class(frame$exposure, "factor") + expect_identical(levels(frame$exposure), c("0", "1")) + expect_identical(frame$frac_nearby, res@results$frac_nearby) + + plot <- ggplot2::autoplot(res, type = "distribution") + expect_identical(plot$data, frame) + + built <- ggplot2::ggplot_build(plot) + # One panel per exposure group, in level order, with every row of the frame + # binned inside the panel of its own group. + expect_identical( + as.character(built$layout$layout$exposure), + levels(frame$exposure) + ) + bins <- layer_drawing(built, "count") + expect_identical( + unname(vapply(split(bins$count, bins$PANEL), sum, numeric(1))), + as.numeric(table(frame$exposure)) + ) +}) + +test_that("the hull frame is the membership the bars draw", { + local_quiet() + skip_if_not_installed("lpSolve") + res <- check_extrapolation( + sim_pull_extrapolation(), + exposure, + tidyselect::starts_with("x"), + hull = TRUE + ) + + frame <- pull_plot_data(res, type = "hull") + + # The view stacks the bars on a reading of `in_hull` rather than on the + # logical itself, so the reading is a column of the frame and reads the way + # the legend does. + expect_identical(names(frame), c(names(res@results), "membership")) + expect_s3_class(frame$exposure, "factor") + expect_s3_class(frame$membership, "factor") + expect_identical(levels(frame$membership), c("Inside hull", "Outside hull")) + expect_identical(frame$membership == "Inside hull", res@results$in_hull) + + plot <- ggplot2::autoplot(res, type = "hull") + expect_identical(plot$data, frame) + + built <- ggplot2::ggplot_build(plot) + bars <- layer_drawing(built, "count") + # One bar per exposure group per membership, counting the frame's rows and no + # others. + expect_identical(nrow(bars), 4L) + expect_identical(sum(bars$count), as.numeric(nrow(frame))) + expect_identical( + unname(vapply(split(bars$count, bars$x), sum, numeric(1))), + as.numeric(table(frame$exposure)) + ) + expect_setequal( + bars$count, + as.numeric(table(frame$exposure, frame$membership)) + ) + + # The key names the memberships the frame carries, in the frame's own order. + guide <- ggplot2::get_guide_data(plot, "fill") + expect_identical(as.character(guide$.label), levels(frame$membership)) +}) + +test_that("pull_plot_data() and autoplot() refuse the hull view alike", { + local_quiet() + res <- check_extrapolation( + sim_pull_extrapolation(), + exposure, + tidyselect::starts_with("x"), + hull = FALSE + ) + + pull_error <- rlang::catch_cnd(pull_plot_data(res, type = "hull")) + plot_error <- rlang::catch_cnd(ggplot2::autoplot(res, type = "hull")) + + expect_identical(conditionMessage(pull_error), conditionMessage(plot_error)) + expect_match(deparsed_call(pull_error), "pull_plot_data", fixed = TRUE) + expect_match(deparsed_call(plot_error), "autoplot", fixed = TRUE) + + expect_snapshot_abort( + pull_plot_data(res, type = "hull"), + class = "positively_hull_absent_error" + ) + expect_snapshot_abort( + ggplot2::autoplot(res, type = "hull"), + class = "positively_hull_absent_error" + ) +}) + +# ---- Positivity check ----------------------------------------------------- + +sim_pull_check <- function(n = 200, seed = 1) { + withr::local_seed(seed) + x1 <- stats::rnorm(n) + tibble::tibble( + exposure = stats::rbinom(n, 1L, stats::plogis(3 * x1)), + x1 = x1 + ) +} + +test_that("a positivity check delegates to the diagnostic it is named", { + local_quiet() + check <- check_positivity( + sim_pull_check(), + exposure, + x1, + diagnostics = c("port", "edp") + ) + + # A container holds children rather than a frame of its own, so naming one is + # the whole of what this does: the frame is the child's, unchanged. + expect_identical(pull_plot_data(check, "port"), pull_plot_data(check$port)) + expect_identical(pull_plot_data(check, "edp"), pull_plot_data(check$edp)) + + # A view name reaches the child through the dots, so a child's menu is + # reachable without extracting it first. + expect_identical( + pull_plot_data(check, "edp", type = "ecdf"), + pull_plot_data(check$edp, type = "ecdf") + ) +}) + +test_that("a positivity check refuses to guess which diagnostic to pull", { + local_quiet() + check <- check_positivity( + sim_pull_check(), + exposure, + x1, + diagnostics = c("port", "edp") + ) + + # A container holds one frame per child and this returns one tibble, so + # without a name there is nothing to return. autoplot() has a panel to fall + # back on; there is no such thing for a frame. + missing_error <- rlang::catch_cnd(pull_plot_data(check)) + expect_s3_class(missing_error, "positively_diagnostic_error") + for (name in names(check)) { + expect_match(conditionMessage(missing_error), name, fixed = TRUE) + } + + # autoplot() spells the panel as a `NULL` diagnostic, so a reader carrying the + # argument over lands on the same refusal rather than on an internal failure. + null_error <- rlang::catch_cnd(pull_plot_data(check, NULL)) + expect_s3_class(null_error, "positively_diagnostic_error") + expect_identical( + conditionMessage(null_error), + conditionMessage(missing_error) + ) + + # extract_named_check() is what `$` and `[[` refuse an unknown name with, so + # every way of naming a child spells the refusal the same, and the reader is + # shown the call they made. + unknown_error <- rlang::catch_cnd(pull_plot_data(check, "bogus")) + extract_error <- rlang::catch_cnd(check[["bogus"]]) + expect_s3_class(unknown_error, "positively_diagnostic_error") + expect_identical( + conditionMessage(unknown_error), + conditionMessage(extract_error) + ) + expect_match(deparsed_call(unknown_error), "pull_plot_data", fixed = TRUE) + + expect_snapshot_abort( + pull_plot_data(check), + class = "positively_diagnostic_error" + ) + expect_snapshot_abort( + pull_plot_data(check, "bogus"), + class = "positively_diagnostic_error" + ) +}) From 2d7ecbd5d56634db4362cc3ccad21fb9968dd4f6 Mon Sep 17 00:00:00 2001 From: Malcolm Barrett Date: Fri, 21 Aug 2026 22:36:44 +0200 Subject: [PATCH 5/7] draw the diagnostic histograms in the stock ggplot2 appearance --- R/plot-density-ratios.R | 2 +- R/plot-edp.R | 2 +- R/plot-eta-bias.R | 2 +- R/plot-extrapolation.R | 7 +- R/plot-hat-values.R | 2 +- .../density-ratios-point-histogram.svg | 60 +-- .../edp-histogram-by-intervention-value.svg | 120 +++--- .../edp-histogram-shift-intervention.svg | 120 +++--- .../eta-bias-bootstrap-distribution.svg | 60 +-- .../eta-bias-categorical-bootstrap-grid.svg | 360 +++++++++--------- ...extrapolation-frac-nearby-distribution.svg | 120 +++--- .../hat-values-null-distribution.svg | 16 +- .../positivity-check-panel.svg | 120 +++--- tests/testthat/helper-plot.R | 24 ++ tests/testthat/test-check-density-ratios.R | 11 + tests/testthat/test-check-edp.R | 22 ++ tests/testthat/test-check-eta-bias.R | 34 ++ tests/testthat/test-check-extrapolation.R | 18 + tests/testthat/test-check-hat-values.R | 15 + tests/testthat/test-check-port.R | 15 + 20 files changed, 632 insertions(+), 498 deletions(-) create mode 100644 tests/testthat/helper-plot.R diff --git a/R/plot-density-ratios.R b/R/plot-density-ratios.R index b31fcb6..c387e30 100644 --- a/R/plot-density-ratios.R +++ b/R/plot-density-ratios.R @@ -51,7 +51,7 @@ method(autoplot, density_ratios_result) <- function( autoplot_density_ratios_histogram <- function(object) { plot_data <- pull_plot_data(object, type = "distribution") ggplot2::ggplot(plot_data, ggplot2::aes(x = .data$ratio)) + - ggplot2::geom_histogram(bins = 30, fill = "grey70", color = "white") + + ggplot2::geom_histogram(bins = 30) + ggplot2::labs( x = "Density ratio", y = "Observations", diff --git a/R/plot-edp.R b/R/plot-edp.R index 0e5540b..adef9db 100644 --- a/R/plot-edp.R +++ b/R/plot-edp.R @@ -144,7 +144,7 @@ autoplot_edp_histogram <- function(object) { measure <- edp_primary_column(object@results) plot_data <- pull_plot_data(object, type = "histogram") ggplot2::ggplot(plot_data, ggplot2::aes(x = .data[[measure]])) + - ggplot2::geom_histogram(bins = 30, fill = "grey70", color = "white") + + ggplot2::geom_histogram(bins = 30) + ggplot2::facet_wrap( ggplot2::vars(.data$intervention), labeller = ggplot2::label_both diff --git a/R/plot-eta-bias.R b/R/plot-eta-bias.R index afc00d4..b82fee3 100644 --- a/R/plot-eta-bias.R +++ b/R/plot-eta-bias.R @@ -146,7 +146,7 @@ autoplot_eta_bias_bootstrap <- function(object) { ) ggplot2::ggplot(boot_data, ggplot2::aes(x = .data$estimate)) + - ggplot2::geom_histogram(bins = 30, fill = "grey70", color = "white") + + ggplot2::geom_histogram(bins = 30) + ggplot2::geom_vline( ggplot2::aes(xintercept = .data$truth), data = truth_data, diff --git a/R/plot-extrapolation.R b/R/plot-extrapolation.R index 4fbf3e5..2b643e7 100644 --- a/R/plot-extrapolation.R +++ b/R/plot-extrapolation.R @@ -56,12 +56,7 @@ autoplot_extrapolation_distribution <- function(object) { plot_data, ggplot2::aes(x = .data$frac_nearby) ) + - ggplot2::geom_histogram( - bins = 30, - boundary = 0, - fill = "grey70", - color = "white" - ) + + ggplot2::geom_histogram(bins = 30, boundary = 0) + ggplot2::facet_wrap( ggplot2::vars(.data$exposure), ncol = 1, diff --git a/R/plot-hat-values.R b/R/plot-hat-values.R index 9adf680..9ebc0cd 100644 --- a/R/plot-hat-values.R +++ b/R/plot-hat-values.R @@ -58,7 +58,7 @@ autoplot_hat_values_null <- function(object) { # capped at 30 so large null distributions are not oversmoothed. bins <- min(30L, max(1L, ceiling(sqrt(nrow(null_data))))) ggplot2::ggplot(null_data, ggplot2::aes(x = .data$phi)) + - ggplot2::geom_histogram(bins = bins, fill = "grey70", color = "white") + + ggplot2::geom_histogram(bins = bins) + ggplot2::geom_vline( xintercept = object@null_quantile, linetype = "dashed" diff --git a/tests/testthat/_snaps/check-density-ratios/density-ratios-point-histogram.svg b/tests/testthat/_snaps/check-density-ratios/density-ratios-point-histogram.svg index a5cfdb3..fe4b0ba 100644 --- a/tests/testthat/_snaps/check-density-ratios/density-ratios-point-histogram.svg +++ b/tests/testthat/_snaps/check-density-ratios/density-ratios-point-histogram.svg @@ -27,36 +27,36 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/testthat/_snaps/check-edp/edp-histogram-by-intervention-value.svg b/tests/testthat/_snaps/check-edp/edp-histogram-by-intervention-value.svg index edd7ce4..9db74e6 100644 --- a/tests/testthat/_snaps/check-edp/edp-histogram-by-intervention-value.svg +++ b/tests/testthat/_snaps/check-edp/edp-histogram-by-intervention-value.svg @@ -27,36 +27,36 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -68,36 +68,36 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/testthat/_snaps/check-edp/edp-histogram-shift-intervention.svg b/tests/testthat/_snaps/check-edp/edp-histogram-shift-intervention.svg index b5812d1..e3bde3d 100644 --- a/tests/testthat/_snaps/check-edp/edp-histogram-shift-intervention.svg +++ b/tests/testthat/_snaps/check-edp/edp-histogram-shift-intervention.svg @@ -27,36 +27,36 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -68,36 +68,36 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/testthat/_snaps/check-eta-bias/eta-bias-bootstrap-distribution.svg b/tests/testthat/_snaps/check-eta-bias/eta-bias-bootstrap-distribution.svg index f957643..aa55e4d 100644 --- a/tests/testthat/_snaps/check-eta-bias/eta-bias-bootstrap-distribution.svg +++ b/tests/testthat/_snaps/check-eta-bias/eta-bias-bootstrap-distribution.svg @@ -27,36 +27,36 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/testthat/_snaps/check-eta-bias/eta-bias-categorical-bootstrap-grid.svg b/tests/testthat/_snaps/check-eta-bias/eta-bias-categorical-bootstrap-grid.svg index b678d14..505a604 100644 --- a/tests/testthat/_snaps/check-eta-bias/eta-bias-categorical-bootstrap-grid.svg +++ b/tests/testthat/_snaps/check-eta-bias/eta-bias-categorical-bootstrap-grid.svg @@ -27,36 +27,36 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -69,36 +69,36 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -111,36 +111,36 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -153,36 +153,36 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -195,36 +195,36 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -237,36 +237,36 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/testthat/_snaps/check-extrapolation/extrapolation-frac-nearby-distribution.svg b/tests/testthat/_snaps/check-extrapolation/extrapolation-frac-nearby-distribution.svg index d29605c..5e1ab16 100644 --- a/tests/testthat/_snaps/check-extrapolation/extrapolation-frac-nearby-distribution.svg +++ b/tests/testthat/_snaps/check-extrapolation/extrapolation-frac-nearby-distribution.svg @@ -27,36 +27,36 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -68,36 +68,36 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/testthat/_snaps/check-hat-values/hat-values-null-distribution.svg b/tests/testthat/_snaps/check-hat-values/hat-values-null-distribution.svg index 45c57c5..bb8bc01 100644 --- a/tests/testthat/_snaps/check-hat-values/hat-values-null-distribution.svg +++ b/tests/testthat/_snaps/check-hat-values/hat-values-null-distribution.svg @@ -27,14 +27,14 @@ - - - - - - - - + + + + + + + + diff --git a/tests/testthat/_snaps/check-positivity/positivity-check-panel.svg b/tests/testthat/_snaps/check-positivity/positivity-check-panel.svg index 4bd3b37..2569478 100644 --- a/tests/testthat/_snaps/check-positivity/positivity-check-panel.svg +++ b/tests/testthat/_snaps/check-positivity/positivity-check-panel.svg @@ -317,36 +317,36 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -375,36 +375,36 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/testthat/helper-plot.R b/tests/testthat/helper-plot.R new file mode 100644 index 0000000..ddf537d --- /dev/null +++ b/tests/testthat/helper-plot.R @@ -0,0 +1,24 @@ +# The layer of a figure that draws a given aesthetic, found by what it carries +# rather than by the position it was added in. +built_layer <- function(plot, aesthetic) { + built <- ggplot2::ggplot_build(plot) + Filter(function(layer) aesthetic %in% names(layer), built$data)[[1]] +} + +# A histogram is the layer carrying a count. +histogram_bars <- function(plot) { + built_layer(plot, "count") +} + +# The fill and the colour a stock geom_histogram draws with, read off a +# reference figure. ggplot2 4.0 resolves a geom's defaults from the theme, so +# `default_aes` holds the rule rather than the value it settles on, and +# `get_geom_defaults()` postdates the declared ggplot2 floor of 3.5.0. Building +# a reference resolves the defaults through the pipeline the views under test go +# through, so what it returns is the stock appearance whatever version draws it. +stock_histogram_aes <- function() { + reference <- ggplot2::ggplot(mapping = ggplot2::aes(x = c(0, 1, 2, 3))) + + ggplot2::geom_histogram(bins = 2) + bars <- histogram_bars(reference) + list(fill = unique(bars$fill), colour = unique(bars$colour)) +} diff --git a/tests/testthat/test-check-density-ratios.R b/tests/testthat/test-check-density-ratios.R index d609156..9370a79 100644 --- a/tests/testthat/test-check-density-ratios.R +++ b/tests/testthat/test-check-density-ratios.R @@ -582,6 +582,17 @@ test_that("autoplot() rejects an unknown type as a classed error", { ) }) +test_that("the point histogram draws in the stock geom_histogram appearance", { + # The bars stand for counts and nothing else, so they carry the fill and the + # outline ggplot2 gives a histogram rather than a pair of their own. + res <- check_density_ratios(gen_lognormal_ratios(200, k = 1)) + + bars <- histogram_bars(ggplot2::autoplot(res, type = "distribution")) + stock <- stock_histogram_aes() + expect_identical(unique(bars$fill), stock$fill) + expect_identical(unique(bars$colour), stock$colour) +}) + test_that("density-ratio autoplot views render as expected", { announce_doppelganger( "Density ratios point histogram", diff --git a/tests/testthat/test-check-edp.R b/tests/testthat/test-check-edp.R index 1a40be0..de2df4e 100644 --- a/tests/testthat/test-check-edp.R +++ b/tests/testthat/test-check-edp.R @@ -1661,6 +1661,10 @@ test_that("the scatter view separates infinite ideal weights into their own laye # The infinite rows carry a single fixed shape distinct from the finite layer. expect_length(unique(infinite_layer$shape), 1) expect_false(unique(infinite_layer$shape) %in% unique(finite_layer$shape)) + # The crosses sit off the gradient the finite rows are read on, so their + # colour is fixed and is not a value that gradient can take. + expect_identical(unique(infinite_layer$colour), "#B2182B") + expect_false(unique(infinite_layer$colour) %in% unique(finite_layer$colour)) }) test_that("the scatter view builds without a colour label when every weight is infinite", { @@ -1714,6 +1718,24 @@ test_that("the histogram view facets on the intervention label", { expect_setequal(as.character(layout$intervention), c("down", "0")) }) +test_that("the histogram view draws in the stock geom_histogram appearance", { + # The bars stand for counts and nothing else, so they carry the fill and the + # outline ggplot2 gives a histogram rather than a pair of their own. + data <- sim_edp_gaussian(60) + res <- check_edp( + data, + exposure, + x1, + values = c(0, 1), + exposure_type = "continuous" + ) + + bars <- histogram_bars(ggplot2::autoplot(res, type = "histogram")) + stock <- stock_histogram_aes() + expect_identical(unique(bars$fill), stock$fill) + expect_identical(unique(bars$colour), stock$colour) +}) + test_that("the ecdf view colors by the intervention label", { data <- sim_edp_gaussian(60) res <- check_edp( diff --git a/tests/testthat/test-check-eta-bias.R b/tests/testthat/test-check-eta-bias.R index 9dc7198..89b1032 100644 --- a/tests/testthat/test-check-eta-bias.R +++ b/tests/testthat/test-check-eta-bias.R @@ -1609,6 +1609,40 @@ test_that("autoplot() rejects an unknown type as a classed error", { ) }) +test_that("the bootstrap view draws in the stock geom_histogram appearance", { + # The bars stand for counts and nothing else, so they carry the fill and the + # outline ggplot2 gives a histogram rather than a pair of their own. + res <- make_eta_bias_result() + + bars <- histogram_bars(ggplot2::autoplot(res, type = "bootstrap")) + stock <- stock_histogram_aes() + expect_identical(unique(bars$fill), stock$fill) + expect_identical(unique(bars$colour), stock$colour) +}) + +test_that("the bootstrap view marks its truth with a dashed line", { + # The line is a reference rather than a reading, and the dash is what tells a + # reader so. + res <- make_eta_bias_result() + + plot <- ggplot2::autoplot(res, type = "bootstrap") + expect_identical(unique(built_layer(plot, "xintercept")$linetype), "dashed") +}) + +test_that("a one-term sweep keeps its grey band and its dashed zero line", { + # One term needs no key to tell terms apart, so the band is neutral rather + # than a fill a legend would have to name, and the line at zero is a reference + # rather than a reading. + res <- make_eta_bias_result( + results = eta_bias_stub_results(truncation_lower = c(0, 0.05, 0.1)), + boot_estimates = rep(list(c(0.3, 0.5)), 3) + ) + + plot <- ggplot2::autoplot(res, type = "sweep") + expect_identical(unique(built_layer(plot, "ymax")$fill), "grey80") + expect_identical(unique(built_layer(plot, "yintercept")$linetype), "dashed") +}) + test_that("ETA bias autoplot views render as expected", { local_quiet() announce_doppelganger( diff --git a/tests/testthat/test-check-extrapolation.R b/tests/testthat/test-check-extrapolation.R index b08f4a8..85c5875 100644 --- a/tests/testthat/test-check-extrapolation.R +++ b/tests/testthat/test-check-extrapolation.R @@ -1101,6 +1101,24 @@ test_that("plot() draws the view and returns the result invisibly", { expect_identical(plot(res), res) }) +test_that("the distribution view draws in the stock geom_histogram appearance", { + local_quiet() + # The bars stand for counts and nothing else, so they carry the fill and the + # outline ggplot2 gives a histogram rather than a pair of their own. + data <- sim_extrap_gaussian(200, p = 4, sep = 0, seed = 1) + res <- check_extrapolation( + data, + exposure, + tidyselect::starts_with("x"), + hull = FALSE + ) + + bars <- histogram_bars(ggplot2::autoplot(res, type = "distribution")) + stock <- stock_histogram_aes() + expect_identical(unique(bars$fill), stock$fill) + expect_identical(unique(bars$colour), stock$colour) +}) + test_that("the extrapolation distribution view renders as expected", { local_quiet() # The distribution view reads only `frac_nearby` and `exposure`, so it renders diff --git a/tests/testthat/test-check-hat-values.R b/tests/testthat/test-check-hat-values.R index 6e92c5f..539b771 100644 --- a/tests/testthat/test-check-hat-values.R +++ b/tests/testthat/test-check-hat-values.R @@ -1078,6 +1078,21 @@ test_that("plot() draws the view and returns the result invisibly", { expect_identical(plot(res), res) }) +test_that("the null view draws in the stock geom_histogram appearance", { + local_quiet() + # The bars stand for counts and nothing else, so they carry the fill and the + # outline ggplot2 gives a histogram rather than a pair of their own. The two + # lines the view draws over them are the only part of it whose appearance + # carries meaning. + data <- sim_hat_linear(150, beta = 1, seed = 1) + res <- check_hat_values(data, dose, x1, null_reps = 20) + + bars <- histogram_bars(ggplot2::autoplot(res, type = "null")) + stock <- stock_histogram_aes() + expect_identical(unique(bars$fill), stock$fill) + expect_identical(unique(bars$colour), stock$colour) +}) + test_that("hat-value autoplot views render as expected", { local_quiet() announce_doppelganger( diff --git a/tests/testthat/test-check-port.R b/tests/testthat/test-check-port.R index 90744ef..f55c8b4 100644 --- a/tests/testthat/test-check-port.R +++ b/tests/testthat/test-check-port.R @@ -863,6 +863,21 @@ test_that("PoRT autoplot renders as expected", { ) }) +test_that("the bars read low support by fill and the thresholds by a dash", { + local_quiet() + # The fill is the whole of what marks a subgroup as low support, so its two + # levels are drawn in colours a reader can tell apart. The threshold lines are + # references rather than readings, and the dash is what says so. + res <- check_port(port_anchor_data(), exposure, g) + + plot <- ggplot2::autoplot(res) + guide <- ggplot2::get_guide_data(plot, "fill") + expect_identical(as.character(guide$.label), c("FALSE", "TRUE")) + expect_identical(guide$fill, c("grey70", "#B2182B")) + + expect_identical(unique(built_layer(plot, "xintercept")$linetype), "dashed") +}) + test_that("duplicate subgroups draw a bar of their own prevalence, not a sum", { local_quiet() res <- check_port(port_anchor_data(), exposure, g) From ec45a73b2ebdfa6e13369ef92802feec7d28c546 Mon Sep 17 00:00:00 2001 From: Malcolm Barrett Date: Fri, 21 Aug 2026 22:57:00 +0200 Subject: [PATCH 6/7] document function interventions and the EDP plot views --- R/check-edp.R | 56 ++++++++++++++++++++++++------ R/plot-edp.R | 7 ++-- R/pull-plot-data.R | 13 +++---- man/autoplot.edp_result.Rd | 7 ++-- man/check_edp.Rd | 56 ++++++++++++++++++++++++------ man/pull_plot_data.Rd | 13 +++---- vignettes/continuous-exposures.Rmd | 41 ++++++++++++++++++++++ vignettes/positively.Rmd | 11 +++--- 8 files changed, 162 insertions(+), 42 deletions(-) diff --git a/R/check-edp.R b/R/check-edp.R index e49f3de..2503f77 100644 --- a/R/check-edp.R +++ b/R/check-edp.R @@ -49,6 +49,19 @@ edp_result <- new_class( #' dimensions; factor and character covariates, together with binary and #' categorical exposures, are categorical dimensions. #' +#' An intervention need not be a single value every observation is moved to. A +#' function supplied through `values` is a rule \eqn{d} applied to the observed +#' exposure, so observation \eqn{i} is moved to \eqn{d(a_i)} and the +#' intervened-on point is \eqn{o^*_i = (l_i, d(a_i))}. Lowering every observed +#' dose by five units is `function(.x) .x - 5`. EDP is then read per natural +#' observation under the shift rather than at one common target; the algorithm +#' is unchanged, and only the intervened-on points come from the rule. Because +#' those points differ across observations, the exposure factor of the kernel +#' compares every observed exposure against a target that moves with the +#' observation rather than against a single \eqn{a^*}. A run therefore holds one +#' \eqn{n \times n} matrix of exposure-kernel values per function intervention, +#' materialized before any summing and held for the duration of the run. +#' #' The half-distances follow the paper's rule of thumb: `1 * sd` for each #' numeric covariate and `0.5 * sd` for a continuous exposure. These are marginal #' scales, so a mid-support gap in a strongly multimodal exposure can be masked @@ -82,17 +95,27 @@ edp_result <- new_class( #' estimator variant, selected with tidyselect. Defaults to `.covariates`. #' @param .treatment_covariates The treatment-model covariate columns for the #' estimator variant, selected with tidyselect. Defaults to `.covariates`. -#' @param values The interventions. `NULL` (the default) uses every level of a -#' factor exposure and every observed value of a binary or character exposure, -#' and the deciles of the observed exposure for continuous exposures. A vector -#' holds one intervention value \eqn{a^*} per element. A list holds the same -#' single values and, for a continuous exposure, functions of the observed -#' exposure such as `list(3, 5, lower = function(.x) .x - 5)`. A function is -#' called once on the whole exposure vector and returns one intervened -#' exposure \eqn{d(a_i)} per observation, or a single value used for every -#' observation. Each function needs a name, which labels the intervention; an -#' unnamed value is labeled by the value itself, and every label in a list -#' must be distinct. +#' @param values The interventions, one per element. `NULL` (the default) uses +#' every level of a factor exposure and every observed value of a binary or +#' character exposure, and the deciles of the observed exposure for continuous +#' exposures. +#' +#' A vector holds one intervention value \eqn{a^*} per element, the value every +#' observation is moved to. A list holds the same single values and, for a +#' continuous exposure, functions of the observed exposure: +#' `list(3, 5, lower = function(.x) .x - 5)`. A function is called once on the +#' whole exposure vector and returns one intervened exposure \eqn{d(a_i)} per +#' observation, or a single value used for every observation. Functions need a +#' continuous exposure; one supplied for a binary or categorical exposure is an +#' error. +#' +#' Every intervention carries a label, which names it in the results and along +#' a plot's axis. In a vector, every element is labeled by `as.character()` of +#' the value itself and any names the vector carries are ignored. In a list, a +#' named element is labeled by its name, an unnamed single value by +#' `as.character()` of the value itself, and an unnamed function is an error, +#' since it has nothing else to be called by. Two elements of a list cannot +#' share a label. #' @param variant One of `"data"` (report a single `edp` per point, the default) #' or `"estimator"` (report `edp_outcome`, `edp_treatment`, and #' `ideal_weight`). @@ -142,6 +165,17 @@ edp_result <- new_class( #' result <- check_edp(df, dose, x1, values = c(0, 1), exposure_type = "continuous") #' result #' +#' # An intervention can also be a function of the observed exposure, read as a +#' # shift away from the dose each observation received +#' shifted <- check_edp( +#' df, +#' dose, +#' x1, +#' values = list(0, 1, lower = function(.x) .x - 1), +#' exposure_type = "continuous" +#' ) +#' tidy(shifted) +#' #' @export check_edp <- function( .data, diff --git a/R/plot-edp.R b/R/plot-edp.R index adef9db..5a0a590 100644 --- a/R/plot-edp.R +++ b/R/plot-edp.R @@ -11,8 +11,11 @@ #' faceted by intervention, `type = "ecdf"` shows its empirical cumulative #' distribution colored by intervention, and `type = "density"` draws one #' density ridgeline per intervention, which needs the \pkg{ggridges} package. -#' For the estimator variant, `type = "scatter"` plots `edp_outcome` against -#' `edp_treatment` colored by `ideal_weight`; it aborts for the data variant. +#' Each of those three draws one measure: `edp` for the data variant and +#' `edp_outcome` for the estimator variant, the measure carrying the exposure +#' dimension and so the one an intervention moves. `type = "scatter"` plots +#' `edp_outcome` against `edp_treatment` colored by `ideal_weight`; it needs +#' the estimator variant and aborts for the data variant. #' [pull_plot_data()] returns the tibble any of these views draws. #' #' @param object An `edp_result` from [check_edp()]. diff --git a/R/pull-plot-data.R b/R/pull-plot-data.R index 97621a4..02178bc 100644 --- a/R/pull-plot-data.R +++ b/R/pull-plot-data.R @@ -25,9 +25,9 @@ #' variant. `ymin` and `ymax` are the fifth and the ninety-fifth percentile, #' `lower`, `middle`, and `upper` are the quartiles, `n` is the number of #' observations summarized, and `outliers` is a list-column of the values -#' lying outside the whiskers, which the view draws as points. The estimator variant covers -#' `edp_outcome` and `edp_treatment`; `ideal_weight` shares no scale with -#' them and is left to the scatter view. +#' lying outside the whiskers, which the view draws as points. The estimator +#' variant covers `edp_outcome` and `edp_treatment`; `ideal_weight` shares no +#' scale with them and is left to the scatter view. #' * `"histogram"`, `"ecdf"`, and `"density"` return the results at their own #' grain, one row per observation and intervention, with `intervention` as a #' factor whose levels follow the order the interventions were given in. The @@ -52,8 +52,9 @@ #' #' * `"bootstrap"` returns one row per bootstrap draw, keyed on the estimand #' `term`, on the truncation `level`, and on the `label` its facet strip -#' reads. Each draw carries the `truth` its term was aimed at alongside the -#' `estimate`, which the view draws as one line per term. +#' reads. The view bins the `estimate`. Each draw also carries the `truth` its +#' term was aimed at, which the view draws as a reference line, one per term +#' rather than one per draw, since a term's truth is the same in every draw. #' * `"sweep"` returns one row per term and truncation level, holding the #' `truncation` the sweep is drawn against, the `bias`, and the `lower` and #' `upper` ends of the band two Monte Carlo standard errors either side of @@ -92,7 +93,7 @@ #' own, so name the diagnostic to pull: `pull_plot_data(check, "port")`. The #' frame is that child's, unchanged, and a `type` passed alongside reaches it. #' -#' @param x A [positivity_diagnostic]. +#' @param x A [positivity_diagnostic] or a [positivity_check]. #' @param ... Passed to methods. A class drawing more than one view takes #' `type`, the name of the view whose data to return, matched against that #' class's [autoplot()] menu. A [check_port()] result takes diff --git a/man/autoplot.edp_result.Rd b/man/autoplot.edp_result.Rd index 27fc339..20e44a9 100644 --- a/man/autoplot.edp_result.Rd +++ b/man/autoplot.edp_result.Rd @@ -25,8 +25,11 @@ facet per intervention. \code{type = "histogram"} shows the distribution of EDP faceted by intervention, \code{type = "ecdf"} shows its empirical cumulative distribution colored by intervention, and \code{type = "density"} draws one density ridgeline per intervention, which needs the \pkg{ggridges} package. -For the estimator variant, \code{type = "scatter"} plots \code{edp_outcome} against -\code{edp_treatment} colored by \code{ideal_weight}; it aborts for the data variant. +Each of those three draws one measure: \code{edp} for the data variant and +\code{edp_outcome} for the estimator variant, the measure carrying the exposure +dimension and so the one an intervention moves. \code{type = "scatter"} plots +\code{edp_outcome} against \code{edp_treatment} colored by \code{ideal_weight}; it needs +the estimator variant and aborts for the data variant. \code{\link[=pull_plot_data]{pull_plot_data()}} returns the tibble any of these views draws. } \examples{ diff --git a/man/check_edp.Rd b/man/check_edp.Rd index 8b69aa1..b96a235 100644 --- a/man/check_edp.Rd +++ b/man/check_edp.Rd @@ -32,17 +32,27 @@ estimator variant, selected with tidyselect. Defaults to \code{.covariates}.} \item{.treatment_covariates}{The treatment-model covariate columns for the estimator variant, selected with tidyselect. Defaults to \code{.covariates}.} -\item{values}{The interventions. \code{NULL} (the default) uses every level of a -factor exposure and every observed value of a binary or character exposure, -and the deciles of the observed exposure for continuous exposures. A vector -holds one intervention value \eqn{a^*} per element. A list holds the same -single values and, for a continuous exposure, functions of the observed -exposure such as \code{list(3, 5, lower = function(.x) .x - 5)}. A function is -called once on the whole exposure vector and returns one intervened -exposure \eqn{d(a_i)} per observation, or a single value used for every -observation. Each function needs a name, which labels the intervention; an -unnamed value is labeled by the value itself, and every label in a list -must be distinct.} +\item{values}{The interventions, one per element. \code{NULL} (the default) uses +every level of a factor exposure and every observed value of a binary or +character exposure, and the deciles of the observed exposure for continuous +exposures. + +A vector holds one intervention value \eqn{a^*} per element, the value every +observation is moved to. A list holds the same single values and, for a +continuous exposure, functions of the observed exposure: +\code{list(3, 5, lower = function(.x) .x - 5)}. A function is called once on the +whole exposure vector and returns one intervened exposure \eqn{d(a_i)} per +observation, or a single value used for every observation. Functions need a +continuous exposure; one supplied for a binary or categorical exposure is an +error. + +Every intervention carries a label, which names it in the results and along +a plot's axis. In a vector, every element is labeled by \code{as.character()} of +the value itself and any names the vector carries are ignored. In a list, a +named element is labeled by its name, an unnamed single value by +\code{as.character()} of the value itself, and an unnamed function is an error, +since it has nothing else to be called by. Two elements of a list cannot +share a label.} \item{variant}{One of \code{"data"} (report a single \code{edp} per point, the default) or \code{"estimator"} (report \code{edp_outcome}, \code{edp_treatment}, and @@ -108,6 +118,19 @@ otherwise. Numeric covariates and continuous exposures are continuous dimensions; factor and character covariates, together with binary and categorical exposures, are categorical dimensions. +An intervention need not be a single value every observation is moved to. A +function supplied through \code{values} is a rule \eqn{d} applied to the observed +exposure, so observation \eqn{i} is moved to \eqn{d(a_i)} and the +intervened-on point is \eqn{o^*_i = (l_i, d(a_i))}. Lowering every observed +dose by five units is \code{function(.x) .x - 5}. EDP is then read per natural +observation under the shift rather than at one common target; the algorithm +is unchanged, and only the intervened-on points come from the rule. Because +those points differ across observations, the exposure factor of the kernel +compares every observed exposure against a target that moves with the +observation rather than against a single \eqn{a^*}. A run therefore holds one +\eqn{n \times n} matrix of exposure-kernel values per function intervention, +materialized before any summing and held for the duration of the run. + The half-distances follow the paper's rule of thumb: \code{1 * sd} for each numeric covariate and \code{0.5 * sd} for a continuous exposure. These are marginal scales, so a mid-support gap in a strongly multimodal exposure can be masked @@ -144,6 +167,17 @@ df <- data.frame(dose = dose, x1 = x1) result <- check_edp(df, dose, x1, values = c(0, 1), exposure_type = "continuous") result +# An intervention can also be a function of the observed exposure, read as a +# shift away from the dose each observation received +shifted <- check_edp( + df, + dose, + x1, + values = list(0, 1, lower = function(.x) .x - 1), + exposure_type = "continuous" +) +tidy(shifted) + } \references{ Ring C, Schomaker M (2026). A Diagnostic to Find and Help Combat Stochastic diff --git a/man/pull_plot_data.Rd b/man/pull_plot_data.Rd index 2d557fb..ee08a34 100644 --- a/man/pull_plot_data.Rd +++ b/man/pull_plot_data.Rd @@ -7,7 +7,7 @@ pull_plot_data(x, ...) } \arguments{ -\item{x}{A \link{positivity_diagnostic}.} +\item{x}{A \link{positivity_diagnostic} or a \link{positivity_check}.} \item{...}{Passed to methods. A class drawing more than one view takes \code{type}, the name of the view whose data to return, matched against that @@ -41,9 +41,9 @@ intervention, or one row per measure and intervention for the estimator variant. \code{ymin} and \code{ymax} are the fifth and the ninety-fifth percentile, \code{lower}, \code{middle}, and \code{upper} are the quartiles, \code{n} is the number of observations summarized, and \code{outliers} is a list-column of the values -lying outside the whiskers, which the view draws as points. The estimator variant covers -\code{edp_outcome} and \code{edp_treatment}; \code{ideal_weight} shares no scale with -them and is left to the scatter view. +lying outside the whiskers, which the view draws as points. The estimator +variant covers \code{edp_outcome} and \code{edp_treatment}; \code{ideal_weight} shares no +scale with them and is left to the scatter view. \item \code{"histogram"}, \code{"ecdf"}, and \code{"density"} return the results at their own grain, one row per observation and intervention, with \code{intervention} as a factor whose levels follow the order the interventions were given in. The @@ -70,8 +70,9 @@ For a \code{\link[=check_eta_bias]{check_eta_bias()}} result: \itemize{ \item \code{"bootstrap"} returns one row per bootstrap draw, keyed on the estimand \code{term}, on the truncation \code{level}, and on the \code{label} its facet strip -reads. Each draw carries the \code{truth} its term was aimed at alongside the -\code{estimate}, which the view draws as one line per term. +reads. The view bins the \code{estimate}. Each draw also carries the \code{truth} its +term was aimed at, which the view draws as a reference line, one per term +rather than one per draw, since a term's truth is the same in every draw. \item \code{"sweep"} returns one row per term and truncation level, holding the \code{truncation} the sweep is drawn against, the \code{bias}, and the \code{lower} and \code{upper} ends of the band two Monte Carlo standard errors either side of diff --git a/vignettes/continuous-exposures.Rmd b/vignettes/continuous-exposures.Rmd index 63731ac..49e31fc 100644 --- a/vignettes/continuous-exposures.Rmd +++ b/vignettes/continuous-exposures.Rmd @@ -112,6 +112,47 @@ Effective data points detect the `a2` gap directly, because the kernel measures local support in the joint exposure-covariate space without assuming a form for the conditional density. +### Shifting each observation away from the dose it received + +A candidate value asks what the support looks like when every observation is +moved to the same dose. An intervention can instead be a rule applied to the +dose each observation already received. Passing `values` a list that mixes +single values with functions of the observed exposure asks both questions in one +run. Each function needs a name, which labels its intervention; a single value +labels itself. + +```{r edp-shift} +shift <- check_edp( + pos_violations_long, + a2, + c(l0, l1), + values = list(1, 5, shifted = function(.x) .x - 1) +) + +tidy(shift) |> + group_by(intervention) |> + summarize(edp = round(mean(edp), 1)) +``` + +Under `shifted`, each observation is moved to its own `a2 - 1` rather than to a +common target, so effective data points are read per natural observation: a row +reports the support surrounding the dose that observation would have received +under the shift. The `value` column holds that intervened dose, which now varies +from row to row, and `intervention` is what tells the three apart. The mean under +the shift falls well below either single value because subtracting one carries +observations into two thinly supported places: the 126 drawn between four and +five land in the `(2, 4)` gap, and the 121 drawn below one fall beneath the +observed range, where support is thinner still. + +The default `autoplot()` view boxes effective data points by intervention. Its +whiskers reach the fifth and the ninety-fifth percentile, following Ring and +Schomaker (2026) rather than Tukey's fences, and the values outside them are +drawn as points. + +```{r edp-shift-plot, fig.alt = "Boxplots of effective data points for the interventions one, five, and shifted. The boxes for one and five sit near sixty-five effective data points, while the shifted box is centered near thirty-six and reaches lower, with points beyond the whiskers at both ends of every box."} +autoplot(shift) +``` + ## Hat values: leverage-based extrapolation risk `check_hat_values()` computes the hat-value diagnostic of Moodie and Schulz diff --git a/vignettes/positively.Rmd b/vignettes/positively.Rmd index b363666..b244545 100644 --- a/vignettes/positively.Rmd +++ b/vignettes/positively.Rmd @@ -147,7 +147,7 @@ gives it, so a panel reads well for diagnostics whose default view is compact. PoRT reports a couple of hundred subgroups here, and its default chart is worth a figure of its own, so the panel below holds the other two. -```{r container-plot, eval = requireNamespace("patchwork", quietly = TRUE), message = FALSE, fig.width = 9, fig.height = 3.5, fig.alt = "Effective data points faceted by candidate exposure value, each a broad distribution running from zero to the sample maximum, beside the fraction of nearby opposite-group units faceted by exposure group, each rising from near zero to a peak around four tenths."} +```{r container-plot, eval = requireNamespace("patchwork", quietly = TRUE), message = FALSE, fig.width = 9, fig.height = 3.5, fig.alt = "Boxplots of effective data points at the two candidate exposure values, both spanning most of the range from zero to the largest observed value with points beyond their whiskers, beside histograms of the fraction of nearby opposite-group units faceted by exposure group, each rising from near zero to a peak around four tenths."} compact <- check_positivity( pos_violations, exposure, @@ -223,9 +223,12 @@ Effective data points are a relative measure. They are meaningful when compared across units or across candidate values for a fixed covariate set, not against a universal threshold. Units in a thinly supported region carry a low value at the exposure level they rarely receive, because few observed neighbors surround the -intervened-on point. The `autoplot()` method shows the distribution of effective -data points, faceted by candidate exposure value; the left tail near zero is the -set of units with little observed support under that intervention. +intervened-on point. The `autoplot()` method reads effective data points by +candidate value. The default view boxes them, with whiskers at the fifth and the +ninety-fifth percentile following Ring and Schomaker (2026) and the values +outside them drawn as points. `type = "histogram"` opens the same numbers up as +a distribution faceted by candidate value; the left tail near zero is the set of +units with little observed support under that intervention. ```{r edp-plot, fig.alt = "Histograms of effective data points faceted by candidate exposure value, each with a left tail of low-support observations near zero."} autoplot(edp, type = "histogram") From 9b5ebf0a2cae590f887dec7c55e07d6a89ea61a5 Mon Sep 17 00:00:00 2001 From: Malcolm Barrett Date: Tue, 1 Sep 2026 10:57:16 -0400 Subject: [PATCH 7/7] pin the EDP density figure on R 4.4 or later --- tests/testthat/test-check-edp.R | 3 +++ 1 file changed, 3 insertions(+) diff --git a/tests/testthat/test-check-edp.R b/tests/testthat/test-check-edp.R index de2df4e..f21aefb 100644 --- a/tests/testthat/test-check-edp.R +++ b/tests/testthat/test-check-edp.R @@ -1903,6 +1903,9 @@ test_that("the EDP density view renders as expected", { local_quiet() announce_doppelganger("EDP density by intervention") skip_if_not_installed("ggridges") + # R 4.4.0 corrected density()'s coordinate grid, so the ridgeline curve + # differs on earlier R at SVG precision; the figure is pinned on R >= 4.4. + skip_if(getRversion() < "4.4.0") data <- sim_edp_gaussian(150) res <- check_edp( data,