diff --git a/R/get_recommended_interventions.R b/R/get_recommended_interventions.R index f210491..44cc3f1 100644 --- a/R/get_recommended_interventions.R +++ b/R/get_recommended_interventions.R @@ -64,8 +64,8 @@ unachievable_goal_message <- function(lower_outcome_goal) { #' @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)", + "No outcome can be estimated from the outcome model, because it 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", diff --git a/R/outcome_model_fitting.R b/R/outcome_model_fitting.R index e9f5082..2b3ed29 100644 --- a/R/outcome_model_fitting.R +++ b/R/outcome_model_fitting.R @@ -58,6 +58,81 @@ outcome_model_fitting <- function( )) } + # refuse a rank-deficient fit up front, but only when the aliasing lands on a + # coefficient the optimization actually reads. glm() returns NA for a + # coefficient it could not estimate -- two predictors carrying the same + # information, or a saturated fit. An NA in a coefficient the optimizer uses + # makes every predicted outcome NA and no optimization can then proceed, so + # it is fatal and raised here rather than left to fail downstream. But an NA + # in an ADDITIONAL COVARIATE is not fatal: rec_int_processor() never reads + # those into any coefficient vector the optimizer sees, so the recommendation + # is exactly the one the drop-that-covariate fit gives. Refusing on any NA + # turned that usable run into an error, so the set below is exactly the + # coefficients the optimizer will read: the intercept and intervention + # components (interaction terms included, since they are part of + # intervention_components), the fixed center effects, the fixed time effects, + # and the center characteristics. It is built with the same + # term-to-coefficient machinery rec_int_processor() uses to read them, so the + # two cannot drift: a factor's coefficient is named after its level, not its + # column, and is resolved through the model's term mapping rather than by + # matching names. The names are also the only place the aliased TERMS can + # still be named -- by the time an outcome is NA the NA has been summed into + # the center-level effects and carries their names instead. + all_coefs <- coef(model) + coef_mapping <- term_coef_names(model) + named_predictors <- claimed_coef_names(model, coef_mapping, c( + "(Intercept)", intervention_components, additional_covariates, + center_characteristics + )) + optimizer_coef_names <- c( + "(Intercept)", + intervention_components, + if (include_center_effects) { + fixed_effect_coef_names( + "center", coef_mapping, names(all_coefs), named_predictors + ) + }, + if (include_time_effects) { + fixed_effect_coef_names( + "period", coef_mapping, names(all_coefs), named_predictors + ) + }, + if (!is.null(center_characteristics)) { + unlist( + lapply(center_characteristics, predictor_coef_names, coef_mapping), + use.names = FALSE + ) + } + ) + optimizer_coefs <- all_coefs[optimizer_coef_names] + aliased_coef_names <- names(optimizer_coefs)[is.na(optimizer_coefs)] + if (length(aliased_coef_names) > 0) { + stop(rank_deficient_outcome_message(aliased_coef_names)) + } + + # an aliased ADDITIONAL COVARIATE is not fatal -- the optimizer never reads + # it, so the recommendation is unchanged from the fit without it -- but glm() + # dropping it means it contributed nothing, which the caller may not have + # intended. Warn, naming the covariate, and let the optimization continue, in + # keeping with the other non-fatal fit diagnostics below. + if (!is.null(additional_covariates)) { + dropped_covariates <- Filter(function(covariate) { + covariate_coefs <- predictor_coef_names(covariate, coef_mapping) + any(is.na(all_coefs[covariate_coefs])) + }, additional_covariates) + if (length(dropped_covariates) > 0) { + warning(paste0( + "The additional covariate(s) ", + paste(dropped_covariates, collapse = ", "), + " could not be estimated by glm() (their coefficient(s) are NA), so ", + "they were dropped from the fit and contribute nothing to the ", + "recommended intervention. This usually means they are collinear with ", + "other predictors. The recommendation is the one the fit without them ", + "gives. If that is not intended, drop or combine the covariate(s)." + )) + } + } + # run non-fatal fit diagnostics. These only warn; LAGO optimization always # continues so the user still gets a recommended intervention, but is told # when the outcome model fit is questionable and the recommendation should diff --git a/tests/testthat/test-outcome-internals.R b/tests/testthat/test-outcome-internals.R index 587a575..5ec3ab1 100644 --- a/tests/testthat/test-outcome-internals.R +++ b/tests/testthat/test-outcome-internals.R @@ -1716,6 +1716,179 @@ test_that("the optimizers are told which coefficients could not be estimated", { }) +test_that("a rank-deficient fit is refused at fitting, naming the terms", { + # A rank-deficient design makes glm() return NA for one or more coefficients. + # Those NAs used to flow into the optimizer, where every predicted outcome + # became NA and the run failed downstream with the past-tense message. The + # refusal is now made UP FRONT, at the fit, before diagnose_model_fit() and + # before any intervention is tried, so the caller is stopped as soon as the + # model cannot support an outcome rather than after a search over it fails. + # + # This is a behavioural change from db9cfc4: there outcome_model_fitting() + # RETURNED the rank-deficient model, so this whole test failed with "did not + # throw the expected error". Asserted on the fit directly, which is the site + # the refusal moved to. + omf <- getFromNamespace("outcome_model_fitting", "LAGO") + bbp <- as.data.frame(BB_proportions) + + # 1) a rescaled duplicate of an intervention component: launch_dup is exactly + # 2 * launch_duration, so the two carry the same information and glm() + # aliases the second. Exactly one coefficient is NA, which is the boundary + # the refusal fires at. + bbp$launch_dup <- bbp$launch_duration * 2 + dup_err <- tryCatch( + suppressWarnings(suppressMessages(omf( + data = bbp, + outcome_name = "EBP_proportions", + family_object = quasibinomial(link = "logit"), + intervention_components = c("launch_duration", "launch_dup"), + weights = rep(1, nrow(bbp)), + center_characteristics = NULL, + additional_covariates = NULL + ))), + error = conditionMessage + ) + expect_match(dup_err, "rank-deficient") + expect_match(dup_err, "launch_dup") + expect_match(dup_err, "[Dd]rop or combine") + # the estimable terms are not named as something to drop: naming a term the + # caller should keep is worse than naming none. launch_duration is the term + # glm() DID estimate here (launch_dup is the aliased rescaling of it), so it + # must not appear in the list of coefficients to drop. + expect_false(grepl("(Intercept)", dup_err, fixed = TRUE)) + expect_false(grepl("launch_duration", dup_err, fixed = TRUE)) + + # 2) a saturated center-level fit: one row per center with fixed center + # effects leaves no residual degrees of freedom, so glm() cannot estimate + # the intervention coefficients and aliases them. This is the saturated + # cause the message describes, distinct from the collinear one above. + cl <- data.frame( + center = factor(paste0("c", 1:6)), + coaching_updt = c(0, 5, 10, 15, 20, 25), + launch_duration = c(1, 2, 3, 4, 5, 6), + proportion = c(0.2, 0.3, 0.4, 0.5, 0.6, 0.7), + center_sample_size = rep(10, 6) + ) + sat_err <- tryCatch( + suppressWarnings(suppressMessages(omf( + data = cl, + input_data_structure = "center_level", + outcome_name = "proportion", + family_object = quasibinomial(link = "logit"), + intervention_components = c("coaching_updt", "launch_duration"), + weights = NULL, + center_characteristics = NULL, + additional_covariates = NULL, + include_center_effects = TRUE + ))), + error = conditionMessage + ) + expect_match(sat_err, "rank-deficient") + expect_match(sat_err, "coaching_updt, launch_duration") + expect_match(sat_err, "saturated") + + # 3) a full-rank fit is unaffected: it returns the model, and nothing about it + # is refused. This is what confines the refusal to rank deficiency rather + # than to fixed effects or center-level data in general. + ok <- suppressWarnings(suppressMessages(omf( + data = bbp, + outcome_name = "EBP_proportions", + family_object = quasibinomial(link = "logit"), + intervention_components = c("coaching_updt", "launch_duration"), + weights = rep(1, nrow(bbp)), + center_characteristics = NULL, + additional_covariates = NULL + ))) + expect_false(anyNA(coef(ok$model))) + expect_s3_class(ok$model, "glm") + + # 4) an aliased ADDITIONAL COVARIATE is NOT refused. The optimizer never reads + # an additional covariate into any coefficient vector it works on -- they + # appear only in the exclusion list rec_int_processor() holds back from the + # fixed effects, never as a coefficient the recommendation is computed from + # -- so glm() aliasing one does not make an outcome NA and the run still + # produces a recommendation. Refusing on ANY NA coefficient turned that + # usable run into an error, which is the regression this fixes. adj_dup is + # exactly 2 * adj, so glm() aliases it, but neither is an intervention + # component, a center characteristic or a fixed effect. + bbp$adj <- as.numeric(seq_len(nrow(bbp)) %% 7) + 0.3 + bbp$adj_dup <- bbp$adj * 2 + covariate_fit <- suppressWarnings(suppressMessages(omf( + data = bbp, + outcome_name = "EBP_proportions", + family_object = quasibinomial(link = "logit"), + intervention_components = c("coaching_updt", "launch_duration"), + weights = rep(1, nrow(bbp)), + center_characteristics = NULL, + additional_covariates = c("adj", "adj_dup") + ))) + # it returned a model rather than refusing, and glm() did alias the duplicate, + # so this is the aliased case and not a full-rank one that never exercised it. + expect_s3_class(covariate_fit$model, "glm") + expect_true(anyNA(coef(covariate_fit$model))) + expect_true(is.na(coef(covariate_fit$model)[["adj_dup"]])) + # the coefficients the optimizer DOES read are all estimable, which is why the + # fit is kept: the intercept and both intervention components are not NA. + expect_false(anyNA(coef(covariate_fit$model)[ + c("(Intercept)", "coaching_updt", "launch_duration") + ])) +}) + + +test_that("an aliased additional covariate warns but is not refused", { + # The counterpart to the refusal above: the covariate case must not stop, but + # glm() dropping the covariate is worth telling the caller about, since the + # covariate then contributes nothing to the recommendation, which the caller + # may not have intended. Precedent on the branch's parent is to warn, not + # error, for a non-fatal fit diagnostic. The warning names the dropped + # covariate and says the recommendation is the drop-that-covariate one. + omf <- getFromNamespace("outcome_model_fitting", "LAGO") + bbp <- as.data.frame(BB_proportions) + bbp$adj <- as.numeric(seq_len(nrow(bbp)) %% 7) + 0.3 + bbp$adj_dup <- bbp$adj * 2 + + expect_warning( + suppressMessages(omf( + data = bbp, + outcome_name = "EBP_proportions", + family_object = quasibinomial(link = "logit"), + intervention_components = c("coaching_updt", "launch_duration"), + weights = rep(1, nrow(bbp)), + center_characteristics = NULL, + additional_covariates = c("adj", "adj_dup") + )), + "adj_dup" + ) + # and the message says the recommendation is the drop-that-covariate one, so + # the caller knows the run is usable and what its result stands for. + expect_warning( + suppressMessages(omf( + data = bbp, + outcome_name = "EBP_proportions", + family_object = quasibinomial(link = "logit"), + intervention_components = c("coaching_updt", "launch_duration"), + weights = rep(1, nrow(bbp)), + center_characteristics = NULL, + additional_covariates = c("adj", "adj_dup") + )), + "fit without them" + ) + + # a full-rank additional covariate does not warn at all, so the diagnostic is + # confined to the case where glm() actually dropped one. adj alone is + # estimable (adj_dup was its aliased rescaling), so it produces no warning. + expect_silent(suppressMessages(omf( + data = bbp, + outcome_name = "EBP_proportions", + family_object = quasibinomial(link = "logit"), + intervention_components = c("coaching_updt", "launch_duration"), + weights = rep(1, nrow(bbp)), + center_characteristics = NULL, + additional_covariates = "adj" + ))) +}) + + 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