diff --git a/R/get_confidence_set.R b/R/get_confidence_set.R index 06c8d94..cdc1336 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,18 @@ #' 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. 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, @@ -561,6 +598,75 @@ 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 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. + # + # 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 + # exactly 0 or 1 means the delta-method bound was outside the range and has + # been truncated to it, not that the interval ends there. + # + # The alternative was to build the interval on the logit scale and map it + # back through expit(), which cannot leave [0, 1] at all and is the textbook + # fix. Rejected here: it is a different interval, so it would move EVERY + # reported binary number, including every one that is already a probability + # and already right, and the binary branch reports the delta-method interval + # deliberately -- it is what glm()'s own binary machinery reports and what + # this package has always reported. The comment on the continuous branch + # below records that choice. Changing the interval is a larger decision than + # this defect, which is that a number reported as a probability is not one; + # 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, 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. + # + # 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(). @@ -887,17 +993,36 @@ get_confidence_set <- function( if (link == "identity") { # For identity link, predictions are on the correct scale already. + # + # NOT clamped to [0, 1], unlike the binary branch above, because for a + # continuous outcome there is no range to clamp to. The outcome is a mean + # on whatever scale the caller's data is on, and nothing here knows what + # that scale is: "continuous" with a proportion in [0, 1] is supported and + # is what BB_proportions is, but so is a count, a duration or a + # difference, and the goal is only numeric. So an out-of-[0, 1] bound + # here is not knowably wrong -- for BB_proportions it would be, and for a + # blood pressure it plainly would not -- and clamping on the chance the + # outcome happens to be a proportion would corrupt every other continuous + # outcome. The binary branch can clamp precisely because "binary" fixes + # the range. ci_prob_all <- cbind(lb_prob_all, ub_prob_all) } else { # For "logit" (or other logistic-like) link, apply expit. + # This one is already confined to [0, 1] without any clamp, since it + # transforms the BOUNDS and expit() maps onto (0, 1). That is the same + # 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, "\"." @@ -926,8 +1051,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 @@ -999,8 +1124,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, @@ -1264,6 +1389,87 @@ 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 + } + # 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) { + 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..78a7ef9 100644 --- a/R/get_outcome.R +++ b/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. 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/R/validate_inputs.R b/R/validate_inputs.R index d9cca70..1d50e07 100644 --- a/R/validate_inputs.R +++ b/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) +} diff --git a/man/get_confidence_set.Rd b/man/get_confidence_set.Rd index 0ee4288..c067227 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,18 @@ List( 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. 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 & 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) + # 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")) { + fit <- call_cs(link) + model <- fit$model + res <- fit$res + + # 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) + 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. + # + # 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) + ) + 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 + # 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) + ) + + # 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 + # 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. 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 + # 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] + ) + + # 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] + ) + } + } +}) + + +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)) +}) diff --git a/tests/testthat/test-covariate-support.R b/tests/testthat/test-covariate-support.R new file mode 100644 index 0000000..d1093a5 --- /dev/null +++ b/tests/testthat/test-covariate-support.R @@ -0,0 +1,253 @@ +# Regression tests for the additional-covariate support warning. +# +# get_confidence_set() holds every additional covariate at 0 to build the +# prediction grid. For a NUMERIC covariate whose observed values never reach 0 +# -- a clinic size in [5, 6], say -- 0 is a value that never occurs in the +# data, so the reported estimated outcome and its interval are read at an +# unobserved covariate value: an extrapolation. The function now warns about +# that, naming the covariate and its observed range. +# +# The warning is DIAGNOSTIC ONLY. It changes no returned value: the covariate +# is still held at 0 either way, so the interval it warns about is exactly the +# interval it would have returned in silence. The tests below pin both halves: +# that it fires in the one case it should, stays silent in the cases it should +# not, and that the returned result is identical whether or not it fires. + + +# A numeric covariate observed in [5, 6], attached to the bundled BB_data. The +# range is built by rescaling an existing column into [5, 6] deterministically, +# so it never touches 0 and does not depend on a random seed. +bb_with_offset_covariate <- function() { + bb <- as.data.frame(BB_data) + x <- bb$coaching_updt + bb$clinic_size <- 5 + (x - min(x)) / (max(x) - min(x)) + bb +} + +# get_confidence_set() over BB_data with one additional covariate, at a coarse +# grid to keep it cheap. The model is fitted on exactly the predictors passed. +covariate_confidence_set <- function(bb, covariate) { + model <- suppressWarnings(glm( + as.formula(paste( + "pp3_oxytocin_mother ~ coaching_updt + launch_duration +", covariate + )), + data = bb, family = binomial() + )) + get_confidence_set( + predictors_data = bb[ + , c("coaching_updt", "launch_duration", covariate), + drop = FALSE + ], + additional_covariates = covariate, + intervention_components = c("coaching_updt", "launch_duration"), + outcome_data = bb$pp3_oxytocin_mother, + fitted_model = model, + link = "logit", + outcome_goal = 0.85, + outcome_type = "binary", + intervention_lower_bounds = c(1, 1), + intervention_upper_bounds = c(40, 5), + confidence_set_grid_step_size = c(10, 2), + cost_list_of_vectors = list(c(0, 1.7), c(0, 8)), + rec_int = c(1, 2.77847) + ) +} + + +test_that("a numeric covariate observed away from 0 warns it is extrapolated", { + bb <- bb_with_offset_covariate() + # precondition: the covariate really is observed in [5, 6], excluding 0 + expect_gt(min(bb$clinic_size), 0) + expect_equal(range(bb$clinic_size), c(5, 6)) + + # the warning names the covariate and its observed range, and says the report + # is an extrapolation because the covariate is held at 0 off its support + expect_warning( + covariate_confidence_set(bb, "clinic_size"), + "clinic_size.*\\[5, 6\\].*held at 0.*extrapolation" + ) +}) + + +test_that("the warning fires exactly once, listing every offending covariate", { + bb <- bb_with_offset_covariate() + # a second numeric covariate also observed away from 0, in a different range + 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 + balance, + data = bb, family = binomial() + )) + run <- function() { + get_confidence_set( + predictors_data = bb[ + , c( + "coaching_updt", "launch_duration", "clinic_size", "staff_count", + "balance" + ), + drop = FALSE + ], + additional_covariates = c("clinic_size", "staff_count", "balance"), + intervention_components = c("coaching_updt", "launch_duration"), + outcome_data = bb$pp3_oxytocin_mother, + fitted_model = model, + link = "logit", + outcome_goal = 0.85, + outcome_type = "binary", + intervention_lower_bounds = c(1, 1), + intervention_upper_bounds = c(40, 5), + confidence_set_grid_step_size = c(10, 2), + cost_list_of_vectors = list(c(0, 1.7), c(0, 8)), + rec_int = c(1, 2.77847) + ) + } + # the grid has many rows; the warning is fired once for the call, not per row + warnings <- character(0) + withCallingHandlers( + run(), + warning = function(w) { + warnings <<- c(warnings, conditionMessage(w)) + invokeRestart("muffleWarning") + } + ) + support_warnings <- grep( + "held at 0 to compute the confidence set", warnings, + value = TRUE + ) + expect_length(support_warnings, 1) + # 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") +}) + + +test_that("a numeric covariate whose range includes 0 does not warn", { + bb <- bb_with_offset_covariate() + # centre the covariate so its observed range straddles 0; 0 is then a value + # it takes in the data, so holding it there is not an extrapolation + bb$clinic_size <- bb$clinic_size - 5.5 + expect_lt(min(bb$clinic_size), 0) + expect_gt(max(bb$clinic_size), 0) + + expect_no_warning(covariate_confidence_set(bb, "clinic_size")) +}) + + +test_that("a numeric covariate whose range touches 0 does not warn", { + bb <- bb_with_offset_covariate() + # 0 sitting at the boundary of the observed range is still an observed value, + # so the boundary case does not warn either + bb$clinic_size <- bb$clinic_size - 5 + expect_equal(min(bb$clinic_size), 0) + + expect_no_warning(covariate_confidence_set(bb, "clinic_size")) +}) + + +test_that("a factor additional covariate does not warn", { + bb <- as.data.frame(BB_data) + # held at 0, a factor covariate sits at its reference level, an OBSERVED + # level, so it is not an extrapolation and must not warn + bb$arm <- factor( + ifelse(bb$pre_post == 0, "pre", "post"), + levels = c("post", "pre") + ) + expect_no_warning(covariate_confidence_set(bb, "arm")) +}) + + +test_that("a character additional covariate does not warn", { + bb <- as.data.frame(BB_data) + # a character covariate is contrast-coded like a factor: held at 0 it is its + # reference level, so it is not extrapolated + bb$arm <- as.character(ifelse(bb$pre_post == 0, "pre", "post")) + expect_no_warning(covariate_confidence_set(bb, "arm")) +}) + + +test_that("no additional covariates means no warning", { + bb <- as.data.frame(BB_data) + model <- suppressWarnings(glm( + pp3_oxytocin_mother ~ coaching_updt + launch_duration, + data = bb, family = binomial() + )) + run <- function() { + get_confidence_set( + predictors_data = bb[ + , c("coaching_updt", "launch_duration"), + drop = FALSE + ], + intervention_components = c("coaching_updt", "launch_duration"), + outcome_data = bb$pp3_oxytocin_mother, + fitted_model = model, + link = "logit", + outcome_goal = 0.85, + outcome_type = "binary", + intervention_lower_bounds = c(1, 1), + intervention_upper_bounds = c(40, 5), + confidence_set_grid_step_size = c(10, 2), + cost_list_of_vectors = list(c(0, 1.7), c(0, 8)), + rec_int = c(1, 2.77847) + ) + } + expect_no_warning(run()) +}) + + +test_that("the warning changes no returned value (it is diagnostic only)", { + # the warning reads predictors_data only to decide whether to fire. The + # binary branch's interval is built from vcov(fitted_model), NOT from + # predictors_data, so the SAME fitted model with a predictors_data whose + # covariate range straddles 0 rather than excluding it toggles the warning + # off while returning the identical object. That isolates the warning as a + # pure side effect: the only thing that changed is whether the condition was + # signalled, and the returned result is byte-for-byte the same. + bb <- bb_with_offset_covariate() + model <- suppressWarnings(glm( + pp3_oxytocin_mother ~ coaching_updt + launch_duration + clinic_size, + data = bb, family = binomial() + )) + run <- function(predictors_data) { + get_confidence_set( + predictors_data = predictors_data, + additional_covariates = "clinic_size", + intervention_components = c("coaching_updt", "launch_duration"), + outcome_data = bb$pp3_oxytocin_mother, + fitted_model = model, + link = "logit", + outcome_goal = 0.85, + outcome_type = "binary", + intervention_lower_bounds = c(1, 1), + intervention_upper_bounds = c(40, 5), + confidence_set_grid_step_size = c(10, 2), + cost_list_of_vectors = list(c(0, 1.7), c(0, 8)), + rec_int = c(1, 2.77847) + ) + } + pd_excludes_zero <- bb[ + , c("coaching_updt", "launch_duration", "clinic_size"), + drop = FALSE + ] + pd_includes_zero <- pd_excludes_zero + pd_includes_zero$clinic_size <- pd_includes_zero$clinic_size - 5.5 + + with_warning <- suppressWarnings(run(pd_excludes_zero)) + no_warning <- run(pd_includes_zero) + + # the second call really is the one that does not warn + expect_no_warning(run(pd_includes_zero)) + # and the returned object is identical whether or not the warning fired + expect_identical(with_warning, no_warning) +}) diff --git a/tests/testthat/test-goal-modes-and-diagnostics.R b/tests/testthat/test-goal-modes-and-diagnostics.R index 8c6f304..86ee0a0 100644 --- a/tests/testthat/test-goal-modes-and-diagnostics.R +++ b/tests/testthat/test-goal-modes-and-diagnostics.R @@ -106,3 +106,254 @@ test_that("providing a grid step size switches method to grid_search", { "grid_search" ) }) + + +# A binary outcome's ESTIMATE reported outside [0, 1]. On the identity link the +# outcome model is a linear probability model, so the estimate is the linear +# predictor and is not confined to anything: intervention bounds reaching past +# the range the components were fitted over extrapolate it out of the range, +# and glm() does not object because every fitted value ON THE DATA is a +# probability. It was reported with no error and no warning. +# +# The fixture is built so the mechanism is the ONLY thing being tested: x1 is +# fitted over 0..4 and the upper bound is 9, so nothing but the extrapolation +# puts the estimate out of range. max(fitted()) < 1 is asserted, so a run that +# reproduces this cannot be dismissed as a fit glm() would already have +# complained about. +# Both components are significant on this seed, deliberately: the fit +# diagnostics warn about a non-significant intervention component, and a +# fixture that tripped that too would emit a second warning that expect_warning +# lets escape and the suite then counts as a stray. The seed is chosen so the +# range warning is the ONLY warning the run emits. +out_of_range_fixture <- function() { + set.seed(152) + n <- 60 + x1 <- rep(0:4, each = 12) + x2 <- rep(1:3, times = 20) + p <- pmin(pmax(0.10 + 0.11 * x1 + 0.10 * x2, 0), 1) + data.frame(y = rbinom(n, 1, p), x1 = x1, x2 = x2) +} + +out_of_range_run <- function(...) { + args <- list( + data = out_of_range_fixture(), + outcome_name = "y", + outcome_type = "binary", + glm_family = "binomial", + link = "identity", + intervention_components = c("x1", "x2"), + intervention_lower_bounds = c(0, 1), + # 9 is well beyond the 0..4 x1 was fitted over, which is what extrapolates + intervention_upper_bounds = c(9, 3), + cost_list_of_vectors = list(c(0, 1), c(0, 1)), + outcome_goal = 1.3, + optimization_method = "grid_search", + optimization_grid_search_step_size = c(1, 1), + confidence_set_grid_step_size = c(1, 1), + quiet = TRUE + ) + ov <- list(...) + for (nm in names(ov)) args[nm] <- list(ov[[nm]]) + do.call(lago_optimization, args) +} + +test_that("a binary estimate outside [0, 1] warns and is not altered", { + d <- out_of_range_fixture() + model <- suppressWarnings(glm( + y ~ x1 + x2, data = d, family = binomial(link = "identity") + )) + # the preconditions, so this cannot go vacuous: glm() converged and every + # fitted value on the DATA is a probability, so the fit is one glm() accepts + # without complaint and only the extrapolated prediction leaves the range. + expect_true(model$converged) + expect_lt(max(fitted(model)), 1) + + expect_warning( + res <- out_of_range_run(), + "outside \\[0, 1\\] and so is not a probability" + ) + + # the estimate is out of range, which is the defect, and it is reported AS + # COMPUTED: nothing is clamped. Held against the model's own linear + # predictor at the recommended intervention, computed here from coef() + # alone, so this is not the package agreeing with itself. + expect_gt(res$est_outcome_goal, 1) + hand <- as.numeric(c(1, res$rec_int) %*% coef(model)) + expect_equal(res$est_outcome_goal, hand, tolerance = 1e-12) +}) + +# every warning one run emits, so a test can pick out the one it is about. The +# fixture also trips the pre-existing significance warning on x2, and that one +# comes FIRST: tryCatch(warning=) would hand back that message instead, which +# is what made the assertions below look like failures against a warning they +# were never about. +out_of_range_warnings <- function(...) { + ws <- character(0) + withCallingHandlers( + out_of_range_run(...), + warning = function(w) { + ws <<- c(ws, conditionMessage(w)) + invokeRestart("muffleWarning") + } + ) + ws +} + +test_that("the range warning names the range, the cause and the bounds", { + # actionable, not merely present: it has to say what is out of range, what + # the range is, why it happened, and that nothing was altered. A warning + # that only says "out of range" leaves the user with no next step. + ws <- out_of_range_warnings() + w <- ws[grepl("not a probability", ws)] + expect_length(w, 1) + # 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") + expect_match(w, "extrapolates outside") + expect_match(w, "link = \"identity\"") + # the mechanism, named as extrapolation beyond the fitted range + expect_match(w, "beyond the range its components were fitted") + # that nothing was altered, which is the whole reason this warns + expect_match(w, "no reported value has been altered") + # and the actionable alternatives + expect_match(w, "narrowing the intervention bounds") + expect_match(w, "link = \"logit\"") + # it counts the affected reported bounds rather than staying silent on them + 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("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 + # assertion. + ws <- out_of_range_warnings() + expect_identical(sum(grepl("not a probability", ws)), 1L) +}) + +test_that("the range warning does not fire where the estimate is in range", { + # the same fixture with bounds INSIDE the fitted range: the estimate is a + # probability, so there is nothing to warn about and the run must be silent + # on this account. + res <- suppressWarnings(out_of_range_run( + intervention_upper_bounds = c(4, 3), outcome_goal = 0.6 + )) + ws <- out_of_range_warnings( + intervention_upper_bounds = c(4, 3), outcome_goal = 0.6 + ) + expect_lte(res$est_outcome_goal, 1) + expect_gte(res$est_outcome_goal, 0) + expect_identical(sum(grepl("not a probability", ws)), 0L) +}) + +test_that("a minimize run whose estimate is a probability does not warn", { + # THE case that decides where this check belongs. get_outcome() is evaluated + # on the FLIPPED outcome scale under "minimize", where the identity-link flip + # is a negation, so every value it returns on such a run is negative even + # when the reported estimate is a perfectly good probability. A guard placed + # inside get_outcome() would fire on all of them. This asserts the check sees + # the value as REPORTED, on the caller's own scale. + min_args <- list( + intervention_upper_bounds = c(4, 3), + outcome_goal = 0.2, + outcome_goal_intention = "minimize", + include_confidence_set = FALSE + ) + res <- suppressWarnings(do.call(out_of_range_run, min_args)) + ws <- do.call(out_of_range_warnings, min_args) + expect_gte(res$est_outcome_goal, 0) + expect_lte(res$est_outcome_goal, 1) + expect_identical(sum(grepl("not a probability", ws)), 0L) +}) + +test_that("a continuous outcome outside [0, 1] does not warn", { + # a continuous outcome's range is not knowable here, which is the same reason + # get_confidence_set() does not confine its interval. mpg is far outside + # [0, 1] and must pass without this warning, which is what a check keyed on + # the link alone, or one placed where outcome_type does not reach, would get + # wrong. + ws <- character(0) + res <- withCallingHandlers( + suppressMessages(lago_optimization( + data = mtcars, outcome_name = "mpg", outcome_type = "continuous", + glm_family = "gaussian", link = "identity", + intervention_components = c("gear", "qsec"), + intervention_lower_bounds = c(0, 0), + intervention_upper_bounds = c(10, 350), + cost_list_of_vectors = list(c(0, 4), c(4, 6)), + outcome_goal = 40, outcome_goal_intention = "maximize", + include_confidence_set = FALSE, quiet = TRUE + )), + warning = function(w) { + ws <<- c(ws, conditionMessage(w)) + invokeRestart("muffleWarning") + } + ) + expect_gt(res$est_outcome_goal, 1) + expect_identical(sum(grepl("not a probability", ws)), 0L) +}) + +test_that("warning or not, the recommendation is the one the model implies", { + # the reason this WARNS rather than clamps. get_outcome() drives every + # optimizer and the goal comparison, so clamping the estimate would change + # which intervention is recommended. Here the recommendation is held against + # the cheapest grid intervention whose own linear predictor meets the goal, + # computed from coef() alone: that is what the optimizer is supposed to + # return, and a clamped objective could not, since every intervention above + # the boundary would tie at 1. + d <- out_of_range_fixture() + model <- suppressWarnings(glm( + y ~ x1 + x2, data = d, family = binomial(link = "identity") + )) + res <- suppressWarnings(out_of_range_run()) + + grid <- expand.grid(x1 = seq(0, 9, by = 1), x2 = seq(1, 3, by = 1)) + grid$est <- apply(grid, 1, function(r) { + as.numeric(c(1, r[["x1"]], r[["x2"]]) %*% coef(model)) + }) + grid$cost <- grid$x1 + grid$x2 + ok <- grid[grid$est >= 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) +})