tests/testthat/test-optimization.R

In LAGOtrials: Learn-as-You-Go Adaptive Trial Optimization

# Integration/regression tests: run the full lago_optimization() pipeline and
# assert the recommended interventions match the published / previously-verified
# results. These guard the core algorithm against regressions.

test_that("BB_data binary reproduces the Nevo et al. recommendation", {
  res <- suppressWarnings(suppressMessages(lago_optimization(
    data = BB_data,
    outcome_name = "pp3_oxytocin_mother",
    outcome_type = "binary",
    glm_family = "binomial",
    intervention_components = c("coaching_updt", "launch_duration"),
    center_characteristics = "birth_volume_100",
    center_characteristics_optimization_values = 1.75,
    intervention_lower_bounds = c(1, 1),
    intervention_upper_bounds = c(40, 5),
    cost_list_of_vectors = list(c(0, 1700), c(0, 8000)),
    outcome_goal = 0.85,
    outcome_goal_intention = "maximize",
    confidence_set_grid_step_size = c(1, 0.5)
  )))
  # Nevo et al.: launch duration ~2.78 days and 1 coaching visit.
  expect_equal(res$rec_int[1], 1, tolerance = 1e-3)
  expect_equal(res$rec_int[2], 2.77847, tolerance = 1e-3)
  expect_equal(res$est_outcome_goal, 0.85, tolerance = 1e-3)
  expect_equal(res$rec_int_cost, 23927.8, tolerance = 1e-2)
  # confidence set size (regression pin from current output; ~10.5% is
  # consistent with the ~10.5% reported in Nevo et al.)
  expect_equal(res$confidence_set_size_percentage, 0.10556, tolerance = 1e-3)
  # the confidence set data.frame carries the expected columns, and the
  # recommended intervention's CI brackets the outcome goal (0.85).
  expect_true(all(
    c("CI_lower_bound", "CI_upper_bound", "cost") %in% names(res$cs)
  ))
})

test_that("the recommendation is invariant to the scale of the costs", {
  # The recommended intervention minimizes cost, so multiplying every cost by a
  # positive constant (e.g. quoting costs in thousands rather than in ones) must
  # not change it. The numerical optimizer NlcOptim::solnl() is not scale-
  # invariant on its own -- the un-normalized x1000 objective converged to a
  # different, costlier point on some platforms -- so get_recommended_
  # interventions() normalizes the objective. This pins that behaviour: the
  # recommendation and confidence set are identical at 1x and 1000x costs, and
  # only the cost is 1000x apart.
  run <- function(cost_list) {
    suppressWarnings(suppressMessages(lago_optimization(
      data = BB_data,
      outcome_name = "pp3_oxytocin_mother",
      outcome_type = "binary",
      intervention_components = c("coaching_updt", "launch_duration"),
      intervention_lower_bounds = c(1, 1),
      intervention_upper_bounds = c(40, 5),
      cost_list_of_vectors = cost_list,
      outcome_goal = 0.85,
      outcome_goal_intention = "maximize",
      confidence_set_grid_step_size = c(1, 1),
      quiet = TRUE
    )))
  }
  base <- run(list(c(0, 1.7), c(0, 8)))
  scaled <- run(list(c(0, 1700), c(0, 8000)))

  expect_equal(scaled$rec_int, base$rec_int, tolerance = 1e-6)
  expect_equal(scaled$rec_int_cost, 1000 * base$rec_int_cost, tolerance = 1e-6)
  # confidence-set membership (the component columns) is unchanged
  comps <- c("coaching_updt", "launch_duration")
  expect_equal(scaled$cs[comps], base$cs[comps])
})

test_that("mtcars continuous (identity link) is stable", {
  res <- suppressWarnings(suppressMessages(lago_optimization(
    data = mtcars,
    outcome_name = "mpg",
    outcome_type = "continuous",
    glm_family = "gaussian",
    link = "identity",
    intervention_components = c("gear", "qsec"),
    intervention_lower_bounds = c(0, 0),
    intervention_upper_bounds = c(10, 350),
    cost_list_of_vectors = list(c(0, 4), c(4, 6)),
    outcome_goal = 40,
    outcome_goal_intention = "maximize",
    confidence_set_grid_step_size = c(1, 1)
  )))
  expect_equal(res$rec_int[1], 10, tolerance = 1e-3)
  expect_equal(res$rec_int[2], 11.95743, tolerance = 1e-3)
  expect_equal(res$rec_int_cost, 115.7446, tolerance = 1e-2)
  expect_equal(res$confidence_set_size_percentage, 0.04248, tolerance = 1e-3)
})

