From 5a26f71ebf732e8f195b7feaaee07be785579d79 Mon Sep 17 00:00:00 2001 From: Ante Bing Date: Fri, 7 Aug 2026 15:43:49 +0000 Subject: [PATCH 1/9] Keep a reported confidence bound inside the range it is compared against Two ways a reported bound left the range of the quantity it describes, both reachable through the exported get_confidence_set() and both giving a probability outside 0 to 1. Center weights are supposed to form a weighted mean over centers, so they have to sum to 1. The exported entry point checked they were non-negative and finite but not that they summed to anything in particular, because the check that does that in validate_inputs() also renormalises, and the confidence set has nothing to renormalise into. Weights summing to 12 then produced an interval of 1.406 to 2.147 for a proportion outcome, and summing to 200, 18.259 to 30.398. They do not scale the outcome, which would at least be undoable: the weights multiply the center-dummy block while the intercept, the intervention contribution and the center characteristics are added once, so the result is a mixture of differently weighted pieces. Such a set is now refused rather than renormalised. The exported function cannot know which of the weights the caller got wrong, and renormalising would also perturb runs that are already correct, since the weights derived from center sample sizes do not always divide to exactly 1. Both entry points refuse the same vectors in the same words, from one shared check. The interval for a binary outcome is computed by the delta method on the probability scale, which is not confined to that scale: a small noisy fit returned a lower bound of -0.106 and an upper bound of 1.049 in the same set. Bounds are now reported within 0 and 1. This is a bound on the report and not a different interval: computing the interval on the linear predictor and transforming it back cannot leave the range at all, but it is a different interval and would move every binary interval the package has ever reported, which is not what a fix for an out-of-range report should do. Both bounds of a row have to be brought into range together. Clamping only the upper bound can leave it below the lower one, and findInterval() refuses a reversed interval outright rather than reporting something odd. Membership in the confidence set is unchanged, except for an outcome goal of exactly 1: findInterval() treats the interval as half-open, so a bound of exactly 1 was already excluded before this change while 1.049 was included. That is pre-existing and not reachable through lago_optimization(), where such a goal is unachievable and the confidence set is suppressed. The stored presentation snapshot contained two of these out-of-range values, a confidence interval ending at 1.005 and a set row with an upper bound of 1.132. It was pinning the defect, and is regenerated. --- R/get_confidence_set.R | 73 +++- R/validate_inputs.R | 88 ++++- man/get_confidence_set.Rd | 13 +- tests/testthat/_snaps/presentation.md | 8 +- tests/testthat/test-confidence-set-contract.R | 317 ++++++++++++++++++ 5 files changed, 485 insertions(+), 14 deletions(-) diff --git a/R/get_confidence_set.R b/R/get_confidence_set.R index 06c8d94..0b72a38 100644 --- a/R/get_confidence_set.R +++ b/R/get_confidence_set.R @@ -12,7 +12,12 @@ #' satisfy the outcome goal for an (weighted) average center. #' The weights need to sum up to 1, and must all be non-negative and finite. #' A weight of 0 is allowed and excludes that center from the average. Only -#' used, and only checked, when include_center_effects is TRUE. +#' used, and only checked, when include_center_effects is TRUE. A vector whose +#' sum is not 1 is refused here rather than renormalised: the interval is +#' computed AT the weights the optimization ran with, so rescaling them would +#' report an interval for a different weighting than the point estimate beside +#' it. lago_optimization() refuses the same vectors, and normalises within the +#' tolerance the ones it accepts, so weights arriving from it always sum to 1. #' @param include_time_effects A boolean. Specifies whether the fixed time #' effects should be included in the outcome model. #' @param time_effect_optimization_value The period the confidence set is @@ -71,7 +76,11 @@ #' rec_int_ci = , +#' For a binary outcome both bounds are confined to [0, 1], since the +#' outcome is a probability; the interval is still the delta-method one and +#' a bound at exactly 0 or 1 is one that has been truncated to the range. +#' No such confinement is applied for a continuous outcome, whose range is +#' not knowable here. NULL when that interval is not computable>, #' cs = = 0.001) { - stop(paste( - "values in center_weights_for_outcome_goal must", - "sum up to 1." - )) - } # The tolerance above says the input was MEANT to be a set of weights. It # does not make it one: the weights multiply the per-center outcomes and are # then summed, so a set summing to 0.999 scales every reported outcome by @@ -1224,3 +1222,81 @@ refuse_invalid_center_weights <- function(center_weights_for_outcome_goal) { } invisible(NULL) } + + +#' refuse_non_unit_weight_sum +#' +#' @description Internal guard for the center weights as a SET: refuses a +#' vector whose sum is not 1, within the tolerance validate_inputs() has always +#' used. +#' +#' @details This is the other way to the wrong number +#' refuse_invalid_center_weights() refuses. That guard rules out a vector no +#' single element of which is a weight; this one rules out a vector every +#' element of which could be, but which is not a set of weights taken together. +#' The reported outcome is sum(weight_i * outcome_i) over the center-level +#' effects, so it is a weighted mean only when the weights sum to 1: a set +#' summing to 12 scales the intervention and center-characteristic contribution +#' of every reported outcome by 12, and 16 weights of 12/16 gave a +#' confidence interval of 1.406 to 2.147 for an outcome that is a proportion. +#' Every weight there is non-negative and finite, so neither of the checks in +#' refuse_invalid_center_weights() applies, and the symptom is the same one it +#' exists to prevent. +#' +#' Shared with the exported get_confidence_set() rather than left at +#' validate_inputs(), which is the same reason the checks in +#' refuse_invalid_center_weights() are shared: get_confidence_set() does not go +#' through validate_inputs(), so a caller reaching it directly had no sum check +#' at all, and the @param text saying the weights must sum to 1 was true of one +#' entry point and not the other. +#' +#' REFUSED rather than renormalised at get_confidence_set(), while +#' validate_inputs() both refuses this and renormalises what it accepts. The two +#' are not inconsistent: both entry points refuse exactly the same vectors, in +#' the same words. What differs is what happens to a vector inside the +#' tolerance, and there validate_inputs() can do more because it OWNS the +#' weights -- it is where they are derived, and every consumer is handed the +#' value it returns. get_confidence_set() is handed weights an optimization has +#' already run with, and is documented as computing the interval at those +#' weights, so renormalising them there would report an interval for a +#' different weighting than the point estimate it is printed beside. That is the +#' same class of defect as reporting the interval of one period beside the +#' estimate of another. +#' +#' Renormalising there would also not be a no-op on correct input, which is what +#' settles it. validate_inputs() stores center_sizes / total_sample_size, and +#' about one such vector in six hundred does not sum to exactly 1 in floating +#' point; dividing by that sum again perturbs every weight by up to one unit in +#' the last place and moves the interval. So renormalising in the name of +#' matching the other entry point would change the numbers of runs that are +#' already correct, and refusing changes none of them: every internal path +#' arrives with a sum of exactly 1, because validate_inputs() has already +#' normalised it. +#' +#' The tolerance is validate_inputs()' own 0.001 rather than a tighter one. +#' It is documented, callers rely on it, and it is what says the input was MEANT +#' to be a set of weights; tightening it here would turn input the package has +#' always accepted into a hard error, and the residual scaling it admits is at +#' most 0.1%, which is the tolerance's own documented meaning. What this closes +#' is the unbounded case. +#' +#' A vector summing to 0, including one that is all zeros, is refused here: it +#' is 1 away from 1. That is also what keeps validate_inputs()' +#' renormalisation from dividing by zero. +#' +#' @param center_weights_for_outcome_goal A numeric vector of center weights, +#' already known to be finite and non-negative. +#' +#' @return Invisibly NULL when the weights sum to 1 within the tolerance. +#' Raises otherwise. +#' +#' @noRd +refuse_non_unit_weight_sum <- function(center_weights_for_outcome_goal) { + if (abs(sum(center_weights_for_outcome_goal) - 1) >= 0.001) { + stop(paste( + "values in center_weights_for_outcome_goal must", + "sum up to 1." + )) + } + invisible(NULL) +} diff --git a/man/get_confidence_set.Rd b/man/get_confidence_set.Rd index 0ee4288..115d081 100644 --- a/man/get_confidence_set.Rd +++ b/man/get_confidence_set.Rd @@ -42,7 +42,12 @@ weights that will be used for calculating recommended interventions that satisfy the outcome goal for an (weighted) average center. The weights need to sum up to 1, and must all be non-negative and finite. A weight of 0 is allowed and excludes that center from the average. Only -used, and only checked, when include_center_effects is TRUE.} +used, and only checked, when include_center_effects is TRUE. A vector whose +sum is not 1 is refused here rather than renormalised: the interval is +computed AT the weights the optimization ran with, so rescaling them would +report an interval for a different weighting than the point estimate beside +it. lago_optimization() refuses the same vectors, and normalises within the +tolerance the ones it accepts, so weights arriving from it always sum to 1.} \item{include_time_effects}{A boolean. Specifies whether the fixed time effects should be included in the outcome model.} @@ -122,7 +127,11 @@ List( rec_int_ci = , + For a binary outcome both bounds are confined to [0, 1], since the + outcome is a probability; the interval is still the delta-method one and + a bound at exactly 0 or 1 is one that has been truncated to the range. + No such confinement is applied for a continuous outcome, whose range is + not knowable here. NULL when that interval is not computable>, cs = = 0 & is.finite(over))) + expect_equal(sum(over), 12) + + # a sum far from 1 is refused, in the words validate_inputs() uses + expect_error(call_cs(over), "sum up to 1") + expect_error(call_cs(rep(200 / n_centers, n_centers)), "sum up to 1") + expect_error(call_cs(rep(0.5 / n_centers, n_centers)), "sum up to 1") + # all-zero weights are a sum of 0, so they are refused here too + expect_error(call_cs(rep(0, n_centers)), "sum up to 1") + + # compliant weights are UNCHANGED, so the guard only removed the refused + # cases. This is the number the scaled vectors were a multiple of. + unit <- call_cs(rep(1 / n_centers, n_centers)) + expect_false(is.null(unit$rec_int_ci)) + expect_identical( + unname(unit$rec_int_ci), + c(0.557, 0.628) + ) + + # the tolerance is validate_inputs()' own 0.001 and is not narrowed: a + # residual sum a hair off 1, which is what renormalised weights can be, is + # still accepted, and so is rounded input the documentation invites + expect_error(call_cs(rep(1 / n_centers, n_centers)), NA) + nudged <- rep(1 / n_centers, n_centers) + nudged[1] <- nudged[1] + 0.0005 + expect_lt(abs(sum(nudged) - 1), 0.001) + expect_error(call_cs(nudged), NA) + + # and the two entry points refuse the same vector in the SAME words, which is + # the property that makes refusing here consistent rather than divergent + primary <- function() { + suppressWarnings(suppressMessages(lago_optimization( + data = pulesa, + outcome_name = "Proportions", + outcome_type = "continuous", + intervention_components = components, + intervention_lower_bounds = c(1, 0.5), + intervention_upper_bounds = c(5, 1), + cost_list_of_vectors = list(c(0, 1), c(0, 1)), + outcome_goal = 0.5, + include_center_effects = TRUE, + center_weights_for_outcome_goal = over, + include_confidence_set = FALSE, + quiet = TRUE + ))) + } + expect_error(primary(), "sum up to 1") + expect_identical( + tryCatch(primary(), error = conditionMessage), + tryCatch(call_cs(over), error = conditionMessage) + ) + + # the weights are only USED when the fixed center effects are included, so a + # caller who is not asking for them passes the scalar default 1 and must not + # be refused for a vector nobody reads. That default sums to 1 anyway, so the + # gate is asserted with a vector that does NOT. + expect_error( + suppressWarnings(get_confidence_set( + predictors_data = pulesa[, components, drop = FALSE], + include_center_effects = FALSE, + center_weights_for_outcome_goal = c(6, 6), + intervention_components = components, + outcome_data = pulesa$Proportions, + fitted_model = glm( + Proportions ~ AccessMedicines + AccessBPMachines, + data = pulesa, family = gaussian(link = "identity") + ), + link = "identity", + outcome_goal = 0.5, + outcome_type = "continuous", + intervention_lower_bounds = c(1, 0.5), + intervention_upper_bounds = c(5, 1), + confidence_set_grid_step_size = c(1, 0.25), + cost_list_of_vectors = list(c(0, 1), c(0, 1)), + rec_int = c(3, 0.75) + )), + NA + ) +}) + + +test_that("a binary outcome's reported interval is confined to [0, 1]", { + # A binary outcome is a probability, so a reported bound of -0.106 or 1.049 is + # not one. The binary branch builds pred +- z*se on the probability scale, + # which is symmetric there and so free to leave [0, 1] on either side, and it + # did on BOTH links: this fixture reported CI_lower_bound -0.106 and + # CI_upper_bound 1.049 in the same returned set. + # + # The fix is a CLAMP, so the assertions below are written against the clamped + # DELTA-METHOD interval and not against a logit-scale one: the interval is + # still the delta-method interval this package has always reported, and + # clamping only bounds what is reported of it. + set.seed(11) + d <- data.frame(a = rep(1:3, each = 6), b = rep_len(c(1, 2), 18)) + d$y <- rbinom(18, 1, 0.5) + components <- c("a", "b") + + call_cs <- function(link, goal = 0.5) { + model <- glm( + y ~ a + b, data = d, family = binomial(link = link) + ) + list(model = model, res = suppressWarnings(get_confidence_set( + predictors_data = d[, components, drop = FALSE], + intervention_components = components, + outcome_data = d$y, + fitted_model = model, + link = link, + outcome_goal = goal, + outcome_type = "binary", + intervention_lower_bounds = c(1, 1), + intervention_upper_bounds = c(3, 2), + confidence_set_grid_step_size = c(1, 1), + cost_list_of_vectors = list(c(0, 1), c(0, 1)), + rec_int = c(2, 1.5) + ))) + } + + critical_value <- qnorm(0.975) + # the hand oracle: the delta-method interval from glm() directly, then + # clamped. It does not go through the package. + hand <- function(model, x, link) { + row <- c(1, x) + eta <- as.numeric(row %*% coef(model)) + se_eta <- sqrt(as.numeric(t(row) %*% vcov(model) %*% row)) + if (link == "identity") { + point <- eta + se_point <- se_eta + } else { + point <- rje::expit(eta) + se_point <- se_eta * point * (1 - point) + } + bounds <- c(point - critical_value * se_point, + point + critical_value * se_point) + round(pmin(pmax(bounds, 0), 1), 3) + } + + for (link in c("logit", "identity")) { + fit <- call_cs(link) + model <- fit$model + res <- fit$res + + # THE assertion: no reported bound is outside [0, 1], on either link + expect_gt(nrow(res$cs), 0) + expect_true(all(res$cs$CI_lower_bound >= 0)) + expect_true(all(res$cs$CI_upper_bound <= 1)) + expect_true(all(res$rec_int_ci >= 0 & res$rec_int_ci <= 1)) + + # the fixture's own precondition: the UNCLAMPED delta-method interval really + # does leave [0, 1] here, on both sides, so the assertion above has + # something to bite on rather than passing vacuously + unclamped <- function(x) { + row <- c(1, x) + eta <- as.numeric(row %*% coef(model)) + se_eta <- sqrt(as.numeric(t(row) %*% vcov(model) %*% row)) + if (link == "identity") { + point <- eta + se_point <- se_eta + } else { + point <- rje::expit(eta) + se_point <- se_eta * point * (1 - point) + } + round(c(point - critical_value * se_point, + point + critical_value * se_point), 3) + } + raw <- vapply( + seq_len(nrow(res$cs)), + function(i) unclamped(c(res$cs$a[i], res$cs$b[i])), + numeric(2) + ) + expect_true(any(raw[1, ] < 0)) + expect_true(any(raw[2, ] > 1)) + + # every row agrees with the hand oracle, at its OWN coordinates, so the + # bounds are the clamped delta-method interval and not some other interval + for (i in seq_len(nrow(res$cs))) { + expect_identical( + unname(c(res$cs$CI_lower_bound[i], res$cs$CI_upper_bound[i])), + hand(model, c(res$cs$a[i], res$cs$b[i]), link) + ) + } + expect_identical( + unname(res$rec_int_ci), hand(model, c(2, 1.5), link) + ) + + # the interval is NOT the logit-scale one transformed back, which is the + # other way to bound it and would have moved every binary number. Asserted + # so the test says which interval is reported, not merely that it is in + # range. + logit_scale <- function(x) { + row <- c(1, x) + eta <- as.numeric(row %*% coef(model)) + se_eta <- sqrt(as.numeric(t(row) %*% vcov(model) %*% row)) + round(rje::expit(c(eta - critical_value * se_eta, + eta + critical_value * se_eta)), 3) + } + expect_false(identical( + unname(res$rec_int_ci), logit_scale(c(2, 1.5)) + )) + + # a clamped bound sits at exactly the boundary, which is what distinguishes + # truncation from an interval that happens to end there + expect_true(any( + res$cs$CI_lower_bound == 0 | res$cs$CI_upper_bound == 1 + )) + + # MEMBERSHIP is unchanged. The goal is a probability, so it lies in the + # range the bounds are clamped to, and clamping a bound to the range it is + # compared against cannot move it across the goal. Asserted directly: + # every returned row still brackets the goal, and every row NOT returned + # still fails to, judged on the UNCLAMPED bounds. + expect_true(all( + res$cs$CI_lower_bound <= 0.5 & res$cs$CI_upper_bound >= 0.5 + )) + grid <- expand.grid(a = 1:3, b = 1:2) + covers_unclamped <- vapply(seq_len(nrow(grid)), function(i) { + bounds <- unclamped(c(grid$a[i], grid$b[i])) + bounds[1] <= 0.5 && 0.5 <= bounds[2] + }, logical(1)) + returned <- paste(res$cs$a, res$cs$b) + expect_setequal( + returned, + paste(grid$a, grid$b)[covers_unclamped] + ) + } +}) + + +test_that("a continuous outcome's interval is NOT clamped to [0, 1]", { + # The counterpart to the test above, pinning the deliberate asymmetry. The + # binary branch can confine its bounds because "binary" fixes the range at + # [0, 1]; a continuous outcome is a mean on whatever scale the caller's data + # is on, and nothing here knows what that is. "continuous" with a proportion + # is supported, and so is a count or a duration, and the outcome goal is only + # required to be numeric -- so an out-of-[0, 1] bound is not knowably wrong + # for a continuous outcome and clamping it would corrupt every outcome that is + # not a proportion. + # + # This asserts the range is left alone on a fixture whose outcome is far + # outside [0, 1], which a clamp would silently destroy. + set.seed(5) + d <- data.frame(x1 = rep(1:6, each = 4), x2 = rep_len(c(1, 2), 24)) + d$y <- 40 + 3 * d$x1 + 2 * d$x2 + rnorm(24, 0, 1.5) + expect_true(all(d$y > 1)) + + res <- suppressWarnings(suppressMessages(lago_optimization( + data = d, + outcome_name = "y", + outcome_type = "continuous", + intervention_components = c("x1", "x2"), + intervention_lower_bounds = c(1, 1), + intervention_upper_bounds = c(6, 2), + cost_list_of_vectors = list(c(0, 1), c(0, 1)), + outcome_goal = 55, + confidence_set_grid_step_size = c(1, 1), + quiet = TRUE + ))) + + # the bounds are on the outcome's own scale, well outside [0, 1], and are + # reported as they are + expect_gt(nrow(res$cs), 0) + expect_true(all(res$cs$CI_lower_bound > 1)) + expect_true(all(res$cs$CI_upper_bound > 1)) + expect_true(all(res$est_outcome_ci > 1)) + # a clamp would have collapsed every one of them to exactly 1 + expect_false(any(res$cs$CI_upper_bound == 1)) +}) From cb884d3faf6e386d408a0fb226b947b2bbed3d35 Mon Sep 17 00:00:00 2001 From: Ante Bing Date: Fri, 7 Aug 2026 16:16:52 +0000 Subject: [PATCH 2/9] Clamp what is reported, not what membership is decided on The clamp was applied to the matrix the confidence set is then computed from, so confining a bound changed which interventions qualify. findInterval() treats an interval as right-open, so for an outcome goal of exactly 1 a bound brought down from 1.049 to 1 puts the goal at the closed end and the row falls out: an intervention whose delta-method interval covers the goal stopped qualifying because of how its bound was reported. The claim in the previous commit that membership was unchanged held only for goals inside the range. Clamping is a statement about the report, so it is now applied to the report and nothing else. The bounds carried on the returned set and in the interval at the recommended intervention come from the confined matrix, while membership and the missing-bound filter read the interval as computed. Verified against the unclamped tree at goals of 0.5 and 1: the same rows qualify in both, while the reported bounds go from -0.106 and 1.049 to 0 and 1. A comment claimed findInterval() rejects an inverted interval outright. It does not, it reports against a vector it was told is sorted. Clamping both bounds together is still right, for that reason rather than the stated one. --- R/get_confidence_set.R | 34 ++++++++++++++++++++++++---------- 1 file changed, 24 insertions(+), 10 deletions(-) diff --git a/R/get_confidence_set.R b/R/get_confidence_set.R index 0b72a38..069ddea 100644 --- a/R/get_confidence_set.R +++ b/R/get_confidence_set.R @@ -610,11 +610,21 @@ get_confidence_set <- function( # the smaller change fixes exactly that and nothing else. # # Both bounds are clamped in ONE step so the pair cannot cross. Clamping - # only the upper bound would turn [1.02, 1.30] into [1.02, 1] -- an inverted - # interval, which findInterval() below rejects outright as unsorted rather - # than merely reporting oddly. pmin/pmax preserve the matrix's shape and - # leave an NA bound NA, so the complete.cases() filter below is unaffected. - ci_prob_all <- pmin(pmax(ci_prob_all, 0), 1) + # only the upper bound would turn [1.02, 1.30] into [1.02, 1], an inverted + # interval, and findInterval() does not reject that, it silently reports + # against a vector it was told is sorted. pmin/pmax preserve the matrix's + # shape and leave an NA bound NA, so the complete.cases() filter below is + # unaffected. + # + # Clamped SEPARATELY from the interval membership is decided on, which is + # the unclamped one. findInterval() treats the interval as right-open, so + # for an outcome goal of exactly 1 a bound clamped from 1.049 down to 1 puts + # the goal at the closed end and the row falls out of the set: an + # intervention whose delta-method interval covers the goal would stop + # qualifying because of how its bound is REPORTED. Clamping is a statement + # about the report, so it is applied to what is reported and to nothing + # else. + reported_ci_prob_all <- pmin(pmax(ci_prob_all, 0), 1) } else if (outcome_type == "continuous") { # link is either "logit" or "identity", the only links the outcome # machinery implements, see supported_outcome_links(). @@ -961,12 +971,16 @@ get_confidence_set <- function( # property the note above describes, arrived at by construction. ci_prob_all <- expit(cbind(lb_prob_all, ub_prob_all)) } + # a continuous outcome has no range to confine a bound to, so what is + # reported is what was computed. See the note above. + reported_ci_prob_all <- ci_prob_all } else { # ci_prob_all is assigned in the two branches above and read by the shared # code below, so an unrecognised outcome type would reach that code with it # undefined and fail on "object 'ci_prob_all' not found" instead of saying # what is wrong. lago_optimization() validates outcome_type, so this is - # reachable only through a direct call. + # reachable only through a direct call. reported_ci_prob_all is defaulted + # below for the same reason: only the binary branch narrows what it reports. stop(paste0( "'outcome_type' must be either \"binary\" or \"continuous\", not \"", outcome_type, "\"." @@ -995,8 +1009,8 @@ get_confidence_set <- function( # rec_int inside the confidence set. NULL when it could not be computed. rec_int_ci <- if (valid_rows[1]) { c( - lower = round(ci_prob_all[1, 1], 3), - upper = round(ci_prob_all[1, 2], 3) + lower = round(reported_ci_prob_all[1, 1], 3), + upper = round(reported_ci_prob_all[1, 2], 3) ) } else { NULL @@ -1068,8 +1082,8 @@ get_confidence_set <- function( # cs_row_indices are grid_x row numbers, and ci_prob_all is indexed by the # same row numbers, so the bounds line up with the cs rows they belong to. cs_output_names <- names(cs) - ci_lower_bound <- round(ci_prob_all[cs_row_indices, 1], 3) - ci_upper_bound <- round(ci_prob_all[cs_row_indices, 2], 3) + ci_lower_bound <- round(reported_ci_prob_all[cs_row_indices, 1], 3) + ci_upper_bound <- round(reported_ci_prob_all[cs_row_indices, 2], 3) cs <- cbind( cs, From 4c17bdf309c10bb8dfd6e4f3f2458dc174a3f33e Mon Sep 17 00:00:00 2001 From: Ante Bing Date: Fri, 7 Aug 2026 16:20:55 +0000 Subject: [PATCH 3/9] Report no interval where the interval lies wholly outside the range Confining a bound to [0, 1] is only meaningful for an interval that has a part inside it. For one that lies entirely above 1 or entirely below 0 the intersection is empty, and clamping reported that emptiness as a non-empty interval: [1.02, 1.30] became [1, 1], a 95% interval of zero width, reported beside a point estimate it excludes. That is the invariant the previous release was fixed to restore, so clamping reintroduced it. It is reachable on an identity link, where a binomial model is a linear probability model and its prediction is the linear predictor, which is not confined to [0, 1] at all. Two reviewers found it independently, one by measuring 78 of 84 runs reporting a point estimate outside its own interval where the unclamped tree reports none, and the existing assertion that a point estimate lies inside its interval did not catch it because its fixture stays well inside the range and never reaches the clamp. Such rows now report no interval. NA is already what a bound that could not be computed carries here, and the consumers already handle it: the row is excluded from the confidence set by the same filter that drops an uncomputable bound, and the estimated outcome is reported with no interval rather than with an impossible one. Rows that straddle a boundary are still clamped, and rows inside the range are untouched. --- R/get_confidence_set.R | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/R/get_confidence_set.R b/R/get_confidence_set.R index 069ddea..09e4fe9 100644 --- a/R/get_confidence_set.R +++ b/R/get_confidence_set.R @@ -624,7 +624,20 @@ get_confidence_set <- function( # qualifying because of how its bound is REPORTED. Clamping is a statement # about the report, so it is applied to what is reported and to nothing # else. + # + # An interval lying WHOLLY outside [0, 1] has no part inside it to report, + # and clamping it would report that empty intersection as a non-empty + # interval: [1.02, 1.30] becomes [1, 1], a 95% interval of zero width which + # excludes the point estimate reported beside it. That is reachable on an + # identity link, where the model is a linear probability model and its + # prediction is the linear predictor, unbounded. Those rows report no + # interval instead. NA is already what a bound that could not be computed + # carries here, and every consumer below already handles it: the row is + # excluded from the set by the same filter, and lago_optimization() reports + # no interval rather than an impossible one. + outside_range <- ci_prob_all[, 1] > 1 | ci_prob_all[, 2] < 0 reported_ci_prob_all <- pmin(pmax(ci_prob_all, 0), 1) + reported_ci_prob_all[which(outside_range), ] <- NA_real_ } else if (outcome_type == "continuous") { # link is either "logit" or "identity", the only links the outcome # machinery implements, see supported_outcome_links(). From 39793e7ca2ed3f82bc360d13a093c792acca1713 Mon Sep 17 00:00:00 2001 From: Ante Bing Date: Fri, 7 Aug 2026 16:36:31 +0000 Subject: [PATCH 4/9] Confine the interval only where the estimate it belongs to is confined Confining the reported bounds on both links was wrong on one of them. A binomial model fitted with an identity link is a linear probability model, and its point estimate is the linear predictor: extrapolate far enough along a positive slope and the estimate itself exceeds 1, which the package already reports as it stands. Confining only the interval around such an estimate reports an interval that excludes it, and where the whole interval is above 1 a zero-width one. On the logit link the estimate is expit of the linear predictor and is inside the range by construction, so confining the interval brings it into agreement with the estimate rather than out of it. Same operation, opposite effect, so it is applied only where it is right. Gating on the link also removes the case an earlier commit here handled by reporting no interval at all: an interval centred on an estimate inside the range always has a part inside it, so at most one of its bounds is ever confined and none is ever confined to nothing. The estimate leaving the range on an identity link is a defect of its own, and a larger one, since nothing constrains or flags it. It is not this one and is left for its own change. The test's oracle confined its bounds on both links, so its identity arm agreed with any implementation that did the same and could not observe the choice. The oracle now follows the link, the identity arm asserts that every reported interval contains the estimate it belongs to, which is the property confining them would break, and the logit arm keeps asserting the bounds are in range. The interval reported for the recommended intervention is asserted at a recommendation whose own interval leaves the range, since the previous fixture's was inside it and only a stored snapshot was holding that half of the fix. --- R/get_confidence_set.R | 46 ++++---- tests/testthat/test-confidence-set-contract.R | 105 ++++++++++++++++-- 2 files changed, 123 insertions(+), 28 deletions(-) diff --git a/R/get_confidence_set.R b/R/get_confidence_set.R index 09e4fe9..25533b0 100644 --- a/R/get_confidence_set.R +++ b/R/get_confidence_set.R @@ -585,12 +585,23 @@ get_confidence_set <- function( # the shared code below turns these bounds into the confidence set and # its size, for both outcome types ci_prob_all <- cbind(lb_prob_all, ub_prob_all) - # a binary outcome is a probability, so a REPORTED bound of -0.106 or 1.049 - # is not a probability and cannot be one. Both bounds are confined to - # [0, 1] here, on both links: pred +- z*se is symmetric on the outcome scale - # and so is free to leave the range on either side, and it did on both links - # -- a small noisy logit fit reported CI_lower_bound -0.106 and - # CI_upper_bound 1.049 in the same returned set. + # a binary outcome fitted on the logit link is reported as a probability, so + # a REPORTED bound of -0.106 or 1.049 is not one and cannot be: the point + # estimate is expit(eta) and is inside [0, 1] by construction, while + # pred +- z*se is symmetric on the outcome scale and so is free to leave it. + # A small noisy logit fit reported CI_lower_bound -0.106 and CI_upper_bound + # 1.049 in the same returned set. Confining the bounds brings the interval + # into agreement with the estimate it belongs to. + # + # NOT done on the identity link, where the same operation would do the + # opposite. There the model is a linear probability model and the point + # estimate is the linear predictor itself, which extrapolates outside [0, 1] + # and is already reported that way: an estimated outcome of 1.0588 is + # printed as it stands. Confining only the interval around it would report + # an interval that excludes its own estimate, and where the whole interval + # is above 1 a zero-width one. The estimate leaving the range is a defect in + # its own right and is not this one, so this branch reports what it computed + # and leaves that alone rather than papering over half of it. # # This is a CLAMP: the interval is still the delta-method interval, and # clamping only bounds what is reported of it. It is not a different @@ -625,19 +636,16 @@ get_confidence_set <- function( # about the report, so it is applied to what is reported and to nothing # else. # - # An interval lying WHOLLY outside [0, 1] has no part inside it to report, - # and clamping it would report that empty intersection as a non-empty - # interval: [1.02, 1.30] becomes [1, 1], a 95% interval of zero width which - # excludes the point estimate reported beside it. That is reachable on an - # identity link, where the model is a linear probability model and its - # prediction is the linear predictor, unbounded. Those rows report no - # interval instead. NA is already what a bound that could not be computed - # carries here, and every consumer below already handles it: the row is - # excluded from the set by the same filter, and lago_optimization() reports - # no interval rather than an impossible one. - outside_range <- ci_prob_all[, 1] > 1 | ci_prob_all[, 2] < 0 - reported_ci_prob_all <- pmin(pmax(ci_prob_all, 0), 1) - reported_ci_prob_all[which(outside_range), ] <- NA_real_ + # Gating on the link also removes the case of an interval lying WHOLLY + # outside the range, which has no part inside it to report and which + # clamping would report as the non-empty [1, 1]. On the logit link the + # estimate is inside [0, 1], so an interval centred on it always has a part + # inside too, and at most one of its bounds is ever confined. + reported_ci_prob_all <- if (link == "logit") { + pmin(pmax(ci_prob_all, 0), 1) + } else { + ci_prob_all + } } else if (outcome_type == "continuous") { # link is either "logit" or "identity", the only links the outcome # machinery implements, see supported_outcome_links(). diff --git a/tests/testthat/test-confidence-set-contract.R b/tests/testthat/test-confidence-set-contract.R index 85e2e4e..3d2ea52 100644 --- a/tests/testthat/test-confidence-set-contract.R +++ b/tests/testthat/test-confidence-set-contract.R @@ -1183,7 +1183,16 @@ test_that("a binary outcome's reported interval is confined to [0, 1]", { } bounds <- c(point - critical_value * se_point, point + critical_value * se_point) - round(pmin(pmax(bounds, 0), 1), 3) + # confined on the logit link only, matching the code: there the point + # estimate is expit(eta) and inside [0, 1], so confining the interval brings + # it into agreement with the estimate. On the identity link the estimate is + # the linear predictor and is itself unbounded, so confining the interval + # around it would report an interval excluding its own estimate. An oracle + # that clamped on both links would agree with any implementation that did. + if (link == "logit") { + bounds <- pmin(pmax(bounds, 0), 1) + } + round(bounds, 3) } for (link in c("logit", "identity")) { @@ -1191,11 +1200,27 @@ test_that("a binary outcome's reported interval is confined to [0, 1]", { model <- fit$model res <- fit$res - # THE assertion: no reported bound is outside [0, 1], on either link + # THE assertion, on the logit link: no reported bound is outside [0, 1]. + # On the identity link the opposite is required, since the estimate itself + # is not confined there and an interval that excluded its own estimate would + # be worse than one leaving the range. Both are asserted, so neither arm + # passes by agreeing with whatever the code happens to do. expect_gt(nrow(res$cs), 0) - expect_true(all(res$cs$CI_lower_bound >= 0)) - expect_true(all(res$cs$CI_upper_bound <= 1)) - expect_true(all(res$rec_int_ci >= 0 & res$rec_int_ci <= 1)) + if (link == "logit") { + expect_true(all(res$cs$CI_lower_bound >= 0)) + expect_true(all(res$cs$CI_upper_bound <= 1)) + expect_true(all(res$rec_int_ci >= 0 & res$rec_int_ci <= 1)) + } else { + # every reported interval contains the estimate it belongs to, which is + # the property confining them on this link would break + for (i in seq_len(nrow(res$cs))) { + point <- as.numeric( + c(1, res$cs$a[i], res$cs$b[i]) %*% coef(model) + ) + expect_gte(round(point, 3), res$cs$CI_lower_bound[i]) + expect_lte(round(point, 3), res$cs$CI_upper_bound[i]) + } + } # the fixture's own precondition: the UNCLAMPED delta-method interval really # does leave [0, 1] here, on both sides, so the assertion above has @@ -1234,6 +1259,35 @@ test_that("a binary outcome's reported interval is confined to [0, 1]", { unname(res$rec_int_ci), hand(model, c(2, 1.5), link) ) + # rec_int_ci at a recommendation whose UNCLAMPED interval leaves the range. + # The recommendation above is inside it, so that assertion cannot observe + # the clamp on this field at all, and rec_int_ci is row 1 of the same matrix + # and is what the printed interval for the estimated outcome comes from: it + # has to be confined too, and only a stored snapshot was holding that. + out_of_range <- suppressWarnings(get_confidence_set( + predictors_data = d[, c("a", "b"), drop = FALSE], + intervention_components = c("a", "b"), + outcome_data = d$y, fitted_model = model, link = link, + outcome_goal = 0.5, outcome_type = "binary", + intervention_lower_bounds = c(1, 1), + intervention_upper_bounds = c(3, 2), + confidence_set_grid_step_size = c(1, 1), + cost_list_of_vectors = list(c(0, 1), c(0, 1)), + rec_int = c(1, 1) + )) + expect_true(unclamped(c(1, 1))[1] < 0) + expect_identical(unname(out_of_range$rec_int_ci), hand(model, c(1, 1), link)) + if (link == "logit") { + expect_true(all( + out_of_range$rec_int_ci >= 0 & out_of_range$rec_int_ci <= 1 + )) + } else { + # not confined here, so it is the interval as computed + expect_identical( + unname(out_of_range$rec_int_ci), unclamped(c(1, 1)) + ) + } + # the interval is NOT the logit-scale one transformed back, which is the # other way to bound it and would have moved every binary number. Asserted # so the test says which interval is reported, not merely that it is in @@ -1250,10 +1304,43 @@ test_that("a binary outcome's reported interval is confined to [0, 1]", { )) # a clamped bound sits at exactly the boundary, which is what distinguishes - # truncation from an interval that happens to end there - expect_true(any( - res$cs$CI_lower_bound == 0 | res$cs$CI_upper_bound == 1 - )) + # truncation from an interval that happens to end there. Exactly 0 and + # exactly 1, not nearly: findInterval() treats the interval as [lower, + # upper), so a bound confined to 1e-6 rather than 0 would report a + # probability that is not one AND drop every row from a goal of exactly 0. + if (link == "logit") { + for (i in seq_len(nrow(res$cs))) { + bounds <- unclamped(c(res$cs$a[i], res$cs$b[i])) + if (bounds[1] < 0) { + expect_identical(res$cs$CI_lower_bound[i], 0) + } + if (bounds[2] > 1) { + expect_identical(res$cs$CI_upper_bound[i], 1) + } + } + expect_true(any( + res$cs$CI_lower_bound == 0 | res$cs$CI_upper_bound == 1 + )) + } else { + # on the identity link nothing is confined, so a bound that left the range + # is reported as computed. Asserted so this arm cannot pass under a design + # that clamps here. + out <- vapply( + seq_len(nrow(res$cs)), + function(i) { + b <- unclamped(c(res$cs$a[i], res$cs$b[i])) + b[1] < 0 || b[2] > 1 + }, + logical(1) + ) + if (any(out)) { + i <- which(out)[1] + expect_identical( + unname(c(res$cs$CI_lower_bound[i], res$cs$CI_upper_bound[i])), + unclamped(c(res$cs$a[i], res$cs$b[i])) + ) + } + } # MEMBERSHIP is unchanged. The goal is a probability, so it lies in the # range the bounds are clamped to, and clamping a bound to the range it is From a60718f140ba3fb03f1282ff08a50f5595db7b4c Mon Sep 17 00:00:00 2001 From: Ante Bing Date: Fri, 7 Aug 2026 16:41:51 +0000 Subject: [PATCH 5/9] Pin membership at the goal where confining a bound would move it Membership was asserted at an interior goal only, where no bound is confined and the two matrices cannot come apart. The case that distinguishes them is an outcome goal of exactly 1: findInterval() treats an interval as right-open, so a bound brought down to 1 puts the goal at the closed end and the row stops qualifying because of how it is reported. That is the defect fixed two commits ago, and nothing was holding it fixed, which is how it went in unnoticed in the first place. Verified by recombining the two matrices in a throwaway tree, deciding membership from the confined one again: the checks fail there. --- tests/testthat/test-confidence-set-contract.R | 21 +++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/tests/testthat/test-confidence-set-contract.R b/tests/testthat/test-confidence-set-contract.R index 3d2ea52..85e41c6 100644 --- a/tests/testthat/test-confidence-set-contract.R +++ b/tests/testthat/test-confidence-set-contract.R @@ -1360,6 +1360,27 @@ test_that("a binary outcome's reported interval is confined to [0, 1]", { returned, paste(grid$a, grid$b)[covers_unclamped] ) + + # AT A GOAL OF EXACTLY 1, which is where confining a bound and deciding + # membership from it come apart. findInterval() treats the interval as + # [lower, upper), so a bound brought down from 1.049 to 1 puts the goal at + # the closed end and the row stops qualifying: an intervention whose + # computed interval covers the goal would drop out because of how its bound + # is reported. Membership therefore reads the interval as computed, and this + # is what fails if the two are ever recombined. The goal 0.5 case above + # cannot see it, since no bound is confined near an interior goal. + at_one <- call_cs(link, goal = 1)$res + covers_one <- vapply(seq_len(nrow(grid)), function(i) { + bounds <- unclamped(c(grid$a[i], grid$b[i])) + bounds[1] <= 1 && 1 <= bounds[2] + }, logical(1)) + if (any(covers_one)) { + expect_false(is.null(at_one$cs)) + expect_setequal( + paste(at_one$cs$a, at_one$cs$b), + paste(grid$a, grid$b)[covers_one] + ) + } } }) From ee2b10491990e1c1df5f5dabb287fd875909507b Mon Sep 17 00:00:00 2001 From: Ante Bing Date: Fri, 7 Aug 2026 17:10:20 +0000 Subject: [PATCH 6/9] Document the confinement as conditional, since that is what it is The return documentation still promised that both bounds are confined to [0, 1] for a binary outcome, unconditionally. Gating the confinement on the link left that claim behind, so the documented contract said a reported bound is a probability while the function returns -0.169 and 1.011 on an identity-link fit. That is the same shape as the defect this branch is about, a stated guarantee about a reported probability that the value does not honour, and R CMD check cannot see it because it is prose. It now says which link confines and why the other does not. The identity arm's assertions that a reported interval contains its own estimate hold on that fixture whether or not the interval is confined, because no estimate there leaves [0, 1]: the largest is 0.594. What pins the choice is the hand oracle, which follows the link and so disagrees with a confined bound. The assertions are kept, because they state the property the choice exists to protect and would bite on a fixture that extrapolates far enough, and the comment now says which of the two is doing the work. --- R/get_confidence_set.R | 14 +++++++++----- man/get_confidence_set.Rd | 14 +++++++++----- tests/testthat/test-confidence-set-contract.R | 9 ++++++++- 3 files changed, 26 insertions(+), 11 deletions(-) diff --git a/R/get_confidence_set.R b/R/get_confidence_set.R index 25533b0..0bb11df 100644 --- a/R/get_confidence_set.R +++ b/R/get_confidence_set.R @@ -76,11 +76,15 @@ #' rec_int_ci = , +#' For a binary outcome on the logit link both bounds are confined to +#' [0, 1], where the estimate is a probability by construction and the +#' interval around it therefore belongs in that range; the interval is +#' still the delta-method one and a bound at exactly 0 or 1 is one that has +#' been truncated to the range. Not confined on the identity link, where the +#' estimate is the linear predictor and is itself unbounded, so confining +#' the interval would report one that excludes its own estimate. Not +#' confined for a continuous outcome either, whose range is not knowable +#' here. NULL when that interval is not computable>, #' cs = , + For a binary outcome on the logit link both bounds are confined to + [0, 1], where the estimate is a probability by construction and the + interval around it therefore belongs in that range; the interval is + still the delta-method one and a bound at exactly 0 or 1 is one that has + been truncated to the range. Not confined on the identity link, where the + estimate is the linear predictor and is itself unbounded, so confining + the interval would report one that excludes its own estimate. Not + confined for a continuous outcome either, whose range is not knowable + here. NULL when that interval is not computable>, cs = = 0 & res$rec_int_ci <= 1)) } else { # every reported interval contains the estimate it belongs to, which is - # the property confining them on this link would break + # the property confining them on this link would break. + # + # On THIS fixture no estimate leaves [0, 1] -- the largest is 0.594 -- so + # these hold whether or not the interval is confined, and are not what + # pins the choice. The hand oracle below is: it follows the link, so a + # confined bound disagrees with it. These are kept because they state the + # property the choice exists to protect, and they bite on a fixture that + # extrapolates far enough for the estimate itself to leave the range. for (i in seq_len(nrow(res$cs))) { point <- as.numeric( c(1, res$cs$a[i], res$cs$b[i]) %*% coef(model) From bec36439271940bc9d5cdf3d3a2c45b0e12e0c00 Mon Sep 17 00:00:00 2001 From: Ante Bing Date: Fri, 7 Aug 2026 19:20:07 +0000 Subject: [PATCH 7/9] Warn where a reported estimate is an extrapolation, on both routes it can be The confidence-bound work stopped at confining what could be confined and left the estimate itself unaddressed on two paths. Both are now flagged, and neither changes a returned value: a caller who was getting a number still gets the same number, with a warning attached where that number is an extrapolation. A binary outcome on the identity link is a linear probability model, so its estimate is the linear predictor and is not confined to [0, 1]. Extrapolating along a positive slope past the fitted range gives an estimate above 1, which glm() does not object to because only the extrapolated grid predictions leave the range while the fitted values stay inside it. On the bundled shape at an unreachable goal the reported estimated outcome was 1.54, a probability that is not one. lago_optimization() now warns when a binary outcome's reported estimate leaves [0, 1]. The warning lives there and not in get_outcome(), which is the one function that computes the estimate but does not receive the outcome type, so it cannot know [0, 1] is the right range, and which runs on the negated scale under minimize, where every value it returns is negative regardless. Only lago_optimization() has both the outcome type and the estimate on the caller's own scale, and it fires once. The outcome is also predicted with any additional covariates held at 0. When a covariate's observed values never reach 0, that is a prediction for a covariate value that never occurred, and the reported outcome and interval are an extrapolation: a covariate observed in the fives, held at 0, shifted the reported interval. get_confidence_set() now warns, once, naming each numeric covariate whose observed range excludes 0 and that range. A factor or character covariate held at 0 is at its reference level, which is observed, so it does not warn. The value the covariate is held at is unchanged, so no returned number moves; the warning only says the number is an extrapolation. Neither the estimate leaving the range nor the covariate default is altered here. Both are reported as before and now flagged, which was the choice: the first because clamping the estimate would change which intervention the optimizer recommends, since get_outcome() drives the optimization and the goal comparison, and the second because rescaling or moving the held-at value would change every reported number for a run using additional covariates. --- R/get_confidence_set.R | 98 ++++++- R/get_outcome.R | 154 +++++++++++ R/lago_optimization.R | 19 ++ man/get_confidence_set.Rd | 9 +- tests/testthat/test-covariate-support.R | 241 ++++++++++++++++++ .../test-goal-modes-and-diagnostics.R | 214 ++++++++++++++++ 6 files changed, 729 insertions(+), 6 deletions(-) create mode 100644 tests/testthat/test-covariate-support.R diff --git a/R/get_confidence_set.R b/R/get_confidence_set.R index 0bb11df..1f2fad9 100644 --- a/R/get_confidence_set.R +++ b/R/get_confidence_set.R @@ -82,9 +82,12 @@ #' still the delta-method one and a bound at exactly 0 or 1 is one that has #' been truncated to the range. Not confined on the identity link, where the #' estimate is the linear predictor and is itself unbounded, so confining -#' the interval would report one that excludes its own estimate. Not -#' confined for a continuous outcome either, whose range is not knowable -#' here. NULL when that interval is not computable>, +#' the interval would report one that excludes its own estimate. That +#' estimate leaving [0, 1] is a defect in its own right and not this one; +#' lago_optimization() now warns when it does, so it is flagged rather than +#' silent, and the interval is still reported as computed for the same +#' reason as before. Not confined for a continuous outcome either, whose +#' range is not knowable here. NULL when that interval is not computable>, #' cs = 0) { + # the covariates are held at 0 for the prediction (see below). A numeric + # covariate whose observed values never reach 0 is then predicted at a + # value that does not occur in the data, so the reported outcome and + # interval are an extrapolation. Warn about that, once, naming each such + # covariate and its range. Diagnostic only: no assembled value is changed. + warn_covariates_held_off_support(additional_covariates, predictors_data) # one column per coefficient a covariate expands into, not one column per # covariate: a factor or character covariate is one coefficient per # non-reference level, so it needs that many columns. Every column is 0, @@ -607,6 +616,14 @@ get_confidence_set <- function( # its own right and is not this one, so this branch reports what it computed # and leaves that alone rather than papering over half of it. # + # That estimate is no longer silent about it: lago_optimization() warns + # when a binary outcome's estimate leaves [0, 1], naming the extrapolation + # that causes it, and counts these bounds in what it says is affected. The + # two decisions are the same statement from two sides -- neither alters a + # number, and the user is told the report is not a probability -- so this + # branch is unchanged by it. The warning triggers on the ESTIMATE and not on + # a bound, precisely so that it does not re-open the choice made here. + # # This is a CLAMP: the interval is still the delta-method interval, and # clamping only bounds what is reported of it. It is not a different # interval and does not have that interval's coverage; a bound sitting at @@ -1372,6 +1389,81 @@ center_characteristic_coef_names <- function(center_characteristics, unlist(resolved, use.names = FALSE) } +#' warn_covariates_held_off_support +#' +#' @description Internal diagnostic. The confidence set prediction holds every +#' additional covariate at 0. Warns, once, when a NUMERIC additional covariate +#' is held there outside the range it was observed over, since the reported +#' outcome and interval are then an extrapolation to a covariate value that +#' never occurs in the data. +#' +#' @details Only a numeric covariate can be held off its support this way. A +#' factor, character or logical covariate held at 0 sits at its reference level +#' (all dummies 0, or FALSE for a logical), which is an OBSERVED level and so +#' not an extrapolation, so those are left alone. A numeric covariate whose +#' observed range CONTAINS 0 -- lower <= 0 <= upper -- takes the value 0 +#' somewhere in the data and is likewise not extrapolated, so it does not warn +#' either. Only a numeric covariate whose observed range excludes 0 warrants +#' the warning. +#' +#' Fires once per call, listing every offending covariate with its observed +#' range, rather than once per grid row. Changes no value the caller assembles: +#' the covariates are still held at 0. A covariate name is matched to its +#' column in predictors_data with backticks stripped, and a covariate with no +#' such column is skipped, since its support cannot be read. +#' +#' @param additional_covariates A character vector, the names of the additional +#' covariate columns held at 0 for the prediction. +#' @param predictors_data A data.frame, the input data whose columns hold the +#' observed values of those covariates. +#' +#' @return Invisibly NULL. Raises a warning as a side effect when at least one +#' numeric covariate is held at 0 outside its observed range. +#' +#' @noRd +warn_covariates_held_off_support <- function(additional_covariates, + predictors_data) { + offenders <- character(0) + for (cov in additional_covariates) { + column <- gsub("`", "", cov) + if (!column %in% names(predictors_data)) { + next + } + values <- predictors_data[[column]] + # a factor, character or logical covariate held at 0 is at its reference + # level, an observed level, so it is not an extrapolation + if (!is.numeric(values)) { + next + } + observed <- range(values, na.rm = TRUE) + if (!all(is.finite(observed))) { + next + } + # 0 inside the observed range is a value the covariate takes, so holding + # it there is not an extrapolation + if (observed[1] <= 0 && observed[2] >= 0) { + next + } + offenders <- c(offenders, paste0( + column, " (observed range [", signif(observed[1], 6), ", ", + signif(observed[2], 6), "])" + )) + } + if (length(offenders) == 0) { + return(invisible(NULL)) + } + warning(paste0( + "The additional covariate(s) ", paste(offenders, collapse = ", "), + " are held at 0 to compute the confidence set, but 0 is outside their ", + "observed range, so the reported estimated outcome and its interval are ", + "an extrapolation to a covariate value that never occurs in the data. No ", + "reported value has been altered: the covariate is still held at 0. To ", + "report the outcome at an observed covariate value, center the covariate ", + "so that 0 falls within its range, or drop it from the outcome model." + )) + invisible(NULL) +} + # the levels of a column in the order glm() builds its dummies from, i.e. with # the reference level first. glm() uses the factor's levels, dropping any that # no row takes, and orders an unclassed character column's levels by sorting it, diff --git a/R/get_outcome.R b/R/get_outcome.R index aacaaff..db3e71b 100644 --- a/R/get_outcome.R +++ b/R/get_outcome.R @@ -82,6 +82,160 @@ get_outcome <- function( } +#' warn_if_outcome_outside_range +#' +#' @description Non-fatal check that the estimated outcome REPORTED for a +#' binary outcome is a probability, i.e. that it lies in [0, 1]. Warns, never +#' raises, and never alters a value. +#' +#' @details A binary outcome's estimate is a probability, and on the "logit" +#' link it is expit() of the linear predictor and is inside [0, 1] by +#' construction. On the "identity" link the model is a linear probability model +#' and the estimate IS the linear predictor, so it is not confined to anything. +#' A fit whose every fitted value on the DATA is a probability still +#' extrapolates outside [0, 1] at an intervention outside the range its +#' components were fitted over, which is exactly what intervention bounds +#' reaching beyond that range ask for. glm() does not object, because it only +#' ever sees the data. So "Estimated outcome: 1.5351" was reported for a +#' binary outcome with no error and no warning. +#' +#' WARN, rather than refuse or clamp, and the two rejected options are worth +#' recording because neither is harmless. +#' +#' Refusing would reject a fit that is legitimate over its own data range: a +#' linear probability model is a defensible choice, lago_optimization() +#' accepts it deliberately, and the estimate is only out of range for +#' interventions outside the observed support. +#' +#' Clamping the estimate would be worse than either. get_outcome() is what +#' every optimizer and the goal comparison are driven by, so clamping it there +#' would change WHICH intervention is recommended and would flatten the +#' objective above the boundary, turning a reporting problem into an +#' optimization one. Clamping only the reported copy would make the reported +#' outcome disagree with the value the recommendation was chosen by, i.e. two +#' wrong numbers instead of one. Nothing here changes a value. +#' +#' NOT placed inside get_outcome(), which is the single place the estimate is +#' produced, for two reasons that are each sufficient. First, outcome_type does +#' not reach it and neither do its callers carry one, so it cannot know that +#' [0, 1] is the right range and a guard there would fire on a continuous +#' outcome, whose range is not knowable. Second, and decisively, get_outcome() +#' is evaluated on the FLIPPED outcome scale under +#' outcome_goal_intention = "minimize", where the identity-link flip is a +#' negation: on a "minimize" run reporting a perfectly valid 0.0422, all +#' eleven values get_outcome() returned were negative. A guard there would have +#' to be threaded outcome_type AND lower_outcome_goal to say anything true, and +#' it is called once per grid point besides, so it could not warn once per run. +#' The check therefore belongs where the reported value exists on the caller's +#' own outcome scale and outcome_type is already in scope, which is +#' lago_optimization(), and it is called there once. +#' +#' The reported interval is mentioned but does not itself trigger the warning. +#' On the identity link get_confidence_set() deliberately does not confine the +#' interval, because the estimate it belongs to is unconfined and confining +#' only the interval would report an interval excluding its own estimate. That +#' decision and this warning are the same statement from two sides: neither +#' alters a number, and the user is told that what is reported is not a +#' probability. Triggering on a bound alone would re-open a decision already +#' taken, so the trigger is the estimate and the bounds are counted only to say +#' how far the report is affected. +#' +#' @param est_outcome A numeric value, the estimated outcome as reported, on +#' the caller's own outcome scale. +#' @param outcome_type A character string, "binary" or "continuous". Only a +#' binary outcome has a knowable range, so a continuous one returns at once. +#' @param link A character string, the link the outcome model was fitted on. +#' Used to name the mechanism, not to decide the condition: the condition is +#' whether the reported value is a probability, which is the property that is +#' violated. +#' @param reported_ci A numeric vector of the reported interval bounds at the +#' recommended intervention, or NULL when no confidence set was requested. +#' @param cs_rows The reported confidence set as a data.frame with +#' CI_lower_bound and CI_upper_bound columns, or NULL when there is none. +#' +#' @return Invisibly NULL. Called for its side effect of issuing one warning. +#' +#' @noRd +warn_if_outcome_outside_range <- function(est_outcome, + outcome_type, + link, + reported_ci = NULL, + cs_rows = NULL) { + # a continuous outcome's range is not knowable here, which is the same reason + # get_confidence_set() does not confine its interval. + if (!identical(outcome_type, "binary")) { + return(invisible(NULL)) + } + # a non-finite value is not a range violation and is reported, or refused, + # elsewhere. NULL and length-zero inputs fall out of this as FALSE. + outside_unit_range <- function(values) { + is.finite(values) & (values < 0 | values > 1) + } + if (!any(outside_unit_range(est_outcome))) { + return(invisible(NULL)) + } + + # how much of the rest of the report is affected, which is cheap to say and + # tells the user whether the headline number is the only one. Counted over + # the bounds as REPORTED, i.e. after rounding, so the count matches what the + # user can see rather than an unrounded value they cannot. + reported_bounds <- c( + reported_ci, + if (!is.null(cs_rows)) { + c(cs_rows$CI_lower_bound, cs_rows$CI_upper_bound) + } + ) + n_bounds_outside <- sum(outside_unit_range(reported_bounds)) + bounds_sentence <- if (n_bounds_outside > 0) { + paste0( + " ", n_bounds_outside, " reported confidence interval bound(s) are ", + "outside [0, 1] as well, and are likewise reported as computed: the ", + "interval is not confined on this link, because confining it around an ", + "unconfined estimate would report an interval excluding its own ", + "estimate." + ) + } else { + "" + } + + # the mechanism, named from the link rather than assumed, so the sentence is + # true for whichever link produced the value. Only "identity" can reach here + # in practice, since expit() cannot leave [0, 1], but a caller passing its + # own link should not be told about a model it did not fit. + mechanism <- if (link == "identity") { + paste0( + "The outcome model was fitted with link = \"identity\", so it is a ", + "linear probability model whose estimated outcome is the linear ", + "predictor itself and is not confined to [0, 1]. A fit whose every ", + "fitted value on the data is a probability still extrapolates outside ", + "[0, 1] at an intervention beyond the range its components were fitted ", + "over, which is what intervention bounds reaching past that range ask ", + "for." + ) + } else { + paste0( + "The outcome model was fitted with link = \"", link, "\", whose inverse ", + "did not confine the estimated outcome to [0, 1]." + ) + } + + warning(paste0( + "The estimated outcome is ", signif(est_outcome, 6), ", which is outside ", + "[0, 1] and so is not a probability, while the outcome is binary. ", + mechanism, + bounds_sentence, + "\nThe LAGO optimization still ran and the recommended intervention is ", + "the one the fitted model implies, so no reported value has been altered ", + "to fit the range. Please do not read the estimated outcome as a ", + "probability. Consider narrowing the intervention bounds to the range ", + "the data covers, or fitting the outcome model with link = \"logit\", ", + "whose estimated outcome is a probability by construction." + )) + + invisible(NULL) +} + + # The "minimize" direction is implemented by negating the fitted coefficients # (see lago_optimization()), which turns "reach an outcome at most as large as # the goal" into the maximization problem every optimizer here already solves. diff --git a/R/lago_optimization.R b/R/lago_optimization.R index 8055f13..8acf81c 100644 --- a/R/lago_optimization.R +++ b/R/lago_optimization.R @@ -553,6 +553,25 @@ lago_optimization <- function( NULL } + # a binary outcome's estimate is a probability, and this is the one place it + # exists as the value that will be REPORTED: on the caller's own outcome + # scale, with outcome_type in scope, once per run. Both matter. get_outcome() + # produces the estimate but is evaluated on the flipped scale under + # outcome_goal_intention = "minimize", where the identity-link flip is a + # negation, so a check there would fire on every value of a "minimize" run + # whose reported outcome is a perfectly good probability; and it is called + # once per grid point, so it could not warn once per run. outcome_type does + # not reach it either. Nothing here alters a value: the estimate is what the + # optimizers were driven by and what is reported, and clamping either would + # change which intervention is recommended or make the two disagree. + warn_if_outcome_outside_range( + est_outcome = est_outcome_goal, + outcome_type = outcome_type, + link = link, + reported_ci = est_outcome_ci, + cs_rows = if (include_confidence_set) cs$cs else NULL + ) + result <- if (!include_confidence_set) { c( list( diff --git a/man/get_confidence_set.Rd b/man/get_confidence_set.Rd index b605422..c067227 100644 --- a/man/get_confidence_set.Rd +++ b/man/get_confidence_set.Rd @@ -133,9 +133,12 @@ List( still the delta-method one and a bound at exactly 0 or 1 is one that has been truncated to the range. Not confined on the identity link, where the estimate is the linear predictor and is itself unbounded, so confining - the interval would report one that excludes its own estimate. Not - confined for a continuous outcome either, whose range is not knowable - here. NULL when that interval is not computable>, + the interval would report one that excludes its own estimate. That + estimate leaving [0, 1] is a defect in its own right and not this one; + lago_optimization() now warns when it does, so it is flagged rather than + silent, and the interval is still reported as computed for the same + reason as before. Not confined for a continuous outcome either, whose + range is not knowable here. NULL when that interval is not computable>, cs = = 1.3, ] + expect_gt(nrow(ok), 1) + best <- ok[which.min(ok$cost), ] + expect_identical(res$rec_int, c(best$x1, best$x2)) + # and the objective is NOT flat above the boundary, which is what clamping + # would have made it: the qualifying interventions have distinct estimates. + expect_gt(length(unique(ok$est)), 1) +}) From 92e89c2003f3b73bf9c53ab75bf78bdec39c0f6a Mon Sep 17 00:00:00 2001 From: Ante Bing Date: Fri, 7 Aug 2026 19:56:34 +0000 Subject: [PATCH 8/9] Silence a spurious warning and stop the range message contradicting itself Two points from review of the two diagnostics, both about the diagnostics themselves rather than any returned value. The covariate-support check called range() with na.rm on the observed values and then skipped a non-finite result, but an all-NA numeric column made range() itself warn about no non-missing arguments before that skip could run: a second, unrelated warning from something meant to emit one clear one. It now drops the non-finite values before range() and skips a column with none, so the base-R warning never fires. Not reachable through lago_optimization(), where the fit fails first, but reachable on a direct get_confidence_set() call. The out-of-range estimate message printed the value with signif() to six figures, so an estimate just past the boundary, 1.0000004, printed as 1 and the message read "the estimated outcome is 1, which is outside [0, 1]" -- a statement and its negation. It now prints at full precision, so the shown value is always one the reader can see is outside the range. The existing test pinned the rounded value and is updated to the fuller one, and a new test drives the message at a barely-out value to hold that it never rounds to look in range. --- R/get_confidence_set.R | 10 +++++++-- R/get_outcome.R | 8 ++++++- .../test-goal-modes-and-diagnostics.R | 22 ++++++++++++++++++- 3 files changed, 36 insertions(+), 4 deletions(-) diff --git a/R/get_confidence_set.R b/R/get_confidence_set.R index 1f2fad9..cdc1336 100644 --- a/R/get_confidence_set.R +++ b/R/get_confidence_set.R @@ -1435,10 +1435,16 @@ warn_covariates_held_off_support <- function(additional_covariates, if (!is.numeric(values)) { next } - observed <- range(values, na.rm = TRUE) - if (!all(is.finite(observed))) { + # finite values only, and taken before range() rather than with na.rm: an + # all-NA column makes range(na.rm = TRUE) itself warn ("no non-missing + # arguments to min") before the is.finite check below could skip it, which + # would be a second, unrelated warning from a diagnostic that is meant to + # emit one clear one. + finite_values <- values[is.finite(values)] + if (length(finite_values) == 0) { next } + observed <- range(finite_values) # 0 inside the observed range is a value the covariate takes, so holding # it there is not an extrapolation if (observed[1] <= 0 && observed[2] >= 0) { diff --git a/R/get_outcome.R b/R/get_outcome.R index db3e71b..78a7ef9 100644 --- a/R/get_outcome.R +++ b/R/get_outcome.R @@ -219,8 +219,14 @@ warn_if_outcome_outside_range <- function(est_outcome, ) } + # format() to enough digits that the printed value cannot round to a number + # inside [0, 1]: signif(1.0000004, 6) is 1, which would read "the estimated + # outcome is 1, which is outside [0, 1]", contradicting itself. The value is + # only just outside the range in that case, but the message must not say a + # thing and its negation. warning(paste0( - "The estimated outcome is ", signif(est_outcome, 6), ", which is outside ", + "The estimated outcome is ", format(est_outcome, digits = 15), + ", which is outside ", "[0, 1] and so is not a probability, while the outcome is binary. ", mechanism, bounds_sentence, diff --git a/tests/testthat/test-goal-modes-and-diagnostics.R b/tests/testthat/test-goal-modes-and-diagnostics.R index 75f8ff3..5249888 100644 --- a/tests/testthat/test-goal-modes-and-diagnostics.R +++ b/tests/testthat/test-goal-modes-and-diagnostics.R @@ -206,7 +206,11 @@ test_that("the range warning names the range, the cause and the bounds", { ws <- out_of_range_warnings() w <- ws[grepl("not a probability", ws)] expect_length(w, 1) - expect_match(w, "estimated outcome is 1\\.44073") + # the value is printed at full precision, not rounded: a value only just + # outside the range, say 1.0000004, would round to 1 and read "the estimated + # outcome is 1, which is outside [0, 1]", contradicting itself. So the printed + # value is asserted to begin with the true digits and not to be a rounded one. + expect_match(w, "estimated outcome is 1\\.4407272") expect_match(w, "outside \\[0, 1\\]") expect_match(w, "binary") expect_match(w, "linear probability model") @@ -223,6 +227,22 @@ test_that("the range warning names the range, the cause and the bounds", { expect_match(w, "reported confidence interval bound\\(s\\) are outside") }) +test_that("the range warning does not round a barely-out value to look in range", { + # an estimate just past the boundary, 1 + 4e-7, must not print as "1", which + # would make the message say "the estimated outcome is 1, which is outside + # [0, 1]" -- a statement and its negation. The message builder is called + # directly, since driving the optimizer to land an estimate exactly there is + # not controllable. + warn_fn <- getFromNamespace("warn_if_outcome_outside_range", "LAGO") + msg <- tryCatch( + warn_fn(1.0000004, "binary", "identity"), + warning = conditionMessage + ) + expect_true(is.character(msg)) + expect_no_match(msg, "outcome is 1, which is outside") + expect_match(msg, "1.0000004", fixed = TRUE) +}) + test_that("the range warning fires once per run, not once per grid point", { # the grid here is 10 x 3 = 30 interventions and the confidence set grid is # the same size, so a per-point warning would flood. Counting them is the From 823359e328e5171752e0484e538c29765e12f7d1 Mon Sep 17 00:00:00 2001 From: Ante Bing Date: Fri, 7 Aug 2026 20:06:55 +0000 Subject: [PATCH 9/9] Pin the closed boundary and the offenders-only list Two survivors from mutation review, both under-pinning rather than wrong code. The out-of-range check is strict: 0 and 1 are probabilities, so an estimate on either is in range and silent. Nothing pinned that, so loosening the test to <= or >= would warn on a legitimate probability, which a saturated fit reports, and the suite would stay green. An estimate at exactly 0 and exactly 1 is now asserted silent on both links, with just-outside each boundary asserted to warn so the silence is the closed boundary and not a dead check. The covariate warning lists the covariates held off their support. The multiple-offender test asserted the offenders are named but not that a non-offender is absent, so a message naming every covariate would have passed while telling the user a covariate held on its support is off it. A third covariate observed spanning 0 is added and asserted absent from the message. Both verified by mutation: loosening the boundary to non-strict fails the first, and dropping the spans-0 gate so every covariate is named fails the second. --- tests/testthat/test-covariate-support.R | 20 +++++++++++++++---- .../test-goal-modes-and-diagnostics.R | 17 ++++++++++++++++ 2 files changed, 33 insertions(+), 4 deletions(-) diff --git a/tests/testthat/test-covariate-support.R b/tests/testthat/test-covariate-support.R index 48dd7cc..d1093a5 100644 --- a/tests/testthat/test-covariate-support.R +++ b/tests/testthat/test-covariate-support.R @@ -75,19 +75,28 @@ test_that("the warning fires exactly once, listing every offending covariate", { x <- bb$launch_duration bb$staff_count <- 20 + (x - min(x)) / (max(x) - min(x)) * 10 expect_gt(min(bb$staff_count), 0) + # a THIRD numeric covariate observed spanning 0, which must NOT be named: it + # is not held off its support. Included so the test pins that the message + # lists the offenders and ONLY the offenders, not every covariate. + bb$balance <- (x - mean(x)) + expect_lt(min(bb$balance), 0) + expect_gt(max(bb$balance), 0) model <- suppressWarnings(glm( pp3_oxytocin_mother ~ coaching_updt + launch_duration + - clinic_size + staff_count, + clinic_size + staff_count + balance, data = bb, family = binomial() )) run <- function() { get_confidence_set( predictors_data = bb[ - , c("coaching_updt", "launch_duration", "clinic_size", "staff_count"), + , c( + "coaching_updt", "launch_duration", "clinic_size", "staff_count", + "balance" + ), drop = FALSE ], - additional_covariates = c("clinic_size", "staff_count"), + additional_covariates = c("clinic_size", "staff_count", "balance"), intervention_components = c("coaching_updt", "launch_duration"), outcome_data = bb$pp3_oxytocin_mother, fitted_model = model, @@ -115,9 +124,12 @@ test_that("the warning fires exactly once, listing every offending covariate", { value = TRUE ) expect_length(support_warnings, 1) - # and the single warning names both covariates + # and the single warning names both offenders and ONLY them: a message that + # listed every covariate would pass an offenders-are-named check while telling + # the user a covariate held on its support is off it. expect_match(support_warnings, "clinic_size") expect_match(support_warnings, "staff_count") + expect_no_match(support_warnings, "balance") }) diff --git a/tests/testthat/test-goal-modes-and-diagnostics.R b/tests/testthat/test-goal-modes-and-diagnostics.R index 5249888..86ee0a0 100644 --- a/tests/testthat/test-goal-modes-and-diagnostics.R +++ b/tests/testthat/test-goal-modes-and-diagnostics.R @@ -243,6 +243,23 @@ test_that("the range warning does not round a barely-out value to look in range" expect_match(msg, "1.0000004", fixed = TRUE) }) +test_that("an estimate at exactly 0 or 1 is a probability and does not warn", { + # the boundary is closed: 0 and 1 are probabilities, so an estimate landing on + # either is in range and must be silent. The range test is strict (< 0, > 1), + # and loosening it to <= or >= would warn on a legitimate probability, which a + # saturated fit can report. Called directly, since landing the estimate + # exactly on a boundary through the optimizer is not controllable. + warn_fn <- getFromNamespace("warn_if_outcome_outside_range", "LAGO") + for (boundary in c(0, 1)) { + expect_no_warning(warn_fn(boundary, "binary", "identity")) + expect_no_warning(warn_fn(boundary, "binary", "logit")) + } + # just outside each boundary does warn, so the silence above is the closed + # boundary and not a dead check + expect_warning(warn_fn(-1e-6, "binary", "identity"), "outside") + expect_warning(warn_fn(1 + 1e-6, "binary", "identity"), "outside") +}) + test_that("the range warning fires once per run, not once per grid point", { # the grid here is 10 x 3 = 30 interventions and the confidence set grid is # the same size, so a per-point warning would flood. Counting them is the