R/validate_inputs.R

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

validate_inputs <- function(
    data,
    input_data_structure = "individual_level",
    outcome_name,
    outcome_type,
    intervention_components,
    intervention_lower_bounds,
    intervention_upper_bounds,
    outcome_goal,
    outcome_goal_intention,
    power_goal,
    power_goal_approach,
    num_centers_in_next_stage = NULL,
    patients_per_center_in_next_stage = NULL,
    icc = NULL,
    power_goal_cluster_id = NULL,
    unit_costs = NULL,
    unit_cost_doses = NULL,
    default_cost_fxn_type = "cubic",
    cost_list_of_vectors = NULL,
    glm_family = "default",
    link = "default",
    optimization_method = "numerical",
    confidence_set_alpha = 0.05,
    weights = NULL,
    center_characteristics = NULL,
    center_characteristics_optimization_values = NULL,
    main_components = NULL,
    time_effect_optimization_value = NULL,
    additional_covariates = NULL,
    center_weights_for_outcome_goal = NULL,
    optimization_grid_search_step_size = NULL,
    confidence_set_grid_step_size = NULL,
    include_confidence_set = TRUE,
    include_center_effects = NULL,
    center_effects_optimization_values = NULL,
    include_time_effects = FALSE,
    include_interaction_terms = FALSE,
    prev_recommended_interventions = NULL,
    shrinkage_threshold = 0.25,
    quiet = FALSE) {
  # check if the input data is null
  if (is.null(data)) {
    stop("The argument 'data' is NULL.")
  }
  # check if the input data is a data frame type
  if (!is.data.frame(data)) {
    stop("The argument 'data' must be a data frame.")
  }
  # check if the input data is empty
  if (nrow(data) == 0 || ncol(data) == 0) {
    stop("The argument 'data' is empty.")
  }

  # check if the input data structure is a character type
  if (!is.character(input_data_structure)) {
    stop("The input_data_structure must be a character string.")
  }
  # check if the input data structure is supported
  if (!(input_data_structure %in% c("individual_level", "center_level"))) {
    stop(paste(
      "\ninput_data_structure:", input_data_structure,
      "\ninput_data_structure must be either 'individual_level'",
      "or 'center_level'."
    ))
  }

  # check if the outcome name is a character type
  if (!is.character(outcome_name)) {
    stop("The outcome name must be a character string.")
  }
  # check if the outcome name is empty
  if (nchar(outcome_name) == 0) {
    stop("The outcome name is empty.")
  }
  # check if the outcome name is one of the columns in the data frame
  if (!(outcome_name %in% names(data))) {
    stop("The outcome name must be presented in the provided input data.")
  }
  # an outcome with no observed value leaves nothing to fit or to compare a
  # goal against (its mean with na.rm = TRUE would be NaN)
  if (all(is.na(data[[outcome_name]]))) {
    stop(paste0(
      "The outcome column '", outcome_name, "' has no observed values."
    ))
  }

  # check if the outcome type is a character type
  if (!is.character(outcome_type)) {
    stop("The outcome type must be a character string.")
  }
  # check if the outcome type is either continuous or binary
  allowed_outcome_types <- c("continuous", "binary")
  if (!(outcome_type %in% allowed_outcome_types)) {
    stop(paste(
      "Outcome type", outcome_type,
      "must be either 'continuous' or 'binary'."
    ))
  }
  # a factor or character outcome has no mean to compare the goal against
  outcome_column <- data[[outcome_name]]
  if (!is.numeric(outcome_column) && !is.logical(outcome_column)) {
    stop(paste0(
      "The outcome column '", outcome_name, "' must be numeric",
      if (outcome_type == "binary") ", with 0 and 1 for a binary outcome",
      ", but it is a ", class(outcome_column)[1], " column. Recode it, for ",
      "example with ",
      if (outcome_type == "binary") {
        paste0("as.numeric(data$", outcome_name, " == \"yes\") for yes and no values, or ")
      },
      "as.numeric(as.character(data$", outcome_name, ")) for values written as numbers."
    ))
  }
  # check if the input data structure is center_level when
  # the outcome type is continuous.
  if (input_data_structure == "center_level") {
    if (outcome_type == "continuous") {
      stop(paste(
        "For continuous outcomes, LAGO requires individual",
        "level data for the optimization step."
      ))
    }
  }

  # check if the input data has all the required columns
  required_columns <- c()
  if (input_data_structure == "center_level") {
    required_columns <- c(
      required_columns,
      "proportion",
      "center_sample_size"
    )
  }

  missing_columns <- setdiff(required_columns, names(data))
  if (length(missing_columns) > 0) {
    stop(sprintf(
      "The following required columns are missing: %s",
      paste(missing_columns, collapse = ", ")
    ))
  }
  # center level data is fitted on its proportion column, so the outcome must be it
  if (input_data_structure == "center_level" && !identical(outcome_name, "proportion")) {
    stop(paste0(
      "For center level data the outcome is the 'proportion' column, so ",
      "outcome_name must be \"proportion\", not \"", outcome_name, "\"."
    ))
  }

  # check if the glm family is a character type
  if (!is.character(glm_family)) {
    stop("The glm family must be a character string.")
  }
  # assign default glm_family based on the outcome type
  if (glm_family == "default") {
    if (outcome_type == "continuous") {
      glm_family <- "gaussian"
    } else if (outcome_type == "binary") {
      glm_family <- "binomial"
    }
  }

  # check if the link option is a character type
  if (!is.character(link)) {
    stop("The link option must be a character string.")
  }
  # check if the provided link option is supported. The set is
  # supported_outcome_links(), the one place the links the outcome machinery
  # implements are written down, plus "default", which is resolved to one of
  # them from the glm family just below.
  # "probit" and "log" used to be accepted here and are not any more. Nothing
  # implemented their inverse links, so get_outcome() fell through to its
  # identity branch and reported the LINEAR PREDICTOR as the estimated
  # probability or mean: a probit fit asked to hold the outcome at most 0.20
  # reported 0.200 as reached while the actual probability was pnorm(0.20) =
  # 0.579, and get_confidence_set() built its interval on the wrong scale too.
  # Accepting a link only some of the machinery handles reports wrong numbers
  # with no sign anything is wrong, so the link is refused up front instead.
  supported_link_options <- c(supported_outcome_links(), "default")
  if (!(link %in% supported_link_options)) {
    stop(unsupported_link_message(link))
  }
  # assign link based on the glm_family
  if (link == "default") {
    # a quasibinomial outcome is a proportion, so it takes the logit link too
    if (glm_family %in% c("binomial", "quasibinomial")) {
      link <- "logit"
    } else if (glm_family == "gaussian") {
      link <- "identity"
    }
  }

  # center level rows are weighted by their center_sample_size
  if (!is.null(weights) && input_data_structure == "center_level") {
    input_warning(paste(
      "'weights' is ignored for center level data, whose rows are weighted",
      "by their center_sample_size."
    ))
  }
  # check if the provided weights option is numeric
  if (!is.null(weights)) {
    if (!is.numeric(weights)) {
      stop("The 'weights' option must be a numeric vector.")
    }
    # check if the provided weights option has the same length as the number of
    # observations. Use nrow(data) rather than length(data[, comp]): the latter
    # returns the observation count for a base data.frame but the column count
    # for a tibble or a single-column drop, so it mis-validates weights.
    if (length(weights) != nrow(data)) {
      stop(paste0(
        "The length of the weights must be the same as the ",
        "number of observations in the provided input data."
      ))
    }
  }

  # check if intervention_components is a character vector type
  if (!(is.vector(intervention_components) &&
    is.character(intervention_components))) {
    stop("Interventions list must be a character vector.")
  }
  # check if intervention_components are all columns in the data frame
  if (!all(intervention_components %in% names(data))) {
    stop(paste0(
      "All elements in intervention_components ",
      "must be columns in the data frame."
    ))
  }
  # refuse any intervention component whose column is entirely NA. glm()'s
  # internal na.omit drops every row in that case, and the fit dies with an
  # opaque error ("Argument mu must be a nonempty numeric vector") that never
  # names the component. This runs after the column-existence check above so
  # data[[component]] is always a real column, and before the interaction-term
  # backticks are added below so the names still match the data columns.
  # all(is.na(col)) treats numeric, factor and character columns alike, so an
  # all-NA component of any type is caught here rather than left to fail
  # cryptically inside model fitting.
  all_na_components <- intervention_components[vapply(
    intervention_components,
    function(component) all(is.na(data[[component]])),
    logical(1)
  )]
  if (length(all_na_components) > 0) {
    stop(paste0(
      "The intervention component(s) ",
      paste(all_na_components, collapse = ", "),
      " are entirely NA, so glm() would drop every row and they ",
      "cannot enter the model. Remove them from intervention_components ",
      "or supply observed values."
    ))
  }
  # a component is a dose, and glm() would code a non-numeric one as levels
  non_numeric_components <- intervention_components[!vapply(
    intervention_components,
    function(component) is.numeric(data[[component]]),
    logical(1)
  )]
  if (length(non_numeric_components) > 0) {
    stop(paste0(
      "The intervention component(s) ",
      paste0("'", non_numeric_components, "'", collapse = ", "),
      " must be numeric columns. Recode a yes or no component as 0 and 1, ",
      "for example with as.numeric()."
    ))
  }

  # check if center_characteristics is a character vector type
  if (length(center_characteristics) > 0) {
    if (!(is.vector(center_characteristics) &&
      is.character(center_characteristics))) {
      stop("center_characteristics must be a character vector.")
    }
    # check if center_characteristics are all columns in the data frame
    if (!all(center_characteristics %in% names(data))) {
      stop(paste0(
        "All elements in center_characteristics ",
        "must be columns in the data frame."
      ))
    }
    # check if center_characteristics_optimization_values
    # is a numeric vector type
    if (!(is.vector(center_characteristics_optimization_values) &&
      is.numeric(center_characteristics_optimization_values))) {
      stop(paste0(
        "center_characteristics_optimization_values ",
        "must be a numeric vector."
      ))
    }
    # check if center_characteristics_optimization_values is not null
    if (is.null(center_characteristics_optimization_values)) {
      stop(paste0(
        "center_characteristics_optimization_values is NULL. ",
        "You decided to include center characteristics in ",
        "the model, please either provide values of the ",
        "center characteristics for LAGO optimization, ",
        "or consider dropping the center characteristics."
      ))
    }
    # check if length of center_characteristics_optimization_values is the same
    # as the length of center_characteristics
    if (length(center_characteristics) !=
      length(center_characteristics_optimization_values)) {
      stop(paste(
        "The length of center_characteristics must be the",
        "same as the length of",
        "center_characteristics_optimization_values."
      ))
    }
  }

  # if include_center_effects is not specified, set the default value: TRUE
  # for center level data, FALSE otherwise. A value the user gave is kept.
  if (is.null(include_center_effects)) {
    include_center_effects <- input_data_structure == "center_level"
  }
  # check whether the include_center_effects indicator is boolean
  if (!is.logical(include_center_effects) || length(include_center_effects) != 1 ||
    is.na(include_center_effects)) {
    stop("The include_center_effects indicator must be a boolean.")
  }

  # check whether center_characteristics and center effects are both included
  if (include_center_effects && !is.null(center_characteristics)) {
    stop(paste(
      "Fixed center effects and center characteristics cannot be",
      "both included in the outcome model. Please",
      "pick one to use."
    ))
  }

  # check if the input data has the "center" column if
  # include_center_effects is TRUE
  if (include_center_effects) {
    if (!"center" %in% names(data)) {
      stop(paste(
        "To include fixed center effects in the outcome model",
        "please make sure that the input data has a 'center'",
        "column to identify facilities."
      ))
    }
    # refused here, as the center weights below would otherwise fail with an unclear message
    if (all(is.na(data$center))) {
      stop(paste0(
        "Fixed center effects need at least two centers, but the data has no ",
        "observed center. Fill in the center column, or set ",
        "include_center_effects = FALSE."
      ))
    }
    # if "center" column of the input data is not a factor,
    # change it to factor and print a warning message
    if (!is.factor(data$center)) {
      data$center <- as.factor(data$center)
      input_note(paste(
        "'center' column is not a factor type. To ensure",
        "the correct model fit, it has been converted",
        "to the factor type."
      ))
    }
    # an unordered factor, so it is treated as center labels and not as a scale
    if (is.ordered(data$center)) {
      class(data$center) <- "factor"
      input_note(paste(
        "'center' column is an ordered factor. It has been converted to an",
        "unordered factor with the same levels, so the model fits one effect",
        "per center."
      ))
    }
    # one level per center present, as glm() drops unused levels (align handles lost ones)
    data$center <- droplevels(data$center)
  }

  # preliminary check for center_weights_for_outcome_goal
  if (!is.null(center_weights_for_outcome_goal)) {
    # check if include center effects is set to FALSE
    if (!include_center_effects) {
      stop(paste(
        "'include_center_effects' must be set to TRUE,",
        "for values of 'center_weights_for_outcome_goal'",
        "to be used."
      ))
    }
    # check if center_weights_for_outcome_goal is a numeric vector
    if (!is.numeric(center_weights_for_outcome_goal)) {
      stop(paste("center_weights_for_outcome_goal must be a numeric vector."))
    }
    # check if center_weights_for_outcome_goal has the same length as the
    # number of facilities provided in the input data
    # the centers present (a missing center id is not a facility)
    n_facilities <- nlevels(data$center)
    if (length(center_weights_for_outcome_goal) != n_facilities) {
      stop(paste(
        "The length of center_weights_for_outcome_goal must match",
        "the number of unique facilities in the input data."
      ))
    }
    if (!is.null(center_effects_optimization_values)) {
      input_warning(paste(
        "Both 'center_weights_for_outcome_goal' and",
        "'center_effects_optimization_values' are given. The center weights",
        "are used and the named center is ignored. Give only one of them."
      ))
    }
  } else {
    # if center_weights_for_outcome_goal is not specified, but
    # include_center_effects is set to TRUE, we calculate its default values
    if (include_center_effects) {
      # if not specifing center_effects_optimization_values,
      # then we calculate the default values for center weights
      if (is.null(center_effects_optimization_values)) {
        if (input_data_structure == "center_level") {
          # each center's size is its patients, the sum over its rows, as glm() weights them
          center_sizes <- tapply(
            data$center_sample_size,
            data$center,
            function(x) if (all(is.na(x))) NA_real_ else sum(x, na.rm = TRUE)
          )
          no_size <- names(center_sizes)[is.na(center_sizes)]
          if (length(no_size) > 0) {
            stop(paste0(
              "center_sample_size is missing for every row of the center(s) ",
              paste0("'", no_size, "'", collapse = ", "), ", so they have no ",
              "size to weight them by. Fill in their sample size, or drop them."
            ))
          }
          total_sample_size <- sum(center_sizes)
        } else if (input_data_structure == "individual_level") {
          total_sample_size <- sum(as.numeric(table(data$center)))
          center_sizes <- as.numeric(table(data$center))
        }
        center_weights_for_outcome_goal <- center_sizes / total_sample_size
      } else {
        # if center_effects_optimization_values is specified,
        # we set the weight for that specific center to 1, rest to 0
        if (!is.character(center_effects_optimization_values)) {
          stop(paste(
            "center_effects_optimization_values must be a character string",
            "naming one center."
          ))
        }
        if (length(center_effects_optimization_values) != 1) {
          stop("center_effects_optimization_values must be a single center name.")
        }
        # the levels, so a missing center without an explicit NA level is not a center
        if (!center_effects_optimization_values %in% levels(data$center)) {
          missing_name <- is.na(center_effects_optimization_values)
          stop(paste0(
            "center_effects_optimization_values (",
            if (missing_name) "NA_character_" else center_effects_optimization_values,
            ") must be one of the centers in the 'center' column.",
            if (missing_name) {
              paste(
                " To select the rows with a missing center, make them an",
                "explicit level with addNA()."
              )
            }
          ))
        }
        # by level, so an explicit NA level keeps its place among the centers
        center_weights_for_outcome_goal <- as.numeric(
          levels(data$center) %in% center_effects_optimization_values
        )
      }
    } else {
      center_weights_for_outcome_goal <- 1
    }
  }

  # check the values of center_weights_for_outcome_goal. Here rather than in the
  # branch above that reads them from the caller, so that the checks cover every
  # way the vector can be arrived at: the caller's weights, the center-size
  # default and the single-named-center indicator alike. The center-size default
  # 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)) {
    # 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)
    # and whether the vector is a set of weights at all as a SET, i.e. sums to
    # 1. Shared with the exported get_confidence_set() for the same reason the
    # checks above are: see refuse_non_unit_weight_sum().
    refuse_non_unit_weight_sum(center_weights_for_outcome_goal)
    weights_sum <- sum(center_weights_for_outcome_goal)
    # The tolerance above says the input was MEANT to be a set of weights. It
    # does not make it one: the weights multiply the per-center outcomes and are
    # then summed, so a set summing to 0.999 scales every reported outcome by
    # 0.999 and one summing to 1.001 by 1.001. On the documented example weights
    # that is a bias of 6.7e-4 on the reported outcome, in the direction of the
    # error and on every outcome the run reports, including the goal comparison
    # the recommendation is chosen against.
    #
    # Renormalising here rather than tightening the tolerance: the tolerance is
    # documented and callers rely on it, so tightening it would turn accepted
    # input into a hard error. Renormalising leaves weights that sum to exactly
    # 1 bit for bit unchanged, since dividing by 1 is exact, and removes the
    # bias for everyone else. It is not quite a no-op on the defaults computed
    # below: a vector of center sizes divided by its own total need not sum to
    # exactly 1, and about one such vector in four hundred sums to one unit in
    # the last place away from it, so those runs move by the same amount. That
    # is the correction doing its job on a genuinely non-unit sum, not an error
    # introduced here.
    #
    # Done once, here, where the weights are validated, so nothing downstream
    # needs to know: this is where lago_optimization() obtains them, whether
    # from the caller, from the center sample sizes, or from a single named
    # center. After the fit, align_center_weights_to_fit() keeps only the
    # centers glm() fitted (renormalizing when it drops one), and that vector
    # is what every optimizer, the shrinking method and the confidence set get.
    #
    # Silently, deliberately. The correction is at most 0.1% of a weight and is
    # what the caller already asked for by passing something the tolerance
    # accepts as a set of weights. Warning on it would fire on rounded input
    # that the documentation invites -- c(0.333, 0.333, 0.334) is the natural
    # way to write three equal weights -- so it would be noise on correct
    # usage while the actual behaviour is now right. A caller who did not intend
    # weights is still refused, by the check above.
    center_weights_for_outcome_goal <-
      center_weights_for_outcome_goal / weights_sum
  }

  # check if include_time_effects is logical
  if (!is.logical(include_time_effects)) {
    stop(paste("'include_time_effects' must be a boolean."))
  }
  family_object <- switch(glm_family,
    "binomial" = binomial(link = link),
    "gaussian" = gaussian(link = link),
    "quasibinomial" = quasibinomial(link = link),
    stop(paste("Unsupported glm_family:", glm_family))
  )


  # if include_time_effects is set to TRUE
  if (include_time_effects) {
    # the input data must have a 'period' column
    if (!"period" %in% names(data)) {
      stop(paste(
        "To include fixed time effects in the outcome model,",
        "please make sure that the input data has a 'period'",
        "column to identify time period."
      ))
    }
    # the time_effect_optimization_value must be provided, one numeric number
    if (is.null(time_effect_optimization_value) ||
      !is.numeric(time_effect_optimization_value) ||
      length(time_effect_optimization_value) != 1) {
      stop(paste(
        "'include_time_effects' is set to TRUE,",
        "'time_effect_optimization_value' must be provided,",
        "and it must be a numeric value."
      ))
    }
    # if "period" column of the input data is not a factor,
    # change it to factor and print a warning message
    if (!is.factor(data$period)) {
      data$period <- as.factor(data$period)
      input_note(paste(
        "'period' column is not a factor type. To ensure",
        "the correct model fit, it has been converted",
        "to the factor type."
      ))
    }
    # an unordered factor, so it is treated as period labels and not as a scale
    if (is.ordered(data$period)) {
      class(data$period) <- "factor"
      input_note(paste(
        "'period' column is an ordered factor. It has been converted to an",
        "unordered factor with the same levels, so the model fits one effect",
        "per period."
      ))
    }
  }

  # if include_interaction_terms is set to TRUE
  if (include_interaction_terms) {
    # check if interaction terms are included in intervention_components
    # and follow the correct naming scheme
    if (!any(grepl(":", intervention_components, fixed = TRUE))) {
      stop(paste(
        "Please make sure the interaction terms are included",
        "as part of the 'intervention_components', and the interaction terms",
        "need to have the format of 'X1:X5' "
      ))
    }

    # if so, check if main_components is null
    if (is.null(main_components)) {
      stop(paste(
        "'main_components' must be defined as a character vector",
        "if interaction terms are",
        "included in 'intervention_components'."
      ))
    }
    # check if main_components is a character vector
    if (!is.character(main_components)) {
      stop(paste("'main_components' must be a character vector."))
    }

    # if not null, check if main_components have the correct naming scheme
    if (any(grepl(":", main_components, fixed = TRUE))) {
      stop(paste(
        "Interaction terms should be included as part of the",
        "'intervention_components', not part of the",
        "'main_components'."
      ))
    }

    # check if main_components are all columns in the data frame
    # need this to calculate the step size if the step size is not specified.
    if (optimization_method == "grid_search" &&
      is.null(optimization_grid_search_step_size)) {
      if (!all(main_components %in% names(data))) {
        stop(paste(
          "When 'optimization_method' is set to 'grid_search',",
          "and optimization_grid_search_step_size is not specified,",
          "all elements in the main_components",
          "must be columns in the data frame for auto calculating",
          "the step size."
        ))
      }
    }

    # check if interaction terms can be found in the main_components
    check_term <- function(term) {
      if (grepl(":", term, fixed = TRUE)) {
        parts <- unlist(strsplit(term, ":", fixed = TRUE))
        return(all(parts %in% main_components))
      } else {
        return(term %in% main_components)
      }
    }
    interaction_check_results <- sapply(intervention_components, check_term)
    if (!all(interaction_check_results)) {
      stop(paste(
        "Each component of 'intervention_components' including",
        "interaction terms (separated by ':')",
        "must be present in 'main_components'."
      ))
    }

    # add backticks to interactions terms in intervention components
    # for glm model fitting
    needs_backticks <- grepl(":", intervention_components)
    intervention_components[needs_backticks] <- paste0(
      "`", intervention_components[needs_backticks], "`"
    )
  } else {
    # if interaction terms are not included, then the intervention
    # components shouldn't have ":" in their names
    if (any(grepl(":", intervention_components, fixed = TRUE))) {
      stop(paste(
        "':' is found in the intervention_components,",
        "did you mean to set include_interaction_terms to TRUE?"
      ))
    }
  }

  # check if additional_covariates is a character vector
  if (length(additional_covariates) > 0) {
    if (!(is.vector(additional_covariates) &&
      is.character(additional_covariates))) {
      stop("additional_covariates must be a character vector.")
    }
    # check if additional_covariates are all columns in the data frame
    if (!all(additional_covariates %in% names(data))) {
      stop(paste0(
        "All elements in additional_covariates ",
        "must be columns in the data frame."
      ))
    }
    # refuse any additional covariate whose column is entirely NA. glm()'s
    # internal na.omit drops every row in that case, and the fit dies with an
    # opaque error that never names the column. This runs after the
    # column-existence check above so data[[covariate]] is always a real
    # column. all(is.na(col)) treats numeric, factor and character columns
    # alike, so an all-NA covariate of any type is caught here rather than
    # left to fail cryptically inside model fitting.
    all_na_covariates <- additional_covariates[vapply(
      additional_covariates,
      function(covariate) all(is.na(data[[covariate]])),
      logical(1)
    )]
    if (length(all_na_covariates) > 0) {
      stop(paste0(
        "The additional covariate(s) ",
        paste(all_na_covariates, collapse = ", "),
        " are entirely NA, so glm() would drop every row and they ",
        "cannot enter the model. Remove them from additional_covariates ",
        "or supply observed values."
      ))
    }
  }

  # the model formula is built from these names, so each must be a valid R
  # name. An interaction term such as a:b is quoted where it is used.
  model_names <- unique(c(
    outcome_name, intervention_components, main_components,
    additional_covariates, center_characteristics
  ))
  model_names <- model_names[!grepl(":", model_names, fixed = TRUE)]
  bad_names <- model_names[make.names(model_names) != model_names]
  if (length(bad_names) > 0) {
    stop(paste0(
      "The column name(s) ", paste0("'", bad_names, "'", collapse = ", "),
      " are not valid R names, so they cannot be used in the outcome model ",
      "formula. Rename them, for example with make.names(), which gives ",
      paste0("'", make.names(bad_names), "'", collapse = ", "), "."
    ))
  }

  # check if intervention_lower_bounds and intervention_upper_bounds are both
  # numerical vectors
  if (!(is.vector(intervention_lower_bounds) &&
    is.numeric(intervention_lower_bounds))) {
    stop("intervention_lower_bounds must be a numeric vector.")
  }
  if (!(is.vector(intervention_upper_bounds) &&
    is.numeric(intervention_upper_bounds))) {
    stop("intervention_upper_bounds must be a numeric vector.")
  }

  # check if the dimension of intervention_lower_bounds matches
  # the dimension of intervention_components or the dimension
  # of main_components
  if (
    (!include_interaction_terms &&
      (length(intervention_lower_bounds) != length(intervention_components))) ||
      (include_interaction_terms &&
        (length(intervention_lower_bounds) != length(main_components)))
  ) {
    stop(paste(
      "Without interaction terms,",
      "the lengths of 'intervention_lower_bounds' and",
      "'intervention_components' must be the same.",
      "With interaction terms,",
      "the lengths of 'intervention_lower_bounds' and",
      "'main_components' must be the same."
    ))
  }

  # check if lower bounds list and upper bounds list have the same length,
  # and if the lower bounds are <= upper bounds respectively
  if (length(intervention_lower_bounds) != length(intervention_upper_bounds)) {
    stop("The lengths of lower and upper bounds must be the same.")
  }
  if (any(intervention_lower_bounds < 0)) {
    stop("The intervention must have non-negative values only.")
  }
  invalid_indices <- which(intervention_upper_bounds <
    intervention_lower_bounds)
  if (length(invalid_indices) > 0) {
    stop(paste(
      "Invalid bounds at position(s):",
      paste(invalid_indices, collapse = ", "),
      paste0(
        "\nUpper bounds for the interventions must be greater ",
        "than or equal to the lower bounds."
      )
    ))
  }

  # check if lower bounds and upper bounds make sense
  # by checking the actual min and max values of components in the data
  if (include_interaction_terms) {
    all_components <- main_components
  } else {
    all_components <- intervention_components
  }
  for (i in seq_along(all_components)) {
    min_val <- min(data[[all_components[i]]], na.rm = TRUE)
    max_val <- max(data[[all_components[i]]], na.rm = TRUE)
    if (intervention_lower_bounds[i] > min_val) {
      warning(paste(
        "The lower bound for the intervention component",
        all_components[i],
        "is greater than the minimum value in the data."
      ))
    }
    if (intervention_upper_bounds[i] < max_val) {
      warning(paste(
        "The upper bound for the intervention component",
        all_components[i],
        "is less than the maximum value in the data."
      ))
    }
  }

  # check if unit_costs is valid
  if (!is.null(unit_costs)) {
    # check if unit_costs is a numeric vector
    if (!is.numeric(unit_costs)) {
      stop("unit_costs must be a numeric vector.")
    }
    if (any(!is.finite(unit_costs)) || any(unit_costs < 0)) {
      stop("unit_costs must be finite and not negative.")
    }
    # check if the dimension of unit_costs matches
    # the dimension of intervention_components or the dimension
    # of main_components
    if (
      (!include_interaction_terms &&
        (length(unit_costs) != length(intervention_components))) ||
        (include_interaction_terms &&
          (length(unit_costs) != length(main_components)))
    ) {
      stop(paste(
        "Without interaction terms,",
        "the lengths of 'unit_costs' and",
        "'intervention_components' must be the same.",
        "With interaction terms,",
        "the lengths of 'unit_costs' and",
        "'main_components' must be the same."
      ))
    }
    # check if default_cost_fxn_type is valid
    if (!is.character(default_cost_fxn_type)) {
      stop("default_cost_fxn_type must be a character string.")
    }
    # check if default_cost_fxn_type is one of the supported types
    supported_cost_fxn_types <- c("linear", "cubic")
    if (!(default_cost_fxn_type %in% supported_cost_fxn_types)) {
      stop(sprintf(
        "default_cost_fxn_type must be one of the supported types. The supported types are: %s.",
        paste(supported_cost_fxn_types, collapse = ", ")
      ))
    }
    # the dose each unit cost was measured at pins the cubic cost there
    if (!is.null(unit_cost_doses)) {
      if (!is.numeric(unit_cost_doses) || length(unit_cost_doses) != length(unit_costs) ||
        any(!is.finite(unit_cost_doses))) {
        stop(paste(
          "unit_cost_doses must be a finite numeric vector with one dose for",
          "each unit cost."
        ))
      }
      if (any(unit_cost_doses < intervention_lower_bounds) ||
        any(unit_cost_doses > intervention_upper_bounds)) {
        stop(paste(
          "Each value of unit_cost_doses must lie within the lower and upper",
          "bounds of its intervention component."
        ))
      }
    } else if (is.null(cost_list_of_vectors) && default_cost_fxn_type == "cubic") {
      input_warning(paste(
        "'unit_cost_doses' is not given, so the cubic cost only matches each",
        "unit cost on average over the bounds, not at the dose it was measured",
        "at. Give the dose each unit cost was measured at to pin the cost there."
      ))
    }
  }

  # check if cost_list_of_vectors is valid
  if (!is.null(cost_list_of_vectors)) {
    # check if cost_list_of_vectors is a list,
    # and each sublist is a numeric vector
    if (!is.list(cost_list_of_vectors)) {
      stop("cost_list_of_vectors must be a list.")
    }
    # Check if all elements of sublists in cost_list_of_vectors are numeric
    all_numeric <- all(sapply(cost_list_of_vectors, function(sublist) {
      all(sapply(sublist, is.numeric))
    }))
    if (!all_numeric) {
      stop(paste(
        "All elements in the sublists of cost_list_of_vectors",
        "must be numeric."
      ))
    }
    # check if the dimension of cost_list_of_vectors matches
    # the dimension of intervention_components or the dimension
    # of main_components
    if (
      (!include_interaction_terms &&
        (length(cost_list_of_vectors) != length(intervention_components))) ||
        (include_interaction_terms &&
          (length(cost_list_of_vectors) != length(main_components)))
    ) {
      stop(paste(
        "Without interaction terms,",
        "the lengths of 'cost_list_of_vectors' and",
        "'intervention_components' must be the same.",
        "With interaction terms,",
        "the lengths of 'cost_list_of_vectors' and",
        "'main_components' must be the same."
      ))
    }
    # check if default_cost_fxn_type is also provided
    if (!is.null(default_cost_fxn_type) && !quiet) {
      input_note(paste(
        "When 'cost_list_of_vectors' is provided,",
        "'default_cost_fxn_type' is ignored."
      ))
    }
  }

  # check if unit_costs and cost_list_of_vectors are both null
  if (is.null(unit_costs) && is.null(cost_list_of_vectors)) {
    stop(paste(
      "Both 'unit_costs' and 'cost_list_of_vectors' are NULL.",
      "Please provide at least one of them."
    ))
  }

  # check if the outcome goal optimization method is a character type
  if (!is.character(optimization_method)) {
    stop("The optimization_method must be a character string.")
  }
  # check if the outcome goal optimization method is one of the defined methods
  optimization_method_list <- c("numerical", "grid_search")
  if (!(optimization_method %in% optimization_method_list)) {
    stop(paste(
      "optimization_method must be one of the supported methods.",
      "The supported methods are:", optimization_method_list, "."
    ))
  }

  # if optimization_grid_search_step_size is provided, validate it before it
  # can influence optimization_method below, so that an invalid step size is
  # rejected up front rather than after a "switched to grid_search" message.
  if (!is.null(optimization_grid_search_step_size)) {
    # it needs to be a numeric vector
    if (!(is.vector(optimization_grid_search_step_size) &&
      is.numeric(optimization_grid_search_step_size))) {
      stop("optimization_grid_search_step_size must be a numeric vector.")
    }
    # it needs to have the same length as the number of intervention components
    if ((!include_interaction_terms &&
      (length(optimization_grid_search_step_size) !=
        length(intervention_components))) ||
      (include_interaction_terms &&
        (length(optimization_grid_search_step_size) !=
          length(main_components)))
    ) {
      stop(paste(
        "Without interaction terms,",
        "The number of step sizes provided for the grid search",
        "algorithm must be the same as the number of intervention",
        "components. Please check the length of",
        "'optimization_grid_search_step_size'.",
        "With interaction terms,",
        "The number of step sizes provided for the grid search",
        "algorithm must be the same as the length of 'main_components'."
      ))
    }
    # each step size must be finite and strictly positive, otherwise the grid
    # search seq() call would error or produce a degenerate grid.
    if (!all(is.finite(optimization_grid_search_step_size)) ||
      any(optimization_grid_search_step_size <= 0)) {
      stop(paste(
        "All elements of 'optimization_grid_search_step_size' must be",
        "finite and greater than 0."
      ))
    }
  }

  # optimization_grid_search_step_size is only used by the grid search
  # method. If it is provided while the method is "numerical", the user
  # clearly intends to use grid search, so switch the method (with a
  # message) instead of silently ignoring the step size.
  if (!is.null(optimization_grid_search_step_size) &&
    optimization_method == "numerical") {
    optimization_method <- "grid_search"
    input_note(paste(
      "'optimization_grid_search_step_size' was provided, so",
      "'optimization_method' has been switched to 'grid_search'.",
      "To use the numerical method instead, remove",
      "'optimization_grid_search_step_size'."
    ))
  }

  # when the optimization_method is grid_search,
  # check the number of intervention components
  if (optimization_method == "grid_search") {
    if (length(intervention_components) > 3) {
      input_note(paste0(
        "There are more than 3 intervention components, and the grid search ",
        "algorithm may take significant amount of time to run. ",
        "Please consider adjusting the step size, or switch to the ",
        "numerical optimization method."
      ))
    }
  }

  # at least one of outcome_goal or power_goal must be provided. LAGO
  # optimization supports three modes: outcome goal alone, power goal alone,
  # or both together (in which case the effective goal is the higher of the
  # two, see get_recommended_interventions()).
  if (is.null(outcome_goal) && is.null(power_goal)) {
    stop(paste(
      "Both 'outcome_goal' and 'power_goal' are NULL.",
      "Please provide at least one of them."
    ))
  }

  # check if the outcome goal is a numeric number (only when provided)
  if (!is.null(outcome_goal) && !is.numeric(outcome_goal)) {
    stop("The outcome goal must be a numeric value.")
  }

  # check if the outcome_goal_intention is a character type
  if (!is.character(outcome_goal_intention)) {
    stop("The outcome_goal_intention must be a character string.")
  }
  # check if the outcome_goal_intention is one of the defined intentions
  allowed_outcome_goal_intentions <- c("minimize", "maximize")
  if (!(outcome_goal_intention %in% allowed_outcome_goal_intentions)) {
    stop(paste(
      "The outcome_goal_intention must be one",
      "of the supported intentions.",
      "The supported intentions are:",
      allowed_outcome_goal_intentions, "."
    ))
  }

  # set the lower_outcome_goal flag accordingly
  # check the average outcome and print warning messages.
  # when only a power goal is provided (outcome_goal is NULL), the goal is
  # always to increase the outcome to reach the desired power, so
  # lower_outcome_goal is FALSE and the intention checks are skipped.
  # the observed mean the goal is compared with, weighted as glm() weights the
  # rows, so center level data is averaged over patients rather than rows
  observed_mean <- local({
    y <- data[[outcome_name]]
    w <- if (input_data_structure == "center_level") {
      data$center_sample_size
    } else if (!is.null(weights)) {
      weights
    } else {
      rep(1, length(y))
    }
    ok <- !is.na(y) & !is.na(w)
    m <- if (any(ok) && sum(w[ok]) > 0) stats::weighted.mean(y[ok], w[ok]) else NaN
    if (is.finite(m)) m else mean(y, na.rm = TRUE)
  })
  if (is.null(outcome_goal)) {
    lower_outcome_goal <- FALSE
  } else if (outcome_goal_intention == "minimize") {
    lower_outcome_goal <- TRUE
    if (outcome_goal > observed_mean) {
      warning(paste(
        "The specified outcome goal intention is",
        "to minimize the outcome.",
        "The observed mean of the outcome is",
        observed_mean,
        "and the specified outcome goal is", outcome_goal, "."
      ))
    }
  } else {
    lower_outcome_goal <- FALSE
    if (outcome_goal <= observed_mean) {
      warning(paste(
        "The specified outcome goal intention",
        "is to maximize the outcome.",
        "The observed mean of the outcome is",
        observed_mean,
        "and the specified outcome goal is", outcome_goal, "."
      ))
    }
  }

  # check whether the include_confidence_set indicator is boolean
  if (!is.logical(include_confidence_set)) {
    stop("The include_confidence_set indicator must be a boolean.")
  }

  # check if confidence set step size is provided
  if (include_confidence_set && !is.null(confidence_set_grid_step_size)) {
    # if confidence_set_grid_step_size is provided,
    # it needs to be a numeric vector
    if (!(is.vector(confidence_set_grid_step_size) &&
      is.numeric(confidence_set_grid_step_size))) {
      stop("confidence_set_grid_step_size must be a numeric vector.")
    }
    # if confidence_set_grid_step_size is provided,
    # it needs to have the same length as the number of intervention components
    if ((!include_interaction_terms &&
      (length(confidence_set_grid_step_size) !=
        length(intervention_components))) ||
      (include_interaction_terms &&
        (length(confidence_set_grid_step_size) !=
          length(main_components)))
    ) {
      stop(paste(
        "Without interaction terms,",
        "the number of step sizes provided for the grid",
        "search algorithm in the confidence set calculation",
        "must be the same as the length of 'intervention_components'.",
        "Please check the length of 'confidence_set_grid_step_size'.",
        "With interaction terms,",
        "the number of step sizes provided for the grid",
        "search algorithm in the confidence set calculation",
        "must be the same as the length of 'main_components'."
      ))
    }
    if (any(!is.finite(confidence_set_grid_step_size)) ||
      any(confidence_set_grid_step_size <= 0)) {
      stop("confidence_set_grid_step_size must be finite and positive.")
    }
  }

  # check if confidence_set_alpha is a numeric type
  if (!is.numeric(confidence_set_alpha)) {
    stop("confidence_set_alpha must be a numeric value.")
  }
  # check if confidence_set_alpha is between 0 and 1
  if (confidence_set_alpha <= 0 || confidence_set_alpha >= 1) {
    stop("confidence_set_alpha must be between 0 and 1.")
  }

  # check if the specified family matches the outcome type
  valid_families <- list(
    "binary" = c("binomial"),
    "continuous" = c("gaussian", "quasibinomial")
  )
  if (!glm_family %in% valid_families[[outcome_type]]) {
    stop(paste("The specified family '", glm_family,
      "' is not valid for the outcome type '", outcome_type, "'.",
      " Please select a compatible family.",
      sep = ""
    ))
  }

  # check if the prev_recommended_interventions is numeric and
  # has the same length as the number of intervention components
  if (!is.null(prev_recommended_interventions)) {
    if (!is.numeric(prev_recommended_interventions)) {
      stop("prev_recommended_interventions must be a numeric vector.")
    }
    if (include_interaction_terms) {
      all_components <- main_components
    } else {
      all_components <- intervention_components
    }
    if (length(prev_recommended_interventions) != length(all_components)) {
      stop(paste(
        "The length of prev_recommended_interventions",
        "must be the same as",
        "the number of intervention components."
      ))
    }

    # pre_recommended_interventions must be within the bounds
    for (i in 1:length(all_components)) {
      if (prev_recommended_interventions[i] <
        intervention_lower_bounds[i] ||
        prev_recommended_interventions[i] >
          intervention_upper_bounds[i]) {
        stop(paste(
          "The pre_recommended_interventions[", i, "]",
          "is outside the bounds of the intervention component",
          all_components[i], "."
        ))
      }
    }
  }

  # check if shrinkage_threshold is a numeric value
  # and is between 0 and 1
  if (!is.numeric(shrinkage_threshold)) {
    stop("shrinkage_threshold must be a numeric value.")
  }
  if (shrinkage_threshold < 0 || shrinkage_threshold > 1) {
    stop("shrinkage_threshold must be between 0 and 1.")
  }

  # check if power_goal is a numeric number if not null
  if (!is.null(power_goal)) {
    if (!is.numeric(power_goal)) {
      stop("The power goal must be a numeric value.")
    }
    # check if power_goal is between 0 and 1
    if (power_goal <= 0 || power_goal >= 1) {
      stop("The power goal must be between 0 and 1.")
    }
    # check if the outcome_type is binary
    if (outcome_type != "binary") {
      stop("For now, the power goal only works with binary outcomes.")
    }
    # a power goal requires increasing the outcome to reach the desired
    # power, so it is incompatible with minimizing the outcome. This holds
    # whether or not an outcome goal is also provided.
    if (outcome_goal_intention == "minimize") {
      stop(paste(
        "A power goal cannot be combined with",
        "outcome_goal_intention = 'minimize', because achieving a power",
        "goal requires increasing the outcome. Please set",
        "outcome_goal_intention = 'maximize', or remove the power goal."
      ))
    }
    # check if there is a 'group' column in the data
    # and it is either "treatment" or "control"
    if (!"group" %in% names(data)) {
      stop(paste(
        "To include a power goal, please make sure that",
        "there is a 'group' column in the input data."
      ))
    }
    if (!all(data$group %in% c("treatment", "control"))) {
      stop(paste(
        "The 'group' column must have values 'treatment'",
        " or 'control'."
      ))
    }
    # the power calculation uses stage-1 rows with an observed outcome, so
    # each arm needs at least one (a missing outcome is not counted)
    observed_group <- data$group[!is.na(data[[outcome_name]])]
    for (grp in c("control", "treatment")) {
      if (!any(observed_group == grp)) {
        stop(paste0(
          "A power goal needs stage-1 data in both arms, but the '", grp,
          "' arm of the 'group' column has no rows with an observed outcome."
        ))
      }
    }
    # num_centers_in_next_stage and patients_per_center_in_next_stage
    # cannot be null for power goal to work
    if (is.null(num_centers_in_next_stage) ||
      is.null(patients_per_center_in_next_stage)) {
      stop(paste(
        "num_centers_in_next_stage and",
        "patients_per_center_in_next_stage must be provided",
        "for the power goal to work."
      ))
    }
    # if num_center_in_next_stage is not null, it must be a numeric value
    if (!is.null(num_centers_in_next_stage)) {
      if (!is.numeric(num_centers_in_next_stage)) {
        stop("num_centers_in_next_stage must be a numeric value.")
      }
      if (num_centers_in_next_stage <= 0) {
        stop("num_centers_in_next_stage must be greater than 0.")
      }
    }
    # if patients_per_center_in_next_stage is not null, it must be a numeric value
    if (!is.null(patients_per_center_in_next_stage)) {
      if (!is.numeric(patients_per_center_in_next_stage)) {
        stop("patients_per_center_in_next_stage must be a numeric value.")
      }
      if (patients_per_center_in_next_stage <= 0) {
        stop("patients_per_center_in_next_stage must be greater than 0.")
      }
    }
    # icc (issue #29): optional intra-cluster correlation used to inflate the
    # power-calculation variance by a design effect. Scalar (shared across arms)
    # or length-2 c(control, treatment). Each value must be in [0, 1).
    if (!is.null(icc)) {
      if (!is.numeric(icc)) {
        stop("icc must be a numeric value (or a length-2 numeric vector).")
      }
      if (length(icc) != 1 && length(icc) != 2) {
        stop("icc must have length 1 (shared) or 2 (control, treatment).")
      }
      if (any(icc < 0) || any(icc >= 1)) {
        stop("icc must be in [0, 1).")
      }
      # a non-zero icc needs a stage-1 cluster id to compute the stage-1 design
      # effect; validate the column exists and is non-degenerate per arm.
      if (any(icc > 0)) {
        if (is.null(power_goal_cluster_id)) {
          stop(paste(
            "A non-zero 'icc' requires 'power_goal_cluster_id', the name of a",
            "column identifying stage-1 centers, so the stage-1 design effect",
            "can be computed."
          ))
        }
        if (!is.character(power_goal_cluster_id) ||
          length(power_goal_cluster_id) != 1) {
          stop("power_goal_cluster_id must be a single column name.")
        }
        if (!power_goal_cluster_id %in% names(data)) {
          stop(paste0(
            "power_goal_cluster_id column '", power_goal_cluster_id,
            "' was not found in the data."
          ))
        }
        # need at least two centers per arm to identify between-cluster
        # variance, counted on the rows the power calculation uses (an observed
        # outcome) and ignoring a missing center id, as its table() does
        for (grp in c("control", "treatment")) {
          ids <- data[
            data$group == grp & !is.na(data[[outcome_name]]),
            power_goal_cluster_id,
            drop = TRUE
          ]
          n_centers <- length(unique(ids[!is.na(ids)]))
          if (n_centers < 2) {
            stop(paste0(
              "The '", grp, "' arm has fewer than two distinct ",
              "power_goal_cluster_id centers with an observed outcome, so ",
              "the stage-1 design effect is not identifiable. Provide data ",
              "with at least two such centers per arm, or set icc = NULL."
            ))
          }
        }
      }
    }
  } else if (!is.null(icc)) {
    # icc only affects the power calculation; without a power goal it is a no-op.
    input_note("icc is ignored because no power_goal was provided.")
  }

  # check if power_goal_approach is a character type
  # and is one of the defined approaches
  if (!is.character(power_goal_approach)) {
    stop("The power_goal_approach must be a character string.")
  }
  allowed_power_goal_approaches <- c("unconditional", "conditional")
  if (!(power_goal_approach %in% allowed_power_goal_approaches)) {
    stop(paste(
      "The power_goal_approach must be one",
      "of the supported approaches.",
      "The supported approaches are:",
      allowed_power_goal_approaches, "."
    ))
  }


  list(
    data = data,
    input_data_structure = input_data_structure,
    outcome_name = outcome_name,
    outcome_type = outcome_type,
    intervention_components = intervention_components,
    intervention_lower_bounds = intervention_lower_bounds,
    intervention_upper_bounds = intervention_upper_bounds,
    outcome_goal = outcome_goal,
    unit_costs = unit_costs,
    unit_cost_doses = unit_cost_doses,
    default_cost_fxn_type = default_cost_fxn_type,
    cost_list_of_vectors = cost_list_of_vectors,
    glm_family = glm_family,
    link = link,
    optimization_method = optimization_method,
    confidence_set_alpha = confidence_set_alpha,
    weights = weights,
    center_characteristics = center_characteristics,
    center_characteristics_optimization_values = center_characteristics_optimization_values,
    main_components = main_components,
    time_effect_optimization_value = time_effect_optimization_value,
    additional_covariates = additional_covariates,
    center_weights_for_outcome_goal = center_weights_for_outcome_goal,
    optimization_grid_search_step_size = optimization_grid_search_step_size,
    confidence_set_grid_step_size = confidence_set_grid_step_size,
    include_confidence_set = include_confidence_set,
    include_center_effects = include_center_effects,
    include_time_effects = include_time_effects,
    include_interaction_terms = include_interaction_terms,
    family_object = family_object,
    link = link,
    lower_outcome_goal = lower_outcome_goal,
    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) {
  # 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",
      "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)
}