test_that("numerical and grid_search agree on the mtcars recommendation", {
  # slow, so it runs in CI and coverage but not in the CRAN check
  skip_on_cran()
  common <- list(
    data = mtcars, outcome_name = "mpg", outcome_type = "continuous",
    glm_family = "gaussian", link = "identity",
    intervention_components = c("gear", "qsec"),
    intervention_lower_bounds = c(0, 0), intervention_upper_bounds = c(10, 350),
    cost_list_of_vectors = list(c(0, 4), c(4, 6)),
    outcome_goal = 40, outcome_goal_intention = "maximize",
    include_confidence_set = FALSE
  )
  num <- suppressWarnings(suppressMessages(do.call(lago_optimization, common)))
  grid <- suppressWarnings(suppressMessages(do.call(
    lago_optimization,
    c(common, list(
      optimization_method = "grid_search",
      optimization_grid_search_step_size = c(1, 1)
    ))
  )))
  # grid search (step 1) quantizes to integers; the numerical optimum should
  # land within one grid step of the grid solution on every component.
  expect_true(all(abs(num$rec_int - grid$rec_int) <= 1))
})

test_that("minimize direction produces a valid recommendation", {
  # exercises the double-negation minimize path end-to-end (not just the
  # error guard). Lower mpg toward 10 using disp/hp.
  res <- suppressWarnings(suppressMessages(lago_optimization(
    data = mtcars, outcome_name = "mpg", outcome_type = "continuous",
    glm_family = "gaussian", link = "identity",
    intervention_components = c("disp", "hp"),
    intervention_lower_bounds = c(0, 0), intervention_upper_bounds = c(500, 350),
    cost_list_of_vectors = list(c(0, 4), c(4, 6)),
    outcome_goal = 10, outcome_goal_intention = "minimize",
    optimization_method = "grid_search",
    optimization_grid_search_step_size = c(10, 10),
    include_confidence_set = FALSE
  )))
  # recommendation reported on the original (un-negated) scale
  expect_equal(res$rec_int, c(500, 230), tolerance = 1e-3)
  expect_equal(res$est_outcome_goal, 9.8495, tolerance = 1e-2)
})

test_that("return value has the expected shape", {
  res <- suppressWarnings(suppressMessages(lago_optimization(
    data = mtcars, outcome_name = "mpg", outcome_type = "continuous",
    glm_family = "gaussian", link = "identity",
    intervention_components = c("gear", "qsec"),
    intervention_lower_bounds = c(0, 0), intervention_upper_bounds = c(10, 350),
    cost_list_of_vectors = list(c(0, 4), c(4, 6)),
    outcome_goal = 40, outcome_goal_intention = "maximize",
    include_confidence_set = FALSE
  )))
  expect_true(all(c("rec_int", "rec_int_cost", "est_outcome_goal") %in% names(res)))
  expect_length(res$rec_int, 2)
})

test_that("the overall test result is returned (not just printed) when a group column is present", {
  # slow, so it runs in CI and coverage but not in the CRAN check
  skip_on_cran()
  bb <- BB_data
  bb$group <- ifelse(bb$pre_post == 0, "control", "treatment")
  common <- list(
    data = bb, outcome_name = "pp3_oxytocin_mother", outcome_type = "binary",
    intervention_components = c("coaching_updt", "launch_duration"),
    center_characteristics = "birth_volume_100",
    center_characteristics_optimization_values = 1.75,
    intervention_lower_bounds = c(1, 1), intervention_upper_bounds = c(40, 5),
    cost_list_of_vectors = list(c(0, 1700), c(0, 8000)),
    outcome_goal = 0.85, outcome_goal_intention = "maximize"
  )
  # returned both when the confidence set is skipped ...
  res_no_cs <- suppressWarnings(suppressMessages(do.call(
    lago_optimization, c(common, list(include_confidence_set = FALSE))
  )))
  expect_false(is.null(res_no_cs$test_results))
  expect_true(all(c("test_stat", "p_val") %in% names(res_no_cs$test_results)))

  # ... and when it is computed.
  res_cs <- suppressWarnings(suppressMessages(do.call(
    lago_optimization,
    c(common, list(
      include_confidence_set = TRUE, confidence_set_grid_step_size = c(1, 0.5)
    ))
  )))
  expect_false(is.null(res_cs$test_results))
  expect_equal(res_cs$test_results$test_stat, res_no_cs$test_results$test_stat)

  # NULL (not missing) when no valid group column is supplied.
  common_no_group <- common
  common_no_group$data <- BB_data # BB_data has no 'group' column
  res_no_group <- suppressWarnings(suppressMessages(do.call(
    lago_optimization, c(common_no_group, list(include_confidence_set = FALSE))
  )))
  expect_null(res_no_group$test_results)
})

