From e046c95912d37599e929c320720c811cb206c65a Mon Sep 17 00:00:00 2001 From: Ante Bing Date: Thu, 6 Aug 2026 22:30:28 +0000 Subject: [PATCH 1/4] Guard the exported entry point, and build a binary interval on its own link Three defects the reviews of the previous change disclosed and left. The message refusing an all-failed numerical optimization recommended grid search and blamed having more than three intervention components. Neither is right for the case that actually reaches it. A model whose terms are not all estimable returns NA coefficients, every outcome computed from them is NA, and every restart fails, so the path is reached through rank deficiency rather than through a hard optimization, and grid search fails on the same input for the same reason. The refusal now names the terms that could not be estimated and says to drop or combine them, keeping the old wording only for the case it was written for. Grid search no longer walks into its own unguarded comparison on the same input: it refuses it with the same message, which fires exactly when the comparison it precedes would have failed. get_confidence_set() is exported and does not go through validate_inputs(), so the weight checks added last time did not protect it. A negative center weight produced an upper confidence bound above 1 for a binary outcome, silently. The checks are now one function that both entry points call, so the two cannot disagree, and the argument documentation that already claimed weights must be non-negative and finite is true for both. The sum-to-one check stays out of the exported path deliberately: validate_inputs() renormalises, and the confidence set has nothing to renormalise into. A binary outcome's interval was built on the logit scale whatever link the model was fitted on, because that branch keyed on the outcome type and applied expit() and the logit delta-method factor unconditionally. A binomial model fitted with an identity link therefore reported an interval on the wrong scale: for one such fit 0.681 to 0.685, where the interval on the model's own scale is 0.759 to 0.779, which does not even contain the estimate. The point estimate was always right; only the interval was wrong. The branch now follows the link, and the logit arithmetic is unchanged term for term, since reassociating it moves the result in the last few digits. --- R/get_confidence_set.R | 74 ++++- R/get_recommended_interventions.R | 162 +++++++++- R/rec_int_processor.R | 16 +- R/validate_inputs.R | 141 ++++++--- man/get_confidence_set.Rd | 4 +- man/get_recommended_interventions.Rd | 12 +- tests/testthat/test-confidence-set-contract.R | 282 ++++++++++++++++++ tests/testthat/test-optimization.R | 75 +++-- tests/testthat/test-outcome-internals.R | 207 ++++++++++++- tests/testthat/test-validate-inputs.R | 58 ++++ 10 files changed, 936 insertions(+), 95 deletions(-) diff --git a/R/get_confidence_set.R b/R/get_confidence_set.R index 4db99bb..06c8d94 100644 --- a/R/get_confidence_set.R +++ b/R/get_confidence_set.R @@ -10,7 +10,9 @@ #' @param center_weights_for_outcome_goal A numeric vector. Specifies the #' 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. +#' 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. #' @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 @@ -185,6 +187,19 @@ get_confidence_set <- function( if (!link %in% supported_outcome_links()) { stop(unsupported_link_message(link)) } + # and for the same reason the weights are checked here too, whenever they are + # used at all. They average the per-center outcomes, so a negative one puts + # the interval outside the range of the intervals it averages: a weight of + # -8 beside 8.5 and 0.5 sums to 1, passed every check this function made, and + # reported a CI_upper_bound of 1.014 for a BINARY outcome, i.e. a + # "probability" above 1. lago_optimization() refuses that in + # validate_inputs(), which this function does not go through, so the guard + # would not otherwise cover a direct caller -- and the @param text saying the + # weights must be non-negative and finite would be true of one entry point + # and not the other. The same guard, so both name the same reason. + if (include_center_effects) { + refuse_invalid_center_weights(center_weights_for_outcome_goal) + } # Create a list to store sequences for each component sequences <- list() # Generate sequences for each intervention component @@ -505,10 +520,40 @@ get_confidence_set <- function( drop = FALSE ] new_data <- as.matrix(new_data) - pred_all <- expit(new_data %*% model_coefs) - se_pred_all <- sqrt( + # the linear predictor and its standard error, which are what the model is + # fitted on whatever the link, and are mapped onto the outcome scale by the + # link below. Computed once so the two link branches can differ ONLY in + # that map, the way get_outcome() is written. + linear_predictor <- new_data %*% model_coefs + se_linear_predictor <- sqrt( diag((new_data) %*% model_vcov %*% t(new_data)) - ) * pred_all * (1 - pred_all) + ) + + # link is either "logit" or "identity", the only links the outcome + # machinery implements, see supported_outcome_links(). Both this and the + # continuous branch below key on it, because the scale the interval belongs + # on is a property of the LINK the model was fitted on and not of the + # outcome's type. This branch used to apply expit() and the logit + # delta-method factor unconditionally, keyed on outcome_type alone, so a + # binomial fit with link = "identity" -- which lago_optimization() accepts + # -- was reported on the logit scale: an interval of 0.636 to 0.655 where + # the identity-scale interval is 1.040 to 1.072, not even containing the + # point estimate the same run reported. The point estimate was always right, + # since get_outcome() has keyed on the link all along; only the interval + # was on the wrong scale. + if (link == "identity") { + # the outcome IS the linear predictor, so the interval is already on the + # outcome scale and the delta-method factor is the derivative of the + # identity map, i.e. 1. + pred_all <- linear_predictor + se_pred_all <- se_linear_predictor + } else { + # the outcome is expit() of the linear predictor, and the delta method + # carries the standard error across by the derivative of that map, + # d expit(eta) / d eta = p * (1 - p). + pred_all <- expit(linear_predictor) + se_pred_all <- se_linear_predictor * pred_all * (1 - pred_all) + } # lower and upper bounds of predictions lb_prob_all <- pred_all - critical_value * se_pred_all @@ -521,6 +566,27 @@ get_confidence_set <- function( # machinery implements, see supported_outcome_links(). # If link == "logit", use the logistic-like approach # If link == "identity", use linear approach. + # + # This branch was already keyed on the link, which is why only the binary + # branch above needed correcting. The two are NOT merged into one, and the + # difference is deliberate rather than drift: this branch builds its own + # variance-covariance matrix from predictors_data -- sandwich-form, + # optionally clustered on up to two dimensions -- because a continuous + # outcome's clustering is not what glm()'s vcov() assumes, while the binary + # branch takes vcov(fitted_model) directly. So the two differ in the + # VARIANCE they use, and only agree on how a variance becomes an interval. + # What they must not differ on is the SCALE, which is the link's business, + # and that is now the same decision on both sides: transform the bounds for + # "logit", leave them alone for "identity". + # + # The order of the two steps differs too, and both are correct. Here the + # interval is built on the linear-predictor scale and its BOUNDS are + # transformed, which cannot leave [0, 1] because expit() cannot. The binary + # branch transforms the POINT and carries the standard error across by the + # delta method, which is the interval glm()'s own binary machinery reports + # and is what this package has always reported for a binary outcome. + # Changing that would move every binary logit number, which is not what is + # wrong with either branch. # define the function to manually calculate var-cov matrix get_vcov <- function(predictors_data, diff --git a/R/get_recommended_interventions.R b/R/get_recommended_interventions.R index 854e606..eaa5088 100644 --- a/R/get_recommended_interventions.R +++ b/R/get_recommended_interventions.R @@ -30,6 +30,55 @@ unachievable_goal_message <- function(lower_outcome_goal) { } +#' rank_deficient_outcome_message +#' +#' @description The message for an optimization that had nothing to compare +#' because the OUTCOME MODEL could not be estimated, i.e. glm() returned NA for +#' at least one coefficient. +#' +#' @details This is the cause of the no-computable-outcome situation, as +#' distinct from the optimizer struggling with a problem it could in principle +#' solve, and the two need opposite advice. An aliased term's coefficient is +#' NA, the outcome computed from it is NA, and nothing that searches over +#' interventions can recover from that: every point it tries is NA whatever the +#' search. So this branch must NOT send the caller to the other optimization +#' method, which is what the single message here used to do for every cause. On +#' the only configuration known to reach the all-failed path at all, following +#' that advice reached the grid search's own goal comparison as the base-R error +#' "missing value where TRUE/FALSE needed", i.e. a second opaque failure after +#' being told to expect a solution. +#' +#' What the caller can act on is the model, so the message names the terms that +#' could not be estimated and says to drop or combine them. +#' +#' The condition is read off the fitted coefficients rather than off the NA +#' outcomes, so the terms can be NAMED: by the time an outcome is NA the NA has +#' been summed into the center-level effects and carries their names, not the +#' aliased term's. +#' +#' @param aliased_coef_names A character vector, the names of the coefficients +#' glm() could not estimate. Non-empty: the caller branches on that. +#' +#' @return A character string, the message to raise. +#' +#' @noRd +rank_deficient_outcome_message <- function(aliased_coef_names) { + paste( + "No outcome could be estimated at any of the interventions tried,", + "because the outcome model is rank-deficient: the coefficient(s)", + paste(aliased_coef_names, collapse = ", "), + "could not be estimated (glm() returned NA for them), which happens when", + "those predictors are collinear with others in the model. An outcome", + "computed from an NA coefficient is NA, so no optimization method can", + "proceed on this fit and changing the 'optimization_method' will not help.", + "Drop or combine the collinear predictor(s) and fit again. Common causes", + "are two intervention components that are rescalings of one another, and", + "fixed center and time effects whose assignments are confounded rather", + "than crossed." + ) +} + + #' refuse_if_all_restarts_failed #' #' @description Internal guard for a multi-start numerical optimization in @@ -40,8 +89,17 @@ unachievable_goal_message <- function(lower_outcome_goal) { #' NlcOptim::solnl() call as NA and carry on, so that one bad starting point #' does not lose the whole optimization. All of them failing is a different #' situation: there is nothing to select, and the only thing to tell the -#' caller is that this optimizer could not solve their problem and which one -#' can. +#' caller is what stopped the optimization and what to do about it. +#' +#' WHICH of those to say depends on the cause, and the two causes want opposite +#' advice, so the message branches on aliased_coef_names. A rank-deficient +#' outcome model gets rank_deficient_outcome_message(), which does not +#' recommend another optimizer, because none can help: every outcome is NA +#' before any search begins. The wording below is kept for the other case, a +#' genuine numerical-optimization failure on an estimable model, which is what +#' a caller in that position still needs. One message for both blamed a cause +#' ("more than three intervention components") that no reachable configuration +#' demonstrated, and recommended a method that fails on the same input. #' #' The condition and the message live here, once, rather than at each loop. #' The max-achievable-outcome loop had no such check at all while the cost loop @@ -54,14 +112,23 @@ unachievable_goal_message <- function(lower_outcome_goal) { #' #' @param results A numeric vector, one entry per restart, holding the value #' that restart converged to. NA marks a restart whose optimization failed. +#' @param aliased_coef_names A character vector, the names of the outcome +#' model's coefficients glm() could not estimate. Empty (the default) means the +#' model is of full rank, so an all-failed restart set is the optimizer's own +#' failure. Defaulted so that a call site with no model to hand still gets the +#' guard, with the wording it had. #' #' @return Invisibly NULL when at least one restart succeeded. Raises #' otherwise, so the callers below can treat returning as "there is something #' to choose from". #' #' @noRd -refuse_if_all_restarts_failed <- function(results) { +refuse_if_all_restarts_failed <- function(results, + aliased_coef_names = character(0)) { if (all(is.na(results))) { + if (length(aliased_coef_names) > 0) { + stop(rank_deficient_outcome_message(aliased_coef_names)) + } stop(paste( "Numerical optimization failed to find a solution.", "Please consider using the 'grid_search' method by", @@ -76,6 +143,59 @@ refuse_if_all_restarts_failed <- function(results) { } +#' refuse_if_no_grid_outcome +#' +#' @description Internal guard for a grid search none of whose grid points has +#' a computable outcome, so there is nothing to compare against the goal. +#' +#' @details The grid search's counterpart to +#' refuse_if_all_restarts_failed(). The numerical optimizer records a failed +#' restart as NA and refuses the all-failed case by name; the grid search had no +#' such check, and an NA outcome flowed into `max(all_outcomes) >= +#' new_outcome_goal`, which is not FALSE for an NA but the R error "missing +#' value where TRUE/FALSE needed". +#' +#' Both guards exist for the same reason and both distinguish the same two +#' causes, so the wording is shared: a rank-deficient outcome model must not be +#' answered by recommending the other optimization method, because the other +#' method fails on the same fit. +#' +#' It refuses on ANY NA rather than on all of them, unlike the restart guard, +#' because the two situations differ. A failed restart is one starting point out +#' of eleven and the survivors are a legitimate answer, so only all-failed is +#' refused. A grid outcome is a function of the model, not of a starting point: +#' one NA there means the model cannot produce an outcome and every other grid +#' point's value is NA too. Refusing on any is also exactly the condition that +#' already failed, since max() propagates a single NA, so this narrows nothing. +#' +#' @param all_outcomes A numeric vector, the estimated outcome at each grid +#' intervention. +#' @param aliased_coef_names A character vector, the names of the outcome +#' model's coefficients glm() could not estimate. Empty means full rank; see +#' refuse_if_all_restarts_failed(). +#' +#' @return Invisibly NULL when every grid outcome is a number. Raises +#' otherwise. +#' +#' @noRd +refuse_if_no_grid_outcome <- function(all_outcomes, + aliased_coef_names = character(0)) { + if (anyNA(all_outcomes)) { + if (length(aliased_coef_names) > 0) { + stop(rank_deficient_outcome_message(aliased_coef_names)) + } + stop(paste( + "The estimated outcome could not be computed at every intervention on", + "the grid, so none of them can be compared against the outcome goal.", + "This means the outcome model produced a missing estimated outcome.", + "Please check the fitted outcome model and the values supplied for the", + "center weights, the center characteristics and the time effect." + )) + } + invisible(NULL) +} + + #' select_restart_within_bounds #' #' @description Internal function that picks the recommended intervention out @@ -247,6 +367,14 @@ select_restart_within_bounds <- function(restart_points, #' @param power_goal_cluster_id A character string. The name of a column in the #' data identifying the stage-1 centers, used to compute the stage-1 design #' effect when icc is non-zero. Default NULL. +#' @param aliased_coef_names A character vector. The names of the outcome +#' model's coefficients glm() could not estimate, which the caller reads off +#' the fitted model. Empty (the default) means the model is of full rank. This +#' function is given coefficient VECTORS rather than the model, so it cannot +#' establish this for itself: an NA it sees has already been summed into the +#' center-level effects and no longer names the term it came from. Used only to +#' say which of the two causes stopped an optimization that could compute no +#' outcome at all, since the two need opposite advice. #' #' @return List( #' est_rec_int = recommended interventions, @@ -289,7 +417,8 @@ get_recommended_interventions <- function( patients_per_center_in_next_stage, outcome_name, icc = NULL, - power_goal_cluster_id = NULL) { + power_goal_cluster_id = NULL, + aliased_coef_names = character(0)) { # check if power goal is null, if not, calculate the desired outcome # value needed to achieve the power goal if (!is.null(power_goal)) { @@ -474,6 +603,19 @@ get_recommended_interventions <- function( all_outcomes <- sapply(all_results, function(x) x$outcome) all_costs <- sapply(all_results, function(x) x$cost) + # an NA among the grid outcomes is refused here, before it reaches the + # goal comparison below. max() of anything holding an NA or NaN is NA, so + # "NA >= goal" is not FALSE but the R error "missing value where + # TRUE/FALSE needed", which is what this used to fail with -- no mention + # of the outcome, the model, or what to do. This is a GUARD and not a + # change of behaviour: anyNA(all_outcomes) holds exactly when + # is.na(max(all_outcomes)) does, i.e. exactly on the inputs that already + # failed, so no grid search that used to return one still raises. + # The numerical optimizer reaches the same situation as an all-failed + # restart set, since solnl() cannot optimize an NA objective either, and + # refuse_if_all_restarts_failed() names it there. + refuse_if_no_grid_outcome(all_outcomes, aliased_coef_names) + # find the maximum outcome max_outcome <- max(all_outcomes) @@ -635,7 +777,10 @@ get_recommended_interventions <- function( # by which.max() itself, which skips NAs: the surviving restarts are # compared and the winner's column of results_int_components is the # column that same restart wrote, since both are indexed by restart. - refuse_if_all_restarts_failed(results) + # The aliased coefficient names are what lets the refusal say WHY, since + # a rank-deficient model and a hard optimization both land here and want + # opposite advice. + refuse_if_all_restarts_failed(results, aliased_coef_names) max_position <- which.max(results) max_achievable_outcome <- results[max_position] @@ -708,9 +853,10 @@ get_recommended_interventions <- function( } # if numerical solution fails to find a solution. The same refusal the - # max-outcome loop above makes, from the same place, so the two loops - # cannot come to describe the same situation differently. - refuse_if_all_restarts_failed(cost_results) + # max-outcome loop above makes, from the same place and with the same + # cause, so the two loops cannot come to describe the same situation + # differently. + refuse_if_all_restarts_failed(cost_results, aliased_coef_names) # Choosing among the restarts, and making the winner implementable, is # one decision and lives in select_restart_within_bounds(): the in-box diff --git a/R/rec_int_processor.R b/R/rec_int_processor.R index 7d94181..a431f72 100644 --- a/R/rec_int_processor.R +++ b/R/rec_int_processor.R @@ -37,6 +37,19 @@ rec_int_processor <- function( coef_mapping <- term_coef_names(model) all_coefs <- coef(model) + # the coefficients glm() could not estimate. A term aliased with another gets + # an NA coefficient, every outcome computed from it is NA, and no optimizer + # can then find anything: the numerical one fails at every restart and the + # grid search cannot compare any grid point against the goal. Read off the + # model HERE, where the model is, because that is the only place the aliased + # TERMS can still be named -- an NA is summed into the center-level effects + # one step below and comes out carrying their names instead. Handed to the + # optimizers so their refusal can distinguish a rank-deficient fit, which no + # optimization method can help with, from a genuine optimization failure, + # which switching method may. Empty for a full-rank fit, which is every fit + # the package's own path produces unless the data makes it otherwise. + aliased_coef_names <- names(all_coefs)[is.na(all_coefs)] + # get coefficients for the intervention components intervention_components_coeff <- model$coefficients[c("(Intercept)", intervention_components)] @@ -199,7 +212,8 @@ rec_int_processor <- function( patients_per_center_in_next_stage = patients_per_center_in_next_stage, outcome_name = outcome_name, icc = icc, - power_goal_cluster_id = power_goal_cluster_id + power_goal_cluster_id = power_goal_cluster_id, + aliased_coef_names = aliased_coef_names ) list( diff --git a/R/validate_inputs.R b/R/validate_inputs.R index 7908eec..c4e02ae 100644 --- a/R/validate_inputs.R +++ b/R/validate_inputs.R @@ -353,55 +353,12 @@ validate_inputs <- function( # is not automatically a set of weights either -- with every # center_sample_size zero it is 0/0, i.e. all NaN. if (!is.null(center_weights_for_outcome_goal)) { - # every weight has to be a number before either of the checks that follow - # can be made at all: any(NA < 0) is not FALSE, it is the R error "missing - # value where TRUE/FALSE needed", and that is what a missing weight used to - # surface as at the sum check below. An Inf was refused there, but for - # summing to Inf rather than for being one, and an Inf beside a -Inf sums to - # NaN and hit the same opaque comparison. None of these is a weight, so each - # is named here instead. - if (!all(is.finite(center_weights_for_outcome_goal))) { - stop(paste( - "values in center_weights_for_outcome_goal must all be finite", - "numbers. NA, NaN, Inf and -Inf are not weights. If the weights were", - "not supplied, they were derived from the 'center_sample_size'", - "column, which cannot be zero for every center." - )) - } - # the weights are a convex combination over the centers: the reported - # outcome is sum(weight_i * outcome_i) over the center-level effects, so it - # is a weighted MEAN of the per-center outcomes and has to lie between the - # smallest and the largest of them. A negative weight breaks exactly that, - # and summing to 1 does not rule one out: weights of -10 and 11 sum to 1 and - # gave 10.95 for a logit outcome, i.e. a reported "probability" outside - # [0, 1] for a binary outcome. That is a wrong number rather than an error, - # so it is refused here and not clamped: a caller who passed a negative - # weight did not mean a weight of 0, and quietly substituting one would - # answer a question they did not ask. - # - # A weight of exactly 0 is allowed, deliberately. It means a center the - # recommendation is not being computed for, which is a meaningful thing to - # ask for and is what the package itself builds from - # center_effects_optimization_values: the named center gets weight 1 and - # every other center 0. Refusing 0 would refuse that documented path. All - # weights being 0 is a different matter, and it is already refused by the - # sum check below, which is what keeps the renormalisation from dividing by - # zero: a vector summing to 0 is 1 away from 1, not within 0.001 of it. - # Compared against a small negative tolerance rather than against 0. A - # weight a caller computes as a residual, one minus the others, can land a - # few floating-point units below zero while the vector still sums to 1, and - # refusing that would contradict allowing a weight of exactly 0. Anything - # further below zero than rounding explains is a weight the caller meant to - # be negative. - if (any(center_weights_for_outcome_goal < -8 * .Machine$double.eps)) { - stop(paste( - "values in center_weights_for_outcome_goal must be non-negative.", - "The weights average the per-center outcomes, so a negative weight", - "makes the reported outcome fall outside the range of the outcomes it", - "averages, and for a binary outcome outside [0, 1]. A weight of 0 is", - "allowed and excludes that center from the average." - )) - } + # what a set of weights has to be at all, i.e. finite and non-negative, + # before the sum check below can be made on it. Shared with the exported + # get_confidence_set(), so the two entry points cannot come to refuse + # different things: see refuse_invalid_center_weights() for why each check + # is there. + refuse_invalid_center_weights(center_weights_for_outcome_goal) weights_sum <- sum(center_weights_for_outcome_goal) if (abs(weights_sum - 1) >= 0.001) { stop(paste( @@ -1169,3 +1126,89 @@ validate_inputs <- function( prev_recommended_interventions = prev_recommended_interventions ) } + + +#' refuse_invalid_center_weights +#' +#' @description Internal guard for the center weights: refuses a vector that +#' is not a set of weights at all, i.e. one holding a non-finite or a negative +#' value. +#' +#' @details Both checks are here rather than at one caller because +#' validate_inputs() is not the only entry point the weights arrive through. +#' The exported get_confidence_set() takes them directly and does not go +#' through validate_inputs(), so a negative weight passed to it reached the +#' interval unchecked and reported a "probability" above 1 for a binary +#' outcome. That is the same wrong number validate_inputs() already refused, +#' from the other door, and the two must refuse it in the same words: a caller +#' who moves between the two entry points should not be told two different +#' things about the same vector. +#' +#' Finiteness first, because neither of the comparisons that follow can be made +#' at all otherwise: any(NA < 0) is not FALSE, it is the R error "missing value +#' where TRUE/FALSE needed", and that is what a missing weight used to surface +#' as at the sum check in validate_inputs(). An Inf was refused there, but for +#' summing to Inf rather than for being one, and an Inf beside a -Inf sums to +#' NaN and hit the same opaque comparison. None of these is a weight, so each +#' is named here instead. +#' +#' Then non-negativity. The weights are a convex combination over the centers: +#' the reported outcome is sum(weight_i * outcome_i) over the center-level +#' effects, so it is a weighted MEAN of the per-center outcomes and has to lie +#' between the smallest and the largest of them. A negative weight breaks +#' exactly that, and summing to 1 does not rule one out: weights of -10 and 11 +#' sum to 1 and gave 10.95 for a logit outcome, i.e. a reported "probability" +#' outside [0, 1] for a binary outcome. That is a wrong number rather than an +#' error, so it is refused and not clamped: a caller who passed a negative +#' weight did not mean a weight of 0, and quietly substituting one would answer +#' a question they did not ask. +#' +#' A weight of exactly 0 is allowed, deliberately. It means a center the +#' recommendation is not being computed for, which is a meaningful thing to ask +#' for and is what the package itself builds from +#' center_effects_optimization_values: the named center gets weight 1 and every +#' other center 0. Refusing 0 would refuse that documented path. All weights +#' being 0 is a different matter, and it is refused by validate_inputs()' sum +#' check, which is what keeps its renormalisation from dividing by zero: a +#' vector summing to 0 is 1 away from 1, not within 0.001 of it. +#' +#' Compared against a small negative tolerance rather than against 0. A weight +#' a caller computes as a residual, one minus the others, can land a few +#' floating-point units below zero while the vector still sums to 1, and +#' refusing that would contradict allowing a weight of exactly 0. Anything +#' further below zero than rounding explains is a weight the caller meant to be +#' negative. +#' +#' The SUM check is deliberately NOT here. validate_inputs() both refuses a sum +#' far from 1 and renormalises what it accepts, and the renormalisation is what +#' makes the tolerance safe; get_confidence_set() has no such step and is +#' documented as taking the weights the optimization ran with, so requiring a +#' unit sum of it would refuse the vector its own caller has already +#' normalised for it. +#' +#' @param center_weights_for_outcome_goal A numeric vector of center weights. +#' +#' @return Invisibly NULL when every weight is finite and non-negative. Raises +#' otherwise. +#' +#' @noRd +refuse_invalid_center_weights <- function(center_weights_for_outcome_goal) { + if (!all(is.finite(center_weights_for_outcome_goal))) { + stop(paste( + "values in center_weights_for_outcome_goal must all be finite", + "numbers. NA, NaN, Inf and -Inf are not weights. If the weights were", + "not supplied, they were derived from the 'center_sample_size'", + "column, which cannot be zero for every center." + )) + } + if (any(center_weights_for_outcome_goal < -8 * .Machine$double.eps)) { + stop(paste( + "values in center_weights_for_outcome_goal must be non-negative.", + "The weights average the per-center outcomes, so a negative weight", + "makes the reported outcome fall outside the range of the outcomes it", + "averages, and for a binary outcome outside [0, 1]. A weight of 0 is", + "allowed and excludes that center from the average." + )) + } + invisible(NULL) +} diff --git a/man/get_confidence_set.Rd b/man/get_confidence_set.Rd index 49587cb..0ee4288 100644 --- a/man/get_confidence_set.Rd +++ b/man/get_confidence_set.Rd @@ -40,7 +40,9 @@ should be included in the outcome model.} \item{center_weights_for_outcome_goal}{A numeric vector. Specifies the 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.} +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.} \item{include_time_effects}{A boolean. Specifies whether the fixed time effects should be included in the outcome model.} diff --git a/man/get_recommended_interventions.Rd b/man/get_recommended_interventions.Rd index 55883cf..2a21aea 100644 --- a/man/get_recommended_interventions.Rd +++ b/man/get_recommended_interventions.Rd @@ -30,7 +30,8 @@ get_recommended_interventions( patients_per_center_in_next_stage, outcome_name, icc = NULL, - power_goal_cluster_id = NULL + power_goal_cluster_id = NULL, + aliased_coef_names = character(0) ) } \arguments{ @@ -145,6 +146,15 @@ get_power_desired_outcome.} \item{power_goal_cluster_id}{A character string. The name of a column in the data identifying the stage-1 centers, used to compute the stage-1 design effect when icc is non-zero. Default NULL.} + +\item{aliased_coef_names}{A character vector. The names of the outcome +model's coefficients glm() could not estimate, which the caller reads off +the fitted model. Empty (the default) means the model is of full rank. This +function is given coefficient VECTORS rather than the model, so it cannot +establish this for itself: an NA it sees has already been summed into the +center-level effects and no longer names the term it came from. Used only to +say which of the two causes stopped an optimization that could compute no +outcome at all, since the two need opposite advice.} } \value{ List( diff --git a/tests/testthat/test-confidence-set-contract.R b/tests/testthat/test-confidence-set-contract.R index 6711ed9..f6dd4e5 100644 --- a/tests/testthat/test-confidence-set-contract.R +++ b/tests/testthat/test-confidence-set-contract.R @@ -717,3 +717,285 @@ test_that("the estimated outcome lies inside its interval at every requested per intervals <- unique(lapply(results, function(r) unname(r$est_outcome_ci))) expect_equal(length(intervals), length(periods)) }) + + +test_that("the exported get_confidence_set() refuses a bad center weight", { + # get_confidence_set() is EXPORTED and does not go through validate_inputs(), + # so the weight guards lago_optimization() has do not cover a direct caller. + # The weights average the per-center outcomes, so a negative one puts the + # interval outside the range of the intervals it averages: c(8.5, -8, 0.5) + # sums to exactly 1, passed every check this function made, and reported a + # CI_upper_bound of 1.014 for a BINARY outcome -- a "probability" above 1, + # silently. This function already validates its own link and its own + # coefficient-to-predictor match for exactly this reason (it cannot trust its + # caller), and the weights are the same kind of argument. + d <- as.data.frame(BB_data) + d$center <- factor(rep_len(paste0("c", 1:3), nrow(d))) + components <- c("coaching_updt", "launch_duration") + model <- glm( + pp3_oxytocin_mother ~ center + coaching_updt + launch_duration, + data = d, family = binomial(link = "logit") + ) + + call_cs <- function(w, include_center_effects = TRUE) { + suppressWarnings(get_confidence_set( + predictors_data = d[, c("center", components), drop = FALSE], + include_center_effects = include_center_effects, + center_weights_for_outcome_goal = w, + intervention_components = components, + outcome_data = d$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(8, 1), + cost_list_of_vectors = list(c(0, 1.7), c(0, 8)), + rec_int = c(20, 3) + )) + } + + # a negative weight is refused, in the words validate_inputs() uses, so a + # caller who moves between the two entry points is told the same thing + expect_error(call_cs(c(8.5, -8, 0.5)), "must be non-negative") + expect_error(call_cs(c(0.75, -0.5, 0.75)), "must be non-negative") + expect_error(call_cs(c(20.5, -20, 0.5)), "must be non-negative") + # summing to 1 is what let it through, so the fixture asserts that it does + expect_identical(sum(c(8.5, -8, 0.5)), 1) + expect_identical(sum(c(0.75, -0.5, 0.75)), 1) + + # the boundary, so the tolerance cannot be widened without a test failing. + # It admits floating-point noise from a residual weight and nothing larger. + expect_error(call_cs(c(0.5, -1e-9, 0.5)), "must be non-negative") + expect_error(call_cs(c(0.5, -1e-12, 0.5)), "must be non-negative") + + # non-finite weights are named rather than reaching the interval, where they + # made every bound NA and rec_int_ci NULL with nothing said about why + for (bad in list( + c(0.5, NA_real_, 0.5), c(0.5, NaN, 0.5), + c(0.5, Inf, 0.5), c(0.5, 0.5, -Inf) + )) { + expect_error(call_cs(bad), "must all be finite") + } + + # a weight of exactly 0 is still ALLOWED, and so is a residual weight a hair + # below zero: the guard must not narrow what the function accepts. The + # single-named-center path the package itself builds is exactly a vector of + # one 1 and the rest 0. + residual <- -.Machine$double.eps / 2 + expect_lt(residual, 0) + for (good in list( + c(1, 1, 1) / 3, c(0.5, 0, 0.5), c(1, 0, 0), c(0, 0, 1), + c(0.5, residual, 0.5) + )) { + res <- call_cs(good) + expect_false(is.null(res$rec_int_ci)) + expect_true(all(res$rec_int_ci >= 0 & res$rec_int_ci <= 1)) + } + + # and the numbers a compliant vector produces are unchanged, so the guard + # only removed the refused cases + uniform <- call_cs(c(1, 1, 1) / 3) + expect_true(all(uniform$cs$CI_upper_bound <= 1)) + expect_true(all(uniform$cs$CI_lower_bound >= 0)) + + # the weights are only USED when the fixed center effects are included, so + # they are only checked then: a caller who is not asking for center effects + # passes the default 1 and must not be refused for a vector nobody reads. + expect_error( + suppressWarnings(get_confidence_set( + predictors_data = d[, components, drop = FALSE], + center_weights_for_outcome_goal = c(-1, 2), + intervention_components = components, + outcome_data = d$pp3_oxytocin_mother, + fitted_model = glm( + pp3_oxytocin_mother ~ coaching_updt + launch_duration, + data = d, family = binomial(link = "logit") + ), + 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(8, 1), + cost_list_of_vectors = list(c(0, 1.7), c(0, 8)), + rec_int = c(20, 3) + )), + NA + ) + + # the same weights through the guarded primary path are refused there too, so + # the two entry points agree rather than one being stricter + primary <- function() { + suppressWarnings(suppressMessages(lago_optimization( + data = d, + outcome_name = "pp3_oxytocin_mother", + outcome_type = "binary", + glm_family = "binomial", + intervention_components = components, + intervention_lower_bounds = c(1, 1), + intervention_upper_bounds = c(40, 5), + cost_list_of_vectors = list(c(0, 1.7), c(0, 8)), + outcome_goal = 0.85, + include_center_effects = TRUE, + center_weights_for_outcome_goal = c(8.5, -8, 0.5), + include_confidence_set = FALSE, + quiet = TRUE + ))) + } + expect_error(primary(), "must be non-negative") + expect_identical( + tryCatch(primary(), error = conditionMessage), + tryCatch(call_cs(c(8.5, -8, 0.5)), error = conditionMessage) + ) +}) + + +test_that("a binary outcome's interval is built on the link it was fitted on", { + # The binary branch applied expit() and the logit delta-method factor + # p * (1 - p) UNCONDITIONALLY, keyed on outcome_type and ignoring link. A + # binomial model fitted with link = "identity" -- which lago_optimization() + # accepts -- therefore got its interval on the logit scale: 0.636 to 0.655 + # where the identity-scale interval is 0.558 to 0.642, an interval that does + # not even contain the point estimate the same run reported. The scale an + # interval belongs on is a property of the LINK, not of the outcome's type, + # which is what get_outcome() has keyed on all along -- hence the point + # estimate being right while the interval was not. + # + # The fixture converges cleanly (3 IRLS iterations, no glm warning), because a + # binomial identity fit that does not converge would leave it open whether the + # interval or the fit was the problem. It keeps the fitted probabilities well + # inside (0, 1): a linear probability model on a grid, replicated ten rows per + # cell so the design is balanced. + cells <- expand.grid(x1 = seq(0, 9, length.out = 40), x2 = 1:5) + p <- 0.30 + 0.02 * cells$x1 + 0.04 * cells$x2 + expect_true(all(p > 0.05 & p < 0.95)) + d <- do.call(rbind, lapply(seq_len(nrow(cells)), function(i) { + ones <- round(p[i] * 10) + data.frame( + x1 = cells$x1[i], x2 = cells$x2[i], + y = c(rep(1L, ones), rep(0L, 10 - ones)) + ) + })) + + # the fixture's own precondition: the fit converges and glm() says nothing + fit_warnings <- character(0) + model <- withCallingHandlers( + glm(y ~ x1 + x2, data = d, family = binomial(link = "identity")), + warning = function(w) { + fit_warnings <<- c(fit_warnings, conditionMessage(w)) + invokeRestart("muffleWarning") + } + ) + expect_true(model$converged) + expect_length(fit_warnings, 0) + expect_false(anyNA(coef(model))) + + res <- suppressWarnings(suppressMessages(lago_optimization( + data = d, + outcome_name = "y", + outcome_type = "binary", + glm_family = "binomial", + link = "identity", + intervention_components = c("x1", "x2"), + intervention_lower_bounds = c(0, 1), + intervention_upper_bounds = c(9, 5), + cost_list_of_vectors = list(c(0, 1), c(0, 2)), + outcome_goal = 0.6, + confidence_set_grid_step_size = c(1, 1), + quiet = TRUE + ))) + + # THE assertion. The interval is recomputed by hand from glm() itself -- the + # linear predictor and its standard error, with no link transformation, + # because on the identity link the outcome IS the linear predictor and the + # delta-method factor is the derivative of the identity map, i.e. 1. This + # oracle does not go through the package. + critical_value <- qnorm(0.975) + identity_interval <- function(x) { + row <- c(1, x) + point <- as.numeric(row %*% coef(model)) + std_error <- sqrt(as.numeric(t(row) %*% vcov(model) %*% row)) + c( + round(point - critical_value * std_error, 3), + round(point + critical_value * std_error, 3) + ) + } + expect_identical( + unname(res$est_outcome_ci), + identity_interval(res$rec_int) + ) + + # and what the defect produced instead, so the test says which scale is wrong + # rather than only that the number changed: expit() of the point with the + # logit delta factor accounts for the OLD interval exactly, and it is a + # different interval from the one now reported. + row <- c(1, res$rec_int) + point <- as.numeric(row %*% coef(model)) + std_error <- sqrt(as.numeric(t(row) %*% vcov(model) %*% row)) + logit_point <- rje::expit(point) + logit_std_error <- std_error * logit_point * (1 - logit_point) + logit_interval <- c( + round(logit_point - critical_value * logit_std_error, 3), + round(logit_point + critical_value * logit_std_error, 3) + ) + expect_false(identical(logit_interval, identity_interval(res$rec_int))) + expect_false(identical(unname(res$est_outcome_ci), logit_interval)) + + # the interval contains the point estimate, which the logit-scale one did not: + # the estimated outcome was 0.6 and the reported interval 0.636 to 0.655 + expect_true( + res$est_outcome_ci[["lower"]] <= res$est_outcome_goal && + res$est_outcome_goal <= res$est_outcome_ci[["upper"]] + ) + + # every confidence-set row, not only rec_int, and each against its own + # coordinates. A set whose intervals are on the wrong scale can still be + # non-empty, so containment of the goal alone would not catch this. + expect_gt(nrow(res$cs), 0) + for (i in seq_len(nrow(res$cs))) { + expect_identical( + c(res$cs$CI_lower_bound[i], res$cs$CI_upper_bound[i]), + identity_interval(c(res$cs$x1[i], res$cs$x2[i])) + ) + } + # and they cover the goal, as membership claims, on the corrected scale + expect_true(all( + res$cs$CI_lower_bound <= 0.6 & res$cs$CI_upper_bound >= 0.6 + )) + + # the LOGIT case is untouched, which is what confines the change to the link + # that was wrong. Same oracle, the delta-method interval glm()'s own binary + # machinery reports, on the same package's default binary configuration. + logit_res <- suppressWarnings(suppressMessages(lago_optimization( + data = as.data.frame(BB_data), + outcome_name = "pp3_oxytocin_mother", + outcome_type = "binary", + glm_family = "binomial", + intervention_components = c("coaching_updt", "launch_duration"), + intervention_lower_bounds = c(1, 1), + intervention_upper_bounds = c(40, 5), + cost_list_of_vectors = list(c(0, 1.7), c(0, 8)), + outcome_goal = 0.85, + confidence_set_grid_step_size = c(8, 1), + quiet = TRUE + ))) + logit_model <- glm( + pp3_oxytocin_mother ~ coaching_updt + launch_duration, + data = as.data.frame(BB_data), family = binomial(link = "logit") + ) + logit_row <- c(1, logit_res$rec_int) + logit_eta <- as.numeric(logit_row %*% coef(logit_model)) + logit_p <- rje::expit(logit_eta) + logit_se <- sqrt( + as.numeric(t(logit_row) %*% vcov(logit_model) %*% logit_row) + ) * logit_p * (1 - logit_p) + expect_identical( + unname(logit_res$est_outcome_ci), + c( + round(logit_p - critical_value * logit_se, 3), + round(logit_p + critical_value * logit_se, 3) + ) + ) +}) diff --git a/tests/testthat/test-optimization.R b/tests/testthat/test-optimization.R index f2fe51b..382d7f6 100644 --- a/tests/testthat/test-optimization.R +++ b/tests/testthat/test-optimization.R @@ -234,25 +234,29 @@ test_that("weights length is validated against the number of observations, not c }) -test_that("an unsolvable numerical optimization says to use grid_search", { +test_that("a rank-deficient fit is named as such, not sent to grid_search", { # END TO END, which is what the unit test of the guard cannot show: that a # configuration a user can actually pass reaches it. Every solnl() restart of # the max-achievable-outcome loop fails on this fit, and the loop records each # failure as NA. With all of them NA, which.max() is integer(0), the outcome # read out of it is numeric(0), and the goal comparison that follows used to # fail with the base-R error "argument is of length zero" -- no mention of the - # optimizer, the goal, or what to do instead. The cost loop 60 lines further - # down already refused the same situation with a message naming grid_search; - # this path simply never reached it. + # optimizer, the goal, or what to do instead. # # What makes every restart fail here: 6 recycled centers against 3 recycled # periods are ALIASED (6 and 3 share a factor, so the period indicators are # linear combinations of the center ones), glm() returns NA for the aliased # period coefficients, and those NAs flow into the center-level effects. Every - # outcome the solver is asked for is then NA, so every restart fails. That is - # a rank-deficient model rather than a hard optimization, and it is the - # smallest reproducer of the all-failed path; see the note at the end of this - # test about the separate defect it also exposes. + # outcome the solver is asked for is then NA, so every restart fails. + # + # That is a rank-deficient MODEL rather than a hard optimization, and it is + # what the message has to say, because the two want opposite advice. The + # message used to recommend 'grid_search' for every cause -- and grid_search + # fails on this same fit, since an NA outcome is an NA outcome whatever + # searches over it. Both halves are asserted below: the numerical path names + # the rank deficiency, and the grid path, which the old message sent the user + # to, now names the same cause instead of failing with the base-R "missing + # value where TRUE/FALSE needed". bbp <- as.data.frame(BB_proportions) bbp$center <- factor(rep_len(paste0("s", 1:6), nrow(bbp))) bbp$period <- factor(rep_len(1:3, nrow(bbp))) @@ -284,27 +288,30 @@ test_that("an unsolvable numerical optimization says to use grid_search", { ))) } - # the message a user can act on, and the recommendation it makes - expect_error(run(), "Numerical optimization failed to find a solution") - expect_error(run(), "'grid_search'") - # and NOT the base-R error the missing guard produced. This is the assertion - # that fails without the fix: the old message was exactly this string. + # the cause a user can act on, and the terms that carry it. Naming the + # coefficients is the actionable part: they are what has to be dropped or + # combined, and nothing else in the run reports them. + expect_error(run(), "rank-deficient") + expect_error(run(), "period2, period3") + expect_error(run(), "[Dd]rop or combine") + # and NOT the base-R error the missing guard produced expect_error(run(), "^(?!.*argument is of length zero).*$", perl = TRUE) + # and NOT the advice that cannot help. The old message recommended + # 'grid_search' here, which fails on this very fit, so a message naming it is + # the defect: it sends the user to a method that also cannot work. This is + # the assertion that fails on the unfixed source, where the message is + # exactly that recommendation. + expect_error(run(), "^(?!.*grid_search).*$", perl = TRUE) + expect_error(run(), "^(?!.*more than\\s+three intervention components).*$", + perl = TRUE) - # NOT asserted here: that following the message's advice succeeds on THIS - # data. It does not, and deliberately so -- the cause on this fixture is the - # rank-deficient fit above, which grid_search cannot help with either. Asked - # for grid_search on the same inputs, the NA outcome reaches that optimizer's - # own goal comparison and it fails with the base-R "missing value where - # TRUE/FALSE needed", i.e. the same class of opaque failure this test fixes on - # the numerical path, in the other optimizer and from a different cause. - # - # That is a SEPARATE defect from the one under test and is left alone here: it - # is an unguarded NA outcome from an aliased model, not an unguarded - # all-failed restart set, and fixing it means refusing a rank-deficient fit, - # which is a wider behaviour change than this one. Pinned as the current - # behaviour so that the day it is fixed, this assertion is what says so. - grid_error <- tryCatch( + # the other half: the method the old message sent the user to. It fails on + # this input too -- the NA outcome reaches its own goal comparison -- and it + # used to fail with the base-R "missing value where TRUE/FALSE needed", the + # same class of opaque failure, in the other optimizer. It now names the same + # cause as the numerical path, because the cause IS the same and neither + # method can proceed on it. + grid_run <- function() { suppressWarnings(suppressMessages(lago_optimization( data = bbp, outcome_name = "EBP_proportions", @@ -323,10 +330,18 @@ test_that("an unsolvable numerical optimization says to use grid_search", { optimization_grid_search_step_size = c(4, 1), include_confidence_set = FALSE, quiet = TRUE - ))), - error = function(e) conditionMessage(e) + ))) + } + expect_error(grid_run(), "rank-deficient") + expect_error(grid_run(), "period2, period3") + expect_error(grid_run(), + "^(?!.*missing value where TRUE/FALSE needed).*$", perl = TRUE) + # both optimizers say the SAME thing about the same fit, which is the point of + # sharing the message: the cause is the model, not the method. + expect_identical( + tryCatch(run(), error = conditionMessage), + tryCatch(grid_run(), error = conditionMessage) ) - expect_match(grid_error, "missing value where TRUE/FALSE needed") # and on a FULL-RANK fit of the same shape, where the aliasing is gone because # 5 centers and 3 periods share no factor, both optimizers succeed. This is diff --git a/tests/testthat/test-outcome-internals.R b/tests/testthat/test-outcome-internals.R index 5c96ff5..7115ccc 100644 --- a/tests/testthat/test-outcome-internals.R +++ b/tests/testthat/test-outcome-internals.R @@ -1417,7 +1417,10 @@ test_that("an all-failed restart set is refused by one guard, for both loops", { "refuse_if_all_restarts_failed", "LAGO" ) - # all failed: refused, and by the message that names the way out + # all failed on a FULL-RANK model: refused, and by the message that names the + # way out. That the model is of full rank is what makes recommending the other + # optimizer sensible, and it is stated by passing no aliased coefficient + # names -- see the test below for the other branch. expect_error( refuse_if_all_restarts_failed(rep(NA_real_, 11)), "Numerical optimization failed to find a solution" @@ -1426,6 +1429,12 @@ test_that("an all-failed restart set is refused by one guard, for both loops", { refuse_if_all_restarts_failed(rep(NA_real_, 11)), "'grid_search'" ) + # passing an EMPTY set of aliased names is the same as passing none, so the + # argument's default cannot drift away from what the callers rely on + expect_error( + refuse_if_all_restarts_failed(rep(NA_real_, 11), character(0)), + "Numerical optimization failed to find a solution" + ) # a single restart, failed, is still every restart expect_error( refuse_if_all_restarts_failed(NA_real_), @@ -1459,6 +1468,202 @@ test_that("an all-failed restart set is refused by one guard, for both loops", { }) +test_that("a rank-deficient fit is not answered by recommending grid_search", { + # The refusal has one CONDITION and two CAUSES, and the causes want opposite + # advice. Recommending the other optimization method is only sensible when the + # model could be estimated: an aliased term's coefficient is NA, every outcome + # computed from it is NA, and no search over interventions can recover from + # that -- the grid search fails on the same fit, with its own unguarded + # comparison. So the message that names grid_search must NOT be what a + # rank-deficient fit gets, which is what a single message for both did. + refuse_if_all_restarts_failed <- getFromNamespace( + "refuse_if_all_restarts_failed", "LAGO" + ) + refuse_if_no_grid_outcome <- getFromNamespace( + "refuse_if_no_grid_outcome", "LAGO" + ) + + aliased <- c("period2", "period3") + + # the rank-deficient branch: the cause, the terms, and what to do about them + expect_error( + refuse_if_all_restarts_failed(rep(NA_real_, 11), aliased), + "rank-deficient" + ) + expect_error( + refuse_if_all_restarts_failed(rep(NA_real_, 11), aliased), + "period2, period3" + ) + expect_error( + refuse_if_all_restarts_failed(rep(NA_real_, 11), aliased), + "[Dd]rop or combine" + ) + # and NOT the advice that cannot help. This is the assertion the shared + # message fails: it recommended grid_search whatever the cause. + expect_error( + refuse_if_all_restarts_failed(rep(NA_real_, 11), aliased), + "^(?!.*grid_search).*$", + perl = TRUE + ) + # nor the cause it used to blame, which the aliased fit does not have: two + # intervention components, not more than three + expect_error( + refuse_if_all_restarts_failed(rep(NA_real_, 11), aliased), + "^(?!.*more than\\s+three intervention components).*$", + perl = TRUE + ) + + # the branch is on the ALIASED NAMES and not on the restarts, so a full-rank + # all-failed set still gets the wording a caller in that position needs. Both + # branches from the same call, so neither can be reached by accident. + full_rank <- tryCatch( + refuse_if_all_restarts_failed(rep(NA_real_, 11)), + error = conditionMessage + ) + deficient <- tryCatch( + refuse_if_all_restarts_failed(rep(NA_real_, 11), aliased), + error = conditionMessage + ) + expect_false(identical(full_rank, deficient)) + expect_match(full_rank, "grid_search") + expect_no_match(deficient, "grid_search") + + # a SOME-failed set still passes through, whatever the model: the branch must + # not turn the guard into one that fires on any NA. A restart is a starting + # point, and the survivors are a legitimate answer. + expect_silent(refuse_if_all_restarts_failed(c(NA_real_, 0.5), aliased)) + expect_null(refuse_if_all_restarts_failed(c(NA_real_, 0.5), aliased)) + + # the grid search's own guard, which is the second half of the same defect: + # an NA grid outcome reached max() >= goal and failed with the base-R + # "missing value where TRUE/FALSE needed". It says the same thing about the + # same cause, since the cause is the model. + expect_error(refuse_if_no_grid_outcome(c(0.4, NA, 0.6), aliased), + "rank-deficient") + expect_error(refuse_if_no_grid_outcome(c(0.4, NA, 0.6), aliased), + "period2, period3") + expect_error(refuse_if_no_grid_outcome(c(NA_real_, NA_real_), aliased), + "rank-deficient") + expect_identical( + tryCatch(refuse_if_no_grid_outcome(rep(NA_real_, 3), aliased), + error = conditionMessage), + deficient + ) + # and on a full-rank model it names the outcome rather than the model, and + # still does not recommend a method + full_rank_grid <- tryCatch( + refuse_if_no_grid_outcome(c(0.4, NA, 0.6)), + error = conditionMessage + ) + expect_match(full_rank_grid, "could not be computed at every intervention") + expect_no_match(full_rank_grid, "rank-deficient") + + # it fires on ANY NA, unlike the restart guard, and that is exactly the + # condition that already failed: max() of a vector holding one NA is NA, so + # "NA >= goal" was already the R error. Refusing here narrows nothing. + expect_true(is.na(max(c(0.4, NA, 0.6)))) + expect_error( + if (max(c(0.4, NA, 0.6)) >= 0.5) TRUE else FALSE, + "missing value where TRUE/FALSE needed" + ) + # a grid with every outcome a number passes through untouched + expect_silent(refuse_if_no_grid_outcome(c(0.4, 0.5, 0.6))) + expect_null(refuse_if_no_grid_outcome(c(0.4, 0.5, 0.6), aliased)) + expect_silent(refuse_if_no_grid_outcome(numeric(0), aliased)) +}) + + +test_that("the optimizers are told which coefficients could not be estimated", { + # The guards above can only distinguish the two causes if something reads the + # aliased names off the FITTED MODEL and hands them over. + # get_recommended_interventions() is given coefficient vectors, not the model, + # and by then an NA has been summed into the center-level effects and carries + # THEIR names -- so the terms could not be named from there even though the + # NA is visible. rec_int_processor() has the model, which is why the read + # happens there. + bbp <- as.data.frame(BB_proportions) + bbp$center <- factor(rep_len(paste0("s", 1:6), nrow(bbp))) + bbp$period <- factor(rep_len(1:3, nrow(bbp))) + model <- suppressWarnings(glm( + EBP_proportions ~ center + period + coaching_updt + launch_duration, + data = bbp, family = quasibinomial(link = "logit") + )) + all_coefs <- coef(model) + # the precondition, asserted rather than assumed + expect_true(anyNA(all_coefs)) + expect_identical(names(all_coefs)[is.na(all_coefs)], c("period2", "period3")) + + # what the NA looks like one step downstream, which is why the read cannot be + # deferred: it is still there, but it now spells the center-level effects + coef_mapping <- getFromNamespace("term_coef_names", "LAGO")(model) + named_predictors <- getFromNamespace("claimed_coef_names", "LAGO")( + model, coef_mapping, + c("(Intercept)", "coaching_updt", "launch_duration") + ) + fecn <- getFromNamespace("fixed_effect_coef_names", "LAGO") + center_coefs <- fecn("center", coef_mapping, names(all_coefs), + named_predictors) + period_coefs <- fecn("period", coef_mapping, names(all_coefs), + named_predictors) + intercept <- all_coefs["(Intercept)"] + center_level <- c(intercept, all_coefs[center_coefs] + intercept) + indicators <- getFromNamespace("time_effect_indicator", "LAGO")( + model, period_coefs, 1 + ) + center_level <- center_level + sum(indicators * all_coefs[period_coefs]) + expect_true(anyNA(center_level)) + # every entry is NA and each is named after a CENTER, so "period2, period3" + # is unrecoverable from here + expect_true(all(is.na(center_level))) + expect_false(any(c("period2", "period3") %in% names(center_level))) + + # and the intervention coefficients, which is the other vector the optimizer + # gets, are not NA at all: nothing there says the fit is rank-deficient + expect_false(anyNA(all_coefs[c("(Intercept)", "coaching_updt", + "launch_duration")])) + + # so the aliased names are passed in, and the refusal names them. Through the + # processor rather than the optimizer, which is the wiring under test. + err <- tryCatch( + suppressWarnings(suppressMessages(getFromNamespace( + "rec_int_processor", "LAGO" + )( + data = bbp, + model = model, + center_characteristics = NULL, + additional_covariates = NULL, + include_center_effects = TRUE, + include_time_effects = TRUE, + include_interaction_terms = FALSE, + main_components = NULL, + intervention_components = c("coaching_updt", "launch_duration"), + optimization_method = "numerical", + optimization_grid_search_step_size = NULL, + link = "logit", + center_weights_for_outcome_goal = rep(1 / 6, 6), + cost_list_of_vectors = list(c(0, 1.7), c(0, 8)), + intervention_lower_bounds = c(1, 1), + intervention_upper_bounds = c(40, 5), + outcome_goal = 0.85, + center_characteristics_optimization_values = NULL, + time_effect_optimization_value = 1, + lower_outcome_goal = FALSE, + prev_recommended_interventions = NULL, + shrinkage_threshold = 0.25, + power_goal = NULL, + power_goal_approach = "unconditional", + num_centers_in_next_stage = NULL, + patients_per_center_in_next_stage = NULL, + outcome_name = "EBP_proportions" + ))), + error = conditionMessage + ) + expect_match(err, "rank-deficient") + expect_match(err, "period2, period3") + expect_no_match(err, "grid_search") +}) + + test_that("with no restart in the box the winner is projected and recosted", { # THE POINT OF THE EXTRACTION. Both the projection and the cost recomputation # are only reachable when every restart left the box, which from the outside diff --git a/tests/testthat/test-validate-inputs.R b/tests/testthat/test-validate-inputs.R index 0ce9344..65351d3 100644 --- a/tests/testthat/test-validate-inputs.R +++ b/tests/testthat/test-validate-inputs.R @@ -359,6 +359,64 @@ test_that("center_weights_for_outcome_goal must be non-negative and finite", { ) }) +test_that("the weight guard is one function, shared with the confidence set", { + # The finiteness and non-negativity checks are a helper both entry points + # call, not a copy in each. That matters because the exported + # get_confidence_set() does not go through validate_inputs() and had NO weight + # checks at all: a negative weight reached its interval and reported a + # "probability" above 1. Two copies would let the two doors drift apart, which + # is the failure mode this pins -- the same vector must be refused in the same + # words wherever it is passed. + refuse_invalid_center_weights <- getFromNamespace( + "refuse_invalid_center_weights", "LAGO" + ) + + # what it refuses, and what it deliberately does not + expect_error(refuse_invalid_center_weights(c(-10, 11, 0)), + "must be non-negative") + expect_error(refuse_invalid_center_weights(c(0.5, -0.01, 0.51)), + "must be non-negative") + expect_error(refuse_invalid_center_weights(c(0.5, -1e-9, 0.5)), + "must be non-negative") + expect_error(refuse_invalid_center_weights(c(NA_real_, 0.5, 0.5)), + "must all be finite") + expect_error(refuse_invalid_center_weights(c(Inf, -Inf, 1)), + "must all be finite") + expect_silent(refuse_invalid_center_weights(c(0.5, 0.5, 0))) + expect_null(refuse_invalid_center_weights(c(1, 0, 0))) + expect_silent(refuse_invalid_center_weights(c(0.5, -.Machine$double.eps / 2, + 0.5))) + + # the SUM check is NOT in the helper, deliberately: validate_inputs() both + # refuses a sum far from 1 and renormalises what it accepts, and the + # renormalisation is what makes that tolerance safe. get_confidence_set() has + # no such step, so requiring a unit sum of it would refuse the vector its own + # caller already normalised. All-zero weights are therefore accepted by the + # helper and refused by validate_inputs(). + expect_silent(refuse_invalid_center_weights(c(0, 0, 0))) + expect_silent(refuse_invalid_center_weights(c(0.3, 0.3, 0.3))) + + # and validate_inputs() really routes through it rather than carrying its own + # copy: the messages are identical strings, not merely both matching a regex. + cw <- function(w) { + a <- vi_args(include_center_effects = TRUE, + center_weights_for_outcome_goal = w) + a$data <- data.frame( + mpg = c(21, 22, 23, 24, 25, 26), + gear = c(3, 4, 3, 4, 3, 4), + qsec = c(16, 17, 18, 19, 20, 21), + center = factor(c(1, 1, 2, 2, 3, 3)) + ) + suppressMessages(do.call(LAGO:::validate_inputs, a)) + } + for (bad in list(c(-10, 11, 0), c(NA_real_, 0.5, 0.5))) { + expect_identical( + tryCatch(cw(bad), error = conditionMessage), + tryCatch(refuse_invalid_center_weights(bad), error = conditionMessage) + ) + } +}) + test_that("the weight checks cover the DEFAULT weights, not just supplied", { # The checks are placed after every branch that can produce the weights, not # inside the one that reads them from the caller, so they cover the From 22f086f1a2a739034ff4d826c7270cb0d6938cc2 Mon Sep 17 00:00:00 2001 From: Ante Bing Date: Thu, 6 Aug 2026 23:07:20 +0000 Subject: [PATCH 2/4] Refuse a non-numeric weight vector, and name the saturated case Two points from review, both about a guard saying something a caller can act on. is.finite() is TRUE for every level of a factor, so a factor whose levels look like numbers passed the finiteness check, and the comparison after it then gave NA and raised the base-R "missing value where TRUE/FALSE needed" that the guard exists to replace. A non-numeric weight vector is now refused for what it is, with the same treatment for a character vector. This was inherited from the change that introduced the guard rather than new here. The rank-deficient refusal told the caller to drop or combine the collinear predictors. That is unfollowable on the route this branch newly names, where center-level data with one row per center is saturated rather than collinear: dropping a predictor leaves nothing to estimate the others from either. The message now covers that case and says what does help, which is more rows per center or fitting without the center effects. --- R/get_recommended_interventions.R | 7 ++++++- R/validate_inputs.R | 12 ++++++++++++ tests/testthat/test-validate-inputs.R | 10 ++++++++++ 3 files changed, 28 insertions(+), 1 deletion(-) diff --git a/R/get_recommended_interventions.R b/R/get_recommended_interventions.R index eaa5088..f210491 100644 --- a/R/get_recommended_interventions.R +++ b/R/get_recommended_interventions.R @@ -74,7 +74,12 @@ rank_deficient_outcome_message <- function(aliased_coef_names) { "Drop or combine the collinear predictor(s) and fit again. Common causes", "are two intervention components that are rescalings of one another, and", "fixed center and time effects whose assignments are confounded rather", - "than crossed." + "than crossed. A model can also be saturated rather than collinear, with", + "as many coefficients as it has rows to fit them from, which is what", + "center-level data with one row per center gives once center effects are", + "included: there dropping a predictor still leaves nothing to estimate the", + "rest from, and what helps is more rows per center or fitting without the", + "center effects." ) } diff --git a/R/validate_inputs.R b/R/validate_inputs.R index c4e02ae..d9cca70 100644 --- a/R/validate_inputs.R +++ b/R/validate_inputs.R @@ -1193,6 +1193,18 @@ validate_inputs <- function( #' #' @noRd refuse_invalid_center_weights <- function(center_weights_for_outcome_goal) { + # Numeric first, because the two checks below cannot see anything else. + # is.finite() is TRUE for every level of a factor, so a factor passes it, and + # the comparison that follows then gives NA and the caller fails on the + # opaque "missing value where TRUE/FALSE needed" that this function exists to + # replace. A character vector fails the comparison the same way. + if (!is.numeric(center_weights_for_outcome_goal)) { + stop(paste( + "center_weights_for_outcome_goal must be a numeric vector. A factor or", + "character vector is not weights, even when its levels look like", + "numbers: use as.numeric() on the values themselves." + )) + } if (!all(is.finite(center_weights_for_outcome_goal))) { stop(paste( "values in center_weights_for_outcome_goal must all be finite", diff --git a/tests/testthat/test-validate-inputs.R b/tests/testthat/test-validate-inputs.R index 65351d3..3fdc49a 100644 --- a/tests/testthat/test-validate-inputs.R +++ b/tests/testthat/test-validate-inputs.R @@ -346,6 +346,16 @@ test_that("center_weights_for_outcome_goal must be non-negative and finite", { } expect_error(cw(c(NA_real_, 0.5, 0.5)), "not weights") + # A factor whose levels look like numbers reaches neither of those checks: + # is.finite() is TRUE for every level, so it passes, and the comparison that + # follows gives NA and produces the same opaque base-R error this guard + # exists to replace. It is refused for what it is instead. + expect_true(all(is.finite(factor(c("0.5", "0.5"))))) + expect_error( + cw(factor(c("0.5", "0.25", "0.25"))), "must be a numeric vector" + ) + expect_error(cw(c("0.5", "0.25", "0.25")), "must be a numeric vector") + # all-zero weights are refused by the sum check, which is also what keeps the # renormalisation from dividing by zero: a vector summing to 0 is 1 away from # 1, not within the 0.001 tolerance of it. From 5f59748568b5845f8ed6dcb08e2888971ba5f8d2 Mon Sep 17 00:00:00 2001 From: Ante Bing Date: Thu, 6 Aug 2026 23:35:52 +0000 Subject: [PATCH 3/4] Pin the boundary the rank-deficient branch is chosen at The branch that reports a rank-deficient fit is entered when at least one coefficient could not be estimated, and the tests only ever passed it two. A condition of more than one aliased term therefore satisfied them while restoring the whole defect for a single aliased coefficient: grid search recommended, and more than three intervention components blamed, on a fit with two. One aliased term is not a corner case, it is what the message's own example produces, since two intervention components that are rescalings of one another alias exactly one. Both guards are pinned at that boundary now. Verified by mutating the condition in a throwaway tree and running the checks there: the assertion fails, naming the wrong message it received. --- tests/testthat/test-outcome-internals.R | 26 +++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/tests/testthat/test-outcome-internals.R b/tests/testthat/test-outcome-internals.R index 7115ccc..a5f52ed 100644 --- a/tests/testthat/test-outcome-internals.R +++ b/tests/testthat/test-outcome-internals.R @@ -1485,6 +1485,27 @@ test_that("a rank-deficient fit is not answered by recommending grid_search", { aliased <- c("period2", "period3") + # ONE aliased term takes the rank-deficient branch too, which is the boundary + # the branch is chosen at. Asserting only the two-term case leaves a test that + # a condition of "more than one aliased term" would satisfy, and a single + # aliased coefficient is what the message's own example produces: two + # intervention components that are rescalings of one another alias exactly + # one. Getting the boundary wrong restores the whole defect for that case, + # grid search recommended and more than three components blamed. + expect_error( + refuse_if_all_restarts_failed(rep(NA_real_, 11), "launch_duration"), + "rank-deficient" + ) + expect_error( + refuse_if_all_restarts_failed(rep(NA_real_, 11), "launch_duration"), + "launch_duration" + ) + expect_error( + refuse_if_all_restarts_failed(rep(NA_real_, 11), "launch_duration"), + "^(?!.*grid_search).*$", + perl = TRUE + ) + # the rank-deficient branch: the cause, the terms, and what to do about them expect_error( refuse_if_all_restarts_failed(rep(NA_real_, 11), aliased), @@ -1544,6 +1565,11 @@ test_that("a rank-deficient fit is not answered by recommending grid_search", { "period2, period3") expect_error(refuse_if_no_grid_outcome(c(NA_real_, NA_real_), aliased), "rank-deficient") + # the same boundary as above: one aliased term takes this branch too + expect_error(refuse_if_no_grid_outcome(c(0.4, NA, 0.6), "launch_duration"), + "rank-deficient") + expect_error(refuse_if_no_grid_outcome(c(0.4, NA, 0.6), "launch_duration"), + "^(?!.*grid_search).*$", perl = TRUE) expect_identical( tryCatch(refuse_if_no_grid_outcome(rep(NA_real_, 3), aliased), error = conditionMessage), From 64e754c41ebc9315c267cbb04f29d08e3c4c6191 Mon Sep 17 00:00:00 2001 From: Ante Bing Date: Thu, 6 Aug 2026 23:48:35 +0000 Subject: [PATCH 4/4] Assert the refusal names only the terms that could not be estimated Every assertion on the rank-deficient message matched a substring, so a message listing every coefficient in the model satisfied all of them. That message tells the caller to drop the intercept, both intervention components and every center, which is advice they can follow and act on, and worse than naming nothing. The terms glm() did estimate are now asserted absent. The assertions are at the call site that decides which names are passed on, because that is the only place the difference shows: passing a hand-built list straight to the guard cannot see a caller that builds the wrong list. Verified by making that caller name every coefficient in a throwaway tree, where the checks then fail. --- tests/testthat/test-outcome-internals.R | 26 +++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/tests/testthat/test-outcome-internals.R b/tests/testthat/test-outcome-internals.R index a5f52ed..587a575 100644 --- a/tests/testthat/test-outcome-internals.R +++ b/tests/testthat/test-outcome-internals.R @@ -1519,6 +1519,21 @@ test_that("a rank-deficient fit is not answered by recommending grid_search", { refuse_if_all_restarts_failed(rep(NA_real_, 11), aliased), "[Dd]rop or combine" ) + # ONLY the aliased terms. Every assertion here matches a substring, so a + # message that named every coefficient in the model would satisfy all of them + # while telling the caller to drop the intercept, both intervention components + # and every center. The estimable terms are asserted absent for that reason: + # naming a term the caller should keep is worse than naming none, since it is + # advice they can follow. + deficient_message <- tryCatch( + refuse_if_all_restarts_failed(rep(NA_real_, 11), aliased), + error = conditionMessage + ) + for (estimable in c( + "(Intercept)", "coaching_updt", "launch_duration", "center2" + )) { + expect_false(grepl(estimable, deficient_message, fixed = TRUE)) + } # and NOT the advice that cannot help. This is the assertion the shared # message fails: it recommended grid_search whatever the cause. expect_error( @@ -1687,6 +1702,17 @@ test_that("the optimizers are told which coefficients could not be estimated", { expect_match(err, "rank-deficient") expect_match(err, "period2, period3") expect_no_match(err, "grid_search") + # and ONLY those terms. Every assertion above matches a substring, so a + # message listing every coefficient in the model would satisfy them all while + # telling the caller to drop the intercept, both intervention components and + # every center. Those are the terms glm() DID estimate, asserted absent here + # because this is the call site that decides which names are passed on, and so + # the only place the difference is visible. + for (estimable in c( + "(Intercept)", "coaching_updt", "launch_duration", "center2" + )) { + expect_false(grepl(estimable, err, fixed = TRUE)) + } })