Fix/confidence bound range

R/get_outcome.R

@@ -82,6 +82,166 @@ 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]."
+    )
+  }
+
+  # 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 ", format(est_outcome, digits = 15),
+    ", 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.

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(

R/validate_inputs.R

@@ -359,13 +359,11 @@ validate_inputs <- function(
     # different things: see refuse_invalid_center_weights() for why each check
     # is there.
     refuse_invalid_center_weights(center_weights_for_outcome_goal)
+    # and whether the vector is a set of weights at all as a SET, i.e. sums to
+    # 1. Shared with the exported get_confidence_set() for the same reason the
+    # checks above are: see refuse_non_unit_weight_sum().
+    refuse_non_unit_weight_sum(center_weights_for_outcome_goal)
     weights_sum <- sum(center_weights_for_outcome_goal)
-    if (abs(weights_sum - 1) >= 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)
+}

tests/testthat/_snaps/presentation.md

@@ -490,7 +490,7 @@
       dose: 1.6476
       Cost: 1.6476
       Estimated outcome: 0.6
-      95% CI for the estimated outcome: 0.195 - 1.005
+      95% CI for the estimated outcome: 0.195 - 1
       Outcome goal: 0.6
 
       -- Confidence set 
@@ -499,7 +499,7 @@
       First rows of the confidence set (use $cs for all):
     Output
         dose CI_lower_bound CI_upper_bound cost
-      4    2          0.464          1.132    2
+      4    2          0.464              1    2
 
 ---
 
@@ -555,7 +555,7 @@
       dose: 1.6476
       Cost: 1.6476
       Estimated outcome: 0.6
-      95% CI for the estimated outcome: 0.195 - 1.005
+      95% CI for the estimated outcome: 0.195 - 1
       Outcome goal: 0.6
 
       -- Confidence set 
@@ -564,7 +564,7 @@
       First rows of the confidence set (use $cs for all):
     Output
         dose CI_lower_bound CI_upper_bound cost
-      4    2          0.464          1.132    2
+      4    2          0.464              1    2
 
 # print/summary snapshot: object without a confidence set