test_that("a single intervention component does not crash the optimizer (#54)", {
  # slow, so it runs in CI and coverage but not in the CRAN check
  skip_on_cran()
  # data[, one_col] used to collapse to a vector, breaking colMeans() (shrinking
  # path) and the confidence-set construction (data[[col]] / apply over rows).
  # Uses a base data.frame (not a tibble) so the [, ] drop actually triggers.
  d <- data.frame(
    dose = rep(0:3, each = 4),
    y = c(0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1)
  )
  base <- list(
    data = d, outcome_name = "y", outcome_type = "binary",
    intervention_components = "dose",
    intervention_lower_bounds = 0, intervention_upper_bounds = 3,
    cost_list_of_vectors = list(c(0, 1)),
    outcome_goal_intention = "maximize"
  )

  # reachable goal -> confidence-set path
  res_cs <- suppressWarnings(suppressMessages(do.call(
    lago_optimization,
    c(base, list(
      outcome_goal = 0.6, include_confidence_set = TRUE,
      confidence_set_grid_step_size = 1
    ))
  )))
  expect_length(res_cs$rec_int, 1)
  expect_true("dose" %in% names(res_cs$cs))

  # unreachable goal -> shrinking path (the colMeans() site)
  res_shrink <- suppressWarnings(suppressMessages(do.call(
    lago_optimization,
    c(base, list(outcome_goal = 0.999, include_confidence_set = FALSE))
  )))
  expect_length(res_shrink$rec_int, 1)

  # continuous outcome + CS exercises get_vcov / prepare_design_matrix on
  # predictors_data, the confidence_set_processor drop = FALSE site that the
  # binary CS path above does not reach.
  set.seed(1)
  dc <- data.frame(x = seq(0, 10, length.out = 30))
  dc$y <- 1.5 * dc$x + rnorm(30, 0, 1)
  res_cts <- suppressWarnings(suppressMessages(lago_optimization(
    data = dc, outcome_name = "y", outcome_type = "continuous",
    glm_family = "gaussian", link = "identity",
    intervention_components = "x",
    intervention_lower_bounds = 0, intervention_upper_bounds = 10,
    cost_list_of_vectors = list(c(0, 1)),
    outcome_goal = 12, outcome_goal_intention = "maximize",
    include_confidence_set = TRUE, confidence_set_grid_step_size = 1
  )))
  expect_length(res_cts$rec_int, 1)
  expect_true("x" %in% names(res_cts$cs))
})