#' refuse_non_unit_weight_sum
#'
#' @description Internal guard for the center weights as a SET: refuses a
#' vector whose sum is not 1, within the tolerance validate_inputs() has always
#' used.
#'
#' @details This is the other way to the wrong number
#' refuse_invalid_center_weights() refuses. That guard rules out a vector no
#' single element of which is a weight; this one rules out a vector every
#' element of which could be, but which is not a set of weights taken together.
#' The reported outcome is sum(weight_i * outcome_i) over the center-level
#' effects, so it is a weighted mean only when the weights sum to 1: a set
#' summing to 12 scales the intervention and center-characteristic contribution
#' of every reported outcome by 12, and 16 weights of 12/16 gave a
#' confidence interval of 1.406 to 2.147 for an outcome that is a proportion.
#' Every weight there is non-negative and finite, so neither of the checks in
#' refuse_invalid_center_weights() applies, and the symptom is the same one it
#' exists to prevent.
#'
#' Shared with the exported get_confidence_set() rather than left at
#' validate_inputs(), which is the same reason the checks in
#' refuse_invalid_center_weights() are shared: get_confidence_set() does not go
#' through validate_inputs(), so a caller reaching it directly had no sum check
#' at all, and the @param text saying the weights must sum to 1 was true of one
#' entry point and not the other.
#'
#' REFUSED rather than renormalised at get_confidence_set(), while
#' validate_inputs() both refuses this and renormalises what it accepts. The two
#' are not inconsistent: both entry points refuse exactly the same vectors, in
#' the same words. What differs is what happens to a vector inside the
#' tolerance, and there validate_inputs() can do more because it OWNS the
#' weights -- it is where they are derived and normalised, and lago_optimization()
#' only narrows them afterwards to the centers glm() fitted, through
#' align_center_weights_to_fit(), which renormalises when it drops one.
#' get_confidence_set() is handed weights an optimization has
#' already run with, and is documented as computing the interval at those
#' weights, so renormalising them there would report an interval for a
#' different weighting than the point estimate it is printed beside. That is the
#' same class of defect as reporting the interval of one period beside the
#' estimate of another.
#'
#' Renormalising there would also not be a no-op on correct input, which is what
#' settles it. validate_inputs() stores center_sizes / total_sample_size, and
#' about one such vector in six hundred does not sum to exactly 1 in floating
#' point; dividing by that sum again perturbs every weight by up to one unit in
#' the last place and moves the interval. So renormalising in the name of
#' matching the other entry point would change the numbers of runs that are
#' already correct, and refusing changes none of them: every internal path
#' arrives with a sum of 1 to within rounding, because validate_inputs() has
#' already normalised it and align_center_weights_to_fit() renormalises it when
#' it drops a center.
#'
#' The tolerance is validate_inputs()' own 0.001 rather than a tighter one.
#' It is documented, callers rely on it, and it is what says the input was MEANT
#' to be a set of weights; tightening it here would turn input the package has
#' always accepted into a hard error, and the residual scaling it admits is at
#' most 0.1%, which is the tolerance's own documented meaning. What this closes
#' is the unbounded case.
#'
#' A vector summing to 0, including one that is all zeros, is refused here: it
#' is 1 away from 1. That is also what keeps validate_inputs()'
#' renormalisation from dividing by zero.
#'
#' @param center_weights_for_outcome_goal A numeric vector of center weights,
#' already known to be finite and non-negative.
#'
#' @return Invisibly NULL when the weights sum to 1 within the tolerance.
#' Raises otherwise.
#'
#' @noRd
refuse_non_unit_weight_sum <- function(center_weights_for_outcome_goal) {
  if (abs(sum(center_weights_for_outcome_goal) - 1) >= 0.001) {
    stop(paste(
      "values in center_weights_for_outcome_goal must",
      "sum up to 1."
    ))
  }
  invisible(NULL)
}

# Signal an input note, a message about how the inputs were read. A sweep shows
# each one once rather than once per run. See show_each_input_note_once().
input_note <- function(text) {
  message(structure(
    class = c("lago_input_note", "message", "condition"),
    list(message = paste0(text, "\n"), call = NULL)
  ))
}

# Signal an input warning, a warning about how the inputs were read. A sweep
# shows each one once rather than once per run, as it does an input note.
input_warning <- function(text) {
  warning(structure(
    class = c("lago_input_warning", "warning", "condition"),
    list(message = text, call = NULL)
  ))
}