Fix/refuse rank deficient fit

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",

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

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