test_that("weights length is validated against the number of observations, not columns (#54)", {
  # length(data[, comp]) counted columns for a tibble / single-column drop, so
  # correctly-sized weights on a single-component tibble were wrongly rejected.
  skip_if_not_installed("tibble")
  tb <- tibble::tibble(
    dose = rep(0:3, each = 4),
    y = c(0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1)
  )
  base <- list(
    data = tb, outcome_name = "y", outcome_type = "binary",
    intervention_components = "dose",
    intervention_lower_bounds = 0, intervention_upper_bounds = 3,
    cost_list_of_vectors = list(c(0, 1)),
    outcome_goal = 0.6, outcome_goal_intention = "maximize",
    include_confidence_set = FALSE
  )
  # correctly-sized weights (one per observation) are accepted
  expect_length(
    suppressWarnings(suppressMessages(do.call(
      lago_optimization, c(base, list(weights = rep(1, nrow(tb))))
    )))$rec_int,
    1
  )
  # wrong-length weights are still rejected
  expect_error(
    suppressWarnings(suppressMessages(do.call(
      lago_optimization, c(base, list(weights = rep(1, 5)))
    ))),
    "number of observations"
  )
})


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.
  #
  # 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
  # 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)))
  aliased <- suppressWarnings(glm(
    EBP_proportions ~ center + period + coaching_updt + launch_duration,
    data = bbp, family = quasibinomial(link = "logit")
  ))
  # the precondition, asserted rather than assumed: if a future R made this fit
  # full rank the restarts would stop failing and this test would be vacuous
  expect_true(anyNA(coef(aliased)))

  run <- function() {
    suppressWarnings(suppressMessages(lago_optimization(
      data = bbp,
      outcome_name = "EBP_proportions",
      outcome_type = "continuous",
      glm_family = "quasibinomial",
      link = "logit",
      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, 1700), c(0, 8000)),
      outcome_goal = 0.85,
      include_center_effects = TRUE,
      include_time_effects = TRUE,
      time_effect_optimization_value = 1,
      include_confidence_set = FALSE,
      quiet = TRUE
    )))
  }

  # 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)

  # 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",
      outcome_type = "continuous",
      glm_family = "quasibinomial",
      link = "logit",
      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, 1700), c(0, 8000)),
      outcome_goal = 0.85,
      include_center_effects = TRUE,
      include_time_effects = TRUE,
      time_effect_optimization_value = 1,
      optimization_method = "grid_search",
      optimization_grid_search_step_size = c(4, 1),
      include_confidence_set = FALSE,
      quiet = TRUE
    )))
  }
  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)
  )

  # 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
  # what confines the failures above to the rank-deficient fixture rather than
  # to fixed center and time effects in general.
  full_rank <- bbp
  full_rank$center <- factor(rep_len(paste0("s", 1:5), nrow(full_rank)))
  expect_false(anyNA(coef(suppressWarnings(glm(
    EBP_proportions ~ center + period + coaching_updt + launch_duration,
    data = full_rank, family = quasibinomial(link = "logit")
  )))))
  ok <- suppressWarnings(suppressMessages(lago_optimization(
    data = full_rank,
    outcome_name = "EBP_proportions",
    outcome_type = "continuous",
    glm_family = "quasibinomial",
    link = "logit",
    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, 1700), c(0, 8000)),
    outcome_goal = 0.85,
    include_center_effects = TRUE,
    include_time_effects = TRUE,
    time_effect_optimization_value = 1,
    include_confidence_set = FALSE,
    quiet = TRUE
  )))
  expect_length(ok$rec_int, 2)
  expect_true(all(ok$rec_int >= c(1, 1)))
  expect_true(all(ok$rec_int <= c(40, 5)))

  # the comment above says BOTH optimizers succeed on the full-rank fit, so both
  # are run: the call above takes the default numerical method, this one the grid
  # search that fails on the aliased fixture. Asserting it is what confines that
  # failure to rank deficiency rather than to the method.
  ok_grid <- suppressWarnings(suppressMessages(lago_optimization(
    data = full_rank,
    outcome_name = "EBP_proportions",
    outcome_type = "continuous",
    glm_family = "quasibinomial",
    link = "logit",
    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, 1700), c(0, 8000)),
    outcome_goal = 0.85,
    optimization_method = "grid_search",
    optimization_grid_search_step_size = c(10, 2),
    include_center_effects = TRUE,
    include_time_effects = TRUE,
    time_effect_optimization_value = 1,
    include_confidence_set = FALSE,
    quiet = TRUE
  )))
  expect_length(ok_grid$rec_int, 2)
  expect_true(all(ok_grid$rec_int >= c(1, 1)))
  expect_true(all(ok_grid$rec_int <= c(40, 5)))
})

test_that("a degenerate all-zero cost falls back to an unnormalized objective", {
  # Exercises the cost-scale guard in get_recommended_interventions(): when the
  # cost at the upper bound is not positive (here every cost coefficient is 0, so
  # cost(x) == 0 everywhere) the normalization divisor falls back to 1 instead of
  # dividing by zero. The optimization must still run and return a feasible
  # recommendation that meets the goal, at cost 0.
  res <- suppressWarnings(suppressMessages(lago_optimization(
    data = BB_data,
    outcome_name = "pp3_oxytocin_mother",
    outcome_type = "binary",
    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, 0), c(0, 0)),
    outcome_goal = 0.85,
    outcome_goal_intention = "maximize",
    include_confidence_set = FALSE,
    quiet = TRUE
  )))
  expect_length(res$rec_int, 2)
  expect_true(all(res$rec_int >= c(1, 1)))
  expect_true(all(res$rec_int <= c(40, 5)))
  expect_equal(res$rec_int_cost, 0)
  expect_gte(res$est_outcome_goal, 0.85 - 1e-6)